Two-subsystem thermodynamics for the mechanics of aging amorphous solids

The effect of physical aging on the mechanics of amorphous solids as well as mechanical rejuvenation is modeled with nonequilibrium thermodynamics, using the concept of two thermal subsystems, namely a kinetic one and a configurational one. Earlier work (Semkiv and Hütter in J Non-Equilib Thermodyn 41(2):79–88, 2016) is extended to account for a fully general coupling of the two thermal subsystems. This coupling gives rise to hypoelastic-type contributions in the expression for the Cauchy stress tensor, that reduces to the more common hyperelastic case for sufficiently long aging. The general model, particularly the reversible and irreversible couplings between the thermal subsystems, is compared in detail with models in the literature (Boyce et al. in Mech Mater 7:15–33, 1988; Buckley et al. in J Mech Phys Solids 52:2355–2377, 2004; Klompen et al. in Macromolecules 38:6997–7008, 2005; Kamrin and Bouchbinder in J Mech Phys Solids 73:269–288 2014; Xiao and Nguyen in J Mech Phys Solids 82:62–81, 2015). It is found that only for the case of Kamrin and Bouchbinder (J Mech Phys Solids 73:269–288, 2014) there is a nontrivial coupling between the thermal subsystems in the reversible dynamics, for which the Jacobi identity is automatically satisfied. Moreover, in their work as well as in Boyce et al. (Mech Mater 7:15–33, 1988), viscoplastic deformation is driven by the deviatoric part of the Cauchy stress tensor, while for Buckley et al. (J Mech Phys Solids 52:2355–2377, 2004) and Xiao and Nguyen (J Mech Phys Solids 82:62–81, 2015) this is not the case.

processing, most prominently on the cooling rate. The difference of amorphous to (semi)crystalline solids is that amorphous solids typically show ongoing changes over extended periods of time. This is because the amorphous solids, being rapidly quenched to temperatures below their glass transition temperature, are manifestly in a nonequilibrium glassy state. While the kinetics in that state is extraordinarily slow, the system still evolves toward lower energy states. With this dynamics, called aging, changes in physical properties of the material can be observed over time. 1 A comprehensive overview over the main physical effects due to aging dynamics in the case of polymers is given in [1,2]. In particular, the effect of physical aging on the mechanical behavior of glassy polymers has been studied in [3][4][5][6][7]. The effects of physical aging can (to some extent) be reversed by the application of a significant nonelastic deformation, which coined the term of mechanical rejuvenation [8].
The formal description of aging dynamics often is based on a concept of internal state variables [3], particularly in the form of configurational temperature, first proposed by Tool [9], and then used intensively, e.g., in [10][11][12][13][14][15][16]. According to this concept, the amorphous solid is described by two, rather than one, thermal degrees of freedom, namely a kinetic one and a configurational one. Generally, the character of these two thermal degrees of freedom is either that of temperature, entropy density, or internal energy density, respectively. While the kinetic (vibrational) degree of freedom accounts for intra-basin thermodynamics, the configurational one describes the inter-basin thermodynamics [17][18][19][20][21][22][23][24]. Models have been developed for the evolution of the kinetic and configurational thermal degrees of freedom [10][11][12][13][14][15], as well as constitutive expressions have been formulated for a stress tensor and for a plastic flow rule for the mechanical deformation of aging solids. However, in these models the mutual interaction of the kinetic and configurational subsystems is in general incomplete, as discussed in the sequel of this paper, and leaves room for further investigation and generalization.
In a recent publication [16], a closed-form dynamic model for the elasto-viscoplastic deformation of aging amorphous solids has been developed, according to nonequilibrium thermodynamic principles and making use of kinetic and configurational entropy densities. Next to evolution equations for the thermal and mechanical degrees of freedom, constitutive relations for the stress tensor as well as for the plastic flow rule have been formulated. As elaborated in detail in [16], the split of the energy dissipated by the viscoplastic deformation into the kinetic and configurational thermal subsystems is nonunique from a thermodynamics perspective. It has also been observed that the choice for that split significantly affects the driving force for plastic flow. All this concerns the irreversible dynamics of the model. With respect to the reversible dynamics, no coupling of the kinetic and configurational subsystems, namely the respective entropy densities, has been considered so far.
The main goal of this paper is, in view of these lessons learned from the irreversible dynamics, to systematically explore the possibilities for letting the kinetic and configurational subsystems interact nontrivially also through the reversible contributions to the dynamics and to study the ramifications of doing so.
The paper is organized as follows. After presenting the nonequilibrium thermodynamics framework (Sect. 2), the model for the elasto-viscoplastic behavior of aging amorphous solids is developed (Sect. 3), specifically including the nontrivial reversible coupling between the kinetic and configurational degrees of freedom. This general model is then compared in detail with other approaches in the literature, i.e., [3,6,10,14,15] (Sect. 4), before drawing conclusions (Sect. 5).
Before starting with the main content, some general comments are made about the notation used in this paper. Greek indices are used for Cartesian components, and Einstein's summation convention is used for indices that occur twice. Furthermore, with respect to operators, subscripts α and (α, β) imply contraction with any vector A α and any tensor A αβ multiplied from the left, respectively, while subscripts γ and (γ , ε) imply contraction with the vector A γ and tensor A γ ε multiplied from the right, respectively. With respect to partial derivatives, the quantity f ,X i stands for the partial derivative of a function f with respect to the variable X i . Unless stated otherwise, it is tacitly assumed that this derivative is taken at constant X \ X i , with X = (m, s, η, F e ). Finally, A d αβ denotes the deviatoric part of the tensor A αβ , i.e., A d αβ ≡ A αβ − (A μμ /3)δ αβ .

GENERIC framework
The framework used in this paper to formulate the dynamics of the system is the general equation for the nonequilibrium reversible-irreversible coupling (GENERIC) [25][26][27], which has also been used in our preceding work [16]. It is a formal thermodynamic procedure to set up evolution equations for a set of variables X . As discussed in [16], the reversible ("rev") and irreversible ("irr") contributions are clearly distinguished, with where ∂ t denotes the partial derivative with respect to time t. The reversible dynamics is driven by the gradient of energy E, mediated by the Poisson operator L, while the irreversible dynamics is driven by the gradient of entropy S, mediated by the friction matrix M. As part of the GENERIC structure, on the one hand, the reversible dynamics must not change the entropy, and the Poisson operator must be skew-symmetric and satisfy the Jacobi identity, i.e., with the Poisson bracket defined by {A, B} = δ A/δX , L δ B/δX based on the inner product ·, · , and functionals A, B, and C. With these conditions, the Hamiltonian structure of the reversible dynamics is guaranteed, and specifically the time-structure invariance is respected by way of the Jacobi identity [27]. On the other hand, the irreversible dynamics must not change the energy, and the friction matrix must be symmetric and be positive semi-definite, i.e., The resulting irreversible dynamics is thus a generalization of the Ginzburg-Landau equation. The reversible (2) and irreversible (3) contributions are coupled by the degeneracy conditions (4) and (7). For completeness, it is mentioned that M generally must be Onsager-Casimir symmetric, but for the application in this paper M is assumed to be (Onsager) symmetric. For more details, the reader is referred to [25][26][27].

Set of variables, and generating functionals
The first step in the model formulation concerns the specification of meaningful dynamic variables, for which evolution equations are sought. To that end, we depart from our earlier results [16], where the elasto-viscoplastic behavior of amorphous materials was studied, with a focus on the effects of physical aging and mechanical rejuvenation. In that work, the set of variables in an Eulerian (i.e., spatial) setting consisted of the momentum density m, the elastic part of the deformation gradient, F e , and the entropy densities (per unit volume) of the kinetic-and configurational subsystems, s and η, respectively. In summary, the complete set of dynamic variables, which shall also be used in this paper, is thus given by X = {m, s, η, F e }.
It should be noted that the case of a single entropy density is equivalent to what has been studied in [28,29] for describing finite-deformation nonisothermal elasto-viscoplasticity. However, in order to account for the effect of physical aging, also called structural relaxation, a split of that thermal degree of freedom into its kinetic and configurational contributions is considered, in analogy to [12,14,15,20,23]. While from a practical (interpretation) point of view, one might be tempted to use two temperatures instead of two (partial) entropies, this would result in unnecessary technical complications, which would rather disguise the essence of the two-subsystem approach. Therefore, it has been decided on purpose in this paper to consider the kinetic and configurational subsystem entropies as fundamental variables.
In view of the set of variables X , it is convenient to write the energy and entropy in the form the total entropy density [19]. Specifically, the energy is split into its macroscopic kinetic energy and internal energy parts. The mass density is given by ρ = ρ 0 / det F e , with ρ 0 the mass density in the undeformed state, in which it is tacitly assumed that all volumetric change is elastic in origin (see further below), i.e., det F e = det F with F the total deformation gradient. Based on (10)- (11), the functional derivatives become with velocity v γ = m γ /ρ, and with ρ ,F e γ ε = −ρ F e,−1 εγ . Furthermore, we have introduced the kinetic (T ) and configurational (θ ) temperatures, respectively, according to which have also been used frequently in the literature for modeling the two-subsystem thermodynamics [9][10][11][12][13][14][15]19]. From a macroscopic perspective, the two temperatures defined in (15) and (16) denote the variables conjugate to the respective subsystem entropies, which can be used to transition from the internal energy density as a thermodynamic potential to the Helmholtz free energy density, by Legendre transformation. However, in a statistical-mechanics setting, one can think of these temperatures as being the (inverse of the) Lagrange parameters related to the respective subsystem energies in a generalized canonical formulation (see [30]). Such statistical-mechanics background is also required to determine the two temperatures concretely. For example, in the classical case of a single temperature only, the equi-partition theorem can be employed to quantify the temperature, i.e., the absolute temperature being directly related to the kinetic energy of the microscopic constituents. This suggests that, for the case of two subsystems discussed in this paper, the kinetic temperature can be addressed in a similar way. As far as the configurational temperature is concerned, one can examine the dependence of the thermodynamic potential on the configurational-subsystem variable. This has been done by computational means [20], by determining the dependence of the configurational entropy on the energy of inherent states [31]. The authors are not aware of any publication that addresses the measurement of the two temperatures with experimental techniques.

Reversible dynamics
In order to specify the reversible dynamics (2), a choice needs to be made for the Poisson operator L. To that end, we depart from our earlier results [16], where the evolution equation for X has been formulated. In that earlier model, there is no reversible exchange between the subsystem entropy densities s and η, since they evolve simply according to ∂ t s| rev = −∇ γ (sv γ ) and ∂ t η| rev = −∇ γ (ηv γ ). However, one can imagine that in general there is a coupling between the kinetic and configurational degrees of freedom during an imposed (macroscopically affine) deformation. In the language of the potential energy landscape, this would mean that when deforming the system, changes in the intra-basin and the inter-basin distributions would result that are nontrivially coupled. It is thus of interest to examine the possibility of a coupling between the subsystem entropy densities upon (reversible) deformation. While the evolution equation for the total entropy density, ∂ t s t | rev = −∇ γ (s t v γ ), will be unaltered, specifically the possibility of an exchange term between s and η is sought.
Since the coupling between s and η relates to the macroscopic velocity field v and in view of (13), it seems natural to incorporate an exchange between s and η via the first column of L. It is thus suggested that the Poisson operator in [16] is modified in its first column and, due to the anti-symmetry condition (5), also in the first row, namely [16] The Y -related contributions in the first column represent the exchange between the subsystem entropies, since the reversible contributions to their evolution (2) with Poisson operator (17) and energy gradient (13) becomes The exchange term is proportional to the velocity gradient, i.e., homogeneous (bulk) translations do not lead to any entropy exchange, as required. Therefore, this serves as a justification of the structure of the Y -related contributions in the first column of (17). Due to the anti-symmetry condition (5), the Y -related contributions in the first row of (17) can be derived. For this Poisson operator, also the evolution equations for the momentum density and for the elastic part of the deformation gradient can be determined according to (2) with Poisson operator (17) and energy gradient (13), which is as expected. However, it should be noted that the Cauchy stress tensor contains a contribution related to the exchange tensor Y , with the abbreviationf = f /ρ for any quantity f . This implies that the Cauchy stress tensor carries a signature of the entropy exchange; however, (22) suggests that this contribution vanishes as the system tends to thermal equilibrium, T θ . This issue will be discussed below in more detail. If the system is not in thermal equilibrium, the stress tensor is not simply related to the derivative of the thermodynamic potential with respect to deformation (i.e., hyperelastic type), but rather it contains also another contribution, akin to hypoelasticity [32]. As far as the conditions (4)-(6) are concerned, it is clear from the above that the anti-symmetry condition (5) is satisfied by construction. Furthermore, it can be shown readily that the degeneracy condition (4) is fulfilled. One can thus concentrate on the Jacobi identity (6), which puts a constraint on the tensor Y αβ , namely For technical convenience, in all partial derivatives in relation to the Jacobi calculation here,ŝ t is used as a thermal variable rather thanŝ, together withη. To proceed, it is noted that the requirement of a symmetric stress tensor implies that also Y αβ is symmetric. Therefore, one can writê with a symmetric tensorŶ αβ . Note that this factorization (24) is without loss of generality, as long as F e αβ is invertible. Inserting this factorization into (23), one obtains While this condition (25) holds in full generality, we now proceed to discuss it by making two assumptions, for the purpose of illustration. The first assumption is thatŶ αβ depends on F e μν only through the right Cauchy-Green strain tensor C e γ δ = F e μγ F e μδ . Making use of the general relation ∂g/∂ F e αβ = 2 F e αγ (∂g/∂C e γβ ) for any function g of the right Cauchy-Green strain tensor, the condition (25) turns into The second assumption concerns the tensorial character ofŶ αβ . Particularly, it is assumed thatŶ αβ is of quasipotential form, i.e.,Ŷ αβ = (1/y 1 )(∂ y 2 /∂C e αβ )|ŝ t ,η with y 1 and y 2 two scalar-valued functions of the variableŝ s t ,η, and C e . Doing so, the condition (26) becomes This condition can be processed further if the functions y 1 and y 2 depend on C e only through the three invariants J 1 = trC e , J 2 = ln det C e , and J 3 = −trC e,−1 of C e , that have the useful property ∂ J k /∂ C e = C e,1−k [27,33]. In this case, the nontrivial part of condition (27) can be cast into the form It is straightforward to show that (28) is indeed a nontrivial condition. For example for y 1 ≡ 1, one concludes from (28) that only if y 2 depends only on a single invariant J k (e.g., on det C e , by way of ρ = ρ 0 / √ det C e , which implies Y αβ ∝ δ αβ ), the Jacobi identity is trivially satisfied. In contrast, if y 2 depends on more than one of the invariants of C e , the Jacobi identity puts a constraint on the function y 2 . It is mentioned that another example where the Jacobi identity imposes nontrivial constraints on the dynamics has been discussed in the context of complex fluids [27,33].
This elaboration of the implications of the Jacobi identity on the choice of the tensor Y is concluded by a comment about the stress tensor. Using the representation (24), ifŶ is of quasi-potential form, and ifê t depends on F e only through C e , the stress tensor (22) can be written as It is important to note that (θ − T )/y 1 is in general a function of C e ,ŝ, andη, also by virtue of the definitions (15)- (16). Therefore, the expression in the parenthesis in (29) cannot always be written as the C e γ δ -derivative of a potential.

Irreversible dynamics
The irreversible dynamics, due to (3), involves the friction matrix M and is driven by the entropy gradient δS/δX . While the previous section explained in detail how to account for entropy exchange due to reversible dynamics, it is pointed out that entropy exchange due to irreversible dynamics has been accounted for already in [16]. For this reason, only the main results about the irreversible dynamics are summarized here as far as they are relevant for this paper, while the reader interested in the full details is referred to [16].
Four different irreversible processes are accounted for in this model. These are specifically: heat conduction in each of two thermal subsystems (kinetic and configurational), heat exchange between the subsystems, and viscoplastic deformation. Correspondingly, for each of these four processes, there are thermodynamic fluxes, namely J s and J η for the bulk subsystem heat fluxes, J for the subsystem exchange heat flux, and J F e for the viscoplastic deformation, related to the irreversible dynamics of F e . It has been shown earlier that these thermodynamic fluxes result in irreversible dynamics of X , according to (3), of the form [16] which satisfies the degeneracy condition (7). The parameter ψ defines how the dissipation caused by viscoplastic deformation is distributed among the kinetic and configurational entropies. The tensor P γ εαβ is a projection defined by P γ εαβ = δ αγ δ βε − (1/3)F e,−1 εγ F e αβ to ensure that the viscoplastic change to F e is indeed isochoric [29].
The evolution equations (30)-(33) highlight the physics of thermodynamic fluxes, for which the constitutive expressions are required in order to complete the model formulation. Specifying suitable constitutive relations for these thermodynamic fluxes (see [16] for details), the irreversible contributions to the evolution equation can be written in the form The thermal conductivities for the two subsystems are denoted by λ s and λ η , both of which must be symmetric and positive semi-definite to ensure the conditions (8) and (9). Similarly, the quantity μ x , which describes the heat transfer between the two subsystems, must be nonnegative. The quantity D p,d αβ is the deviatoric part of driven by the stress-like quantity The fact that d γ δ is deviatoric implies that also Z γ δ must be a traceless tensor, which is tacitly used throughout this paper. In order to comply with the condition for symmetry (8), the tensor αβγ δ must be symmetric, i.e., αβγ δ = γ δαβ , and furthermore, we impose αβγ δ = βαγ δ = βαδγ ; the positive semi-definiteness condition (9) is guaranteed by requiring A αβ αβγ δ A γ δ ≥ 0 for any tensor A αβ . Based on the evolution equation (37), it can be concluded that D p,d αβ stands for the plastic strain-rate tensor. The quantity Z μν (which contains the effect of ψ discussed earlier) in the evolution equations for the kinetic and configurational entropy densities, (35) and (36), describes the split of the dissipation due to viscoplastic deformation among the two (partial) entropies. While the specification of this split is the main purpose of Z μν , it is worthy to note that nonequilibrium thermodynamics (i.e., the conditions (7)-(9)) require that Z μν also occurs in the driving force for viscoplastic deformation d λ , see (39).

Complete dynamic model, and its temperature and energy reformulations
Combining the results of the previous subsections, the complete set of evolution equations, (1) together with (2) and (3), one finds with the Cauchy stress tensor σ αγ given by (22), and the traceless plastic strain-rate tensor D p,d αμ , which is the deviatoric part of (38), with driving force (39). The fact that the plastic strain-rate tensor is traceless serves, in hindsight, as a justification of the assumption that the viscoplastic deformation is isochoric, and therefore one can indeed use ρ = ρ 0 / det F e for the mass density, as discussed in Sect. 3.1.
To make the proposed model suitable for practical applications, constitutive choices need to be made for the static aspects (namely e t ), for the transport/relaxation coefficients (namely λ s αγ , λ η αγ , μ x , and αβγ ε ), and for the quantities that regulate the exchange and coupling between the thermodynamic subsystems (namely Y αβ and Z αβ ), in terms of the system variables X . Particularly, Y αβ and Z αβ are of interest in this paper and will be discussed in the following section in more detail, when comparing our model with other models in the literature.
It is pointed out that the complete model (40)-(43) with (22) and (38)-(39) is symmetric with respect to the two subsystems, keeping in mind that the sign convention for the exchange tensors Y αβ and Z αβ is arbitrary. This fundamental symmetry of the model is also preserved when transforming to the other sets of variables discussed in the remainder of this section.
Similar to [16], to simplify the comparison with the literature, it is useful to transition from the entropy densities s and η as dynamic variables to other two equivalent variables, for example the temperatures, (15)- (16). There are two main consequences of this change of variables for the above model. First, the stress tensor (22) can be expressed in terms of the Helmholtz free energy per unit volume a = a(T, θ, F e ), which is the Legendre transform of the internal energy density e t = e t (s, η, F e ) [19]. If in addition a split of the internal energy e t into its kinetic (e) and configurational ( ) parts, analogous to (12), is introduced for later convenience, one finds for the Helmholtz free energy With this, it can be shown that the stress tensor (22) can be written in the form The second major consequence of the transition (s, η) → (T, θ) is the replacement of (41) and (42) by evolution equations for the two temperatures. Specifically, it can be shown that with D t = ∂ t + v γ ∇ γ the material (substantial) time-derivative, and with the irreversible contributions Finally, it is noted that the prefactors on the right-hand side of (50) and (51) are generalizations of inverse heat capacities. Instead of replacing the entropy densities by the two temperatures, one could also think of other replacements. In view of the literature, it is useful to consider instead the energy densities (per unit mass) of the kinetic and configurational subsystems,ê andˆ , respectively. From the above, it can be shown that the evolution equations for the partial energy densities become This completes the formulation of the model for physical aging and mechanical rejuvenation of glasses using subsystem entropies, and of alternative formulations of it in terms of other variables for the kinetic and configurational subsystems, namely temperatures or kinetic energies.

Discussion: comparison with the literature
In Sect. 3, we have presented a general model to describe physical aging and mechanical rejuvenation of amorphous solids, by considering the evolution of two (rather than one) thermal degrees of freedom, i.e., kinetic and configurational entropies, or the corresponding temperatures or internal energies, respectively. In particular, as shown in (41) and (42), these entropies are coupled nontrivially in the reversible dynamics through the tensor-valued function Y αβ , while the distribution of dissipated energy (due to viscoplastic deformation) between the two thermal subsystems is controlled by another tensor-valued function, Z αβ , in the irreversible dynamics. In this section, the model developed in Sect. 3 is compared with approaches in the literature, with the particular goal to learn about the common choices made for the quantities Y αβ and Z αβ . In this comparison that follows below, the models of the other references are always (transformed to and) written in an Eulerian formulation, using as much as possible the same symbols and definitions as in our approach, to simplify the comparison. If other symbols and definitions are used as compared to our approach, this is explicitly stated.

Comparison with Boyce et al. [3]
We start our discussion with considering the work of Boyce et al. [3], where an additional variable is introduced to account for the effects of aging and mechanical rejuvenation. However, [3] does not rely on any thermal (temperature, energy density, or entropy density) variable to describe a configurational subsystem. Rather, a general structural variable is employed, here denoted bys. For this quantity, used in the yielding kinetics, a differential equation of the following form is studied (eqn. (28) in [3]) with the functionā(s) describing the physical aging, whiler (s) describes how the plastic strain rateγ p gives rise to mechanical rejuvenation. Since both of these effects are also present in our modeling approach (as well as in the other models discussed below), one can say that they are in this sense equivalent. We now aim at relating the approach of [3] to ours, by identifying the quantities Y αβ and Z αβ in the following way. First, it is noted that in [3] the Cauchy stress tensor is based on the derivative of the Helmholtz free energy with respect to deformation [see their discussion around their Eq. (20) and (21)]. This implies in view of our relations Second, it is concluded from the explanations in [3] that the Helmholtz free energy does not depend on the structural variables. In the context of our model with effective temperature, this means that the Helmholtz free energy density (45) does not depend on the configurational temperature θ , and hence in the stress tensor expression (46)-(48) one can set σ η αβ = 0, since η = a ,θ | T,F e αβ = 0. Third, it is pointed out that in our model, we do not account for strain-hardening. Therefore, we should also not consider the effects of entropic strainhardening, which are actually included in [3]. While neglecting strain-hardening does in principle not imply that the stress tensor is completely free from (other) entropic contributions, the discussion around eqn. (21) in [3] suggests that one might indeed assume σ s αβ = 0, and so it follows Using this result in relation to the temperature evolution equation discussed in Sect. 2.5 in [3], and then comparing this with our temperature evolution Eq. (50) with (52), one obtains where we have used that the thermodynamic potential does not depend on the structural variables, which in our context implies that the thermodynamic potential does not depend on η, leading to T ,η = 0. Inserting expression (60) into our driving force for viscoplastic deformation (39) with (54), one obtains which corresponds to eqn. (7) used in the viscoplastic constitutive relation (11) in [3].

Comparison with Buckley et al. [10]
The work of Buckley et al. [10] formulates a model in terms of temperature evolution. In that model, the Cauchy stress tensor σ αβ is written as a sum of two contributions, specifically a (traceless) "bond-stretching" contribution σ b αβ and a "molecular-conformation" contribution σ c αβ . The thermal evolution equation (eqn. (30) in [10]) can then be written in the form withc the heat capacity and c the heat capacity step across the glass transition. Based on the two temperature equations in our model, (50)-(51), under the assumptions of no cross-coupling effects, T ,η = θ ,s = 0, and no heat fluxes, one can derive an evolution equation analogous to (62). By comparing (62) with our counterpart, particularly with respect to the D p αβ -terms (or conversely the (∇ β v α )-terms), one finds σ b αβ = σ e αβ + σ αβ . This result is valid for any Y αβ and Z αβ . Both of these quantities describe the exchange between the two thermal subsystems, and therefore they cannot be determined from a single evolution equation, as (62). However, both of these quantities can actually be determined by considering the constitutive relations for the Cauchy stress as well as for the driving force for viscoplastic flow. As far as the Cauchy stress is concerned, it is explained in [10] that the stress tensor is the sum of a mechanical visco-elastic Maxwell element and of a purely hyperelastic part. Therefore, the entire stress tensor can be derived from a Helmholtz free energy, and one thus concludes that Furthermore, let us consider the driving force for viscoplastic flow. By combining eqns. (15) and (26) in [10], one finds that the driving force for viscoplastic deformation is given by the (traceless) stress contribution σ b αβ = σ e αβ + σ αβ . Equating this driving force to T d αβ in our approach, with T a positive prefactor with units of temperature and d αβ given by (39) with (54), one can write Clearly, the second contribution (and for T = T also the first one) in this expression diverges as |θ − T | approaches zero. In view of the expression for the driving force for viscoplastic flow, αβ , this is not a problem, since Z αβ is multiplied there by a prefactor that vanishes as |θ − T | vanishes. However, we note that the split of the viscoplastic contribution between the kinetic and configurational subsystems, see (50) and (51) with (52) and (53), does suffer from this divergence in Z αβ as |θ − T | tends to zero. This issue deserves further attention in a future study.

Comparison with Kamrin and Bouchbinder [14]
In the work of Kamrin and Bouchbinder [14], a two-temperature-based continuum thermomechanics model was formulated to describe physical aging and mechanical rejuvenation of deforming amorphous solids. Specifically, by assuming weakly interacting kinetic and configurational subsystems, the evolution equations for corresponding energy densities (their eqns. (9) and (10)) read as follows [14] ρ with the split of the Cauchy stress tensor σ αβ = σ store αβ +σ dis αβ . The plastic stress-power is split into a contribution that represents the "stored plastic power" and affects the configurational subsystem, while the "plastic dissipation" gives a contribution to the kinetic subsystem [14]. Further, to provide a direct comparison between their and our approaches, we neglect in our approach the cross-coupling terms in the energy evolution equations (55) and (56), i.e.,ê ,η =ˆ ,ŝ = 0, leading toê ,ŝ = T andˆ ,η = θ . Doing so, the comparison of their (65)-(66) with our (55)-(56) is straightforward. Specifically, with the choice one can match the (∇ β v α )-terms in D tê (or alternatively in D tˆ ). On the other hand, by analyzing the D p αβterms in D tê (or alternatively in D tˆ ), in view of the Cauchy stress tensor together with (67), one can establish the connection between σ dis αβ and Z αβ , explicitly In [14], the split of σ αβ into σ store αβ and σ dis αβ stands for the split of the effect of viscoplastic dissipation into the kinetic and configurational subsystems. In our approach, this corresponds equivalently to the choice of Z αβ , as exemplified by (68). However, as we can see from the discussion above, from the evolution equations of the energies one can only determine Z αβ in relation to σ dis αβ , but not in absolute terms. To make further progress in this direction, one can now proceed by filling this expression for Z αβ into our driving force for viscoplastic deformation (39) with (54), leading to For a most general comparison of the driving force for viscoplastic flow used in [14] (eqn. (24) therein), σ dis,d αβ , with our approach, one needs to solve σ dis,d where T is a positive prefactor with units of temperature. Doing so, one obtains the specific relation Together with (68), this leads to In other words, [14] uses nontrivial choices for both Y αβ and Z αβ , in general. We conclude by noting that for the specific case discussed in Sect. 3 of [14] (see specifically assumptions in eqn. (17) therein), one finds σ e αβ = σ s αβ = σ η αβ = 0, which implies Y αβ = 0, and Z αβ becomes proportional to σ ,d αβ , by virtue of (67) and (71), respectively. [15] For the comparison of the work of Xiao and Nguyen [15] with our approach, we start by considering the evolution equations for the temperatures. In [15], it is assumed that the kinetic (e and s) and configurational ( and η) subsystem quantities do not depend on the temperature of other subsystem (θ and T , respectively). In our context, in view of (50)-(51), this implies T ,η = θ ,s = 0. Expressing the temperature equations (29) and (30) in [15], one obtains

Comparison with Xiao and Nguyen
These equations are now discussed from the perspective of our corresponding equations (50)-(51). One can naturally make the identifications for the subsystem heat capacities c g = T (T ,s ) −1 and c = θ (θ ,η ) −1 . Furthermore, by comparing the ∇ α v β -terms in the evolution equations for T (or alternatively θ ), one finds By comparing the D p αβ -terms in the evolution equations for T (or alternatively θ ), one finds It should be noted that the zero-choice for Y αβ can also be obtained by comparing the expression for the Cauchy stress in [15] [therein: based on eqn. (18) which corresponds to eqn. (25) in [15].

Comparison with Klompen et al. [6]
In the work of Klompen et al. [6], a model has been developed for describing the mechanical response of glassy polymers. Like in all approaches discussed above, also in [6] an additional variable, the so-called state parameterS, is introduced to describe the effects of physical aging and mechanical rejuvenation, by letting the yielding kinetics depend on this parameter. However, in contrast to all approaches described earlier in this section, the state parameterS in [6] is not specified by way of an evolution equation, but rather in an explicit form. Specifically, the state parameter is considered as a product of two explicit state functions, i.e., S =S aRγ , whereS a =S a (t, T ) represents the aging kinetics dependent on aging time t and temperature T , whereas the mechanical rejuvenation is described byR γ =R γ (γ p ), where γ p is an equivalent plastic strain. Specific procedures are developed and explained in [6,7] to determine the explicit functionsS a (·) andR γ (·). However, in order to establish links with the above discussed thermodynamic approaches to physical aging and mechanical rejuvenation, we briefly comment on how the approach in [6] can be transferred into an evolution equation for the state parameter. Simply differentiatingS with respect to time results iṅ One observes readily that the rate consists of two additive terms: one term representing physical aging and the other representing mechanical rejuvenation. The evolution equation (77) for the state parameterS strongly resembles that of the state parameters, i.e., (57). Particularly, the mechanical rejuvenation term is proportional to the rate of the (equivalent) plastic strain. However, while the rateγ p can instantaneously be expressed in terms of the plastic strain-rate tensor D p αβ , in the approach of Klompen et al. [6] there is an explicit reference to the plastic strain γ p itself. In practical numerical calculation, one can obtain that readily by time-integration of the rateγ p . However, from a fundamental modeling perspective it implies that γ p is elevated to the level of a dynamic variable. Since our thermodynamic approach does not treat γ p as an independent dynamic variable, we do not go into more detail about comparing our approach with that in [6].
Despite the fact that Klompen et al. [6] do not use an evolution equation for the state parameter, one can still strive to make a more close connection with the dynamic approach in this paper. As an illustrative example of how to achieve that task, we consider in the following physical aging in the absence of mechanical deformation, i.e.,R γ = 1, in which caseS =S a . In agreement with literature, also Klompen et al. [6] use a logarithmic dependence ofS a on waiting time t w , i.e.,S a ∝ ln(t w ). In contrast, in the approach in this paper, the physical aging in the absence of mechanical deformation is described by the μ x -related contribution in evolution of the configurational temperature (51) and (53), where we have assumed for simplicity θ ,s = 0. The goal is to make a certain choice for the kinetic function μ x such that, with the resulting solution θ = θ(t w ), one obtainsS a (θ (t w )) ∝ ln(t w ). If one assumes for illustration purposes thatS a (θ ) =S 0 +S 1 θ , one thus looks for a solution θ(t w ) ∝ ln(t w ). As it is shown in the following, this solution is induced by the choice which is always positive as required, μ x > 0, if a > 0, b > 0, and the configurational heat capacity θ/θ ,η > 0. Considering constant kinetic temperature T , the evolution equation for the temperature difference can be derived from (78) and (79), leading to which is identical to the evolution equation (1) discussed in [34]. Particularly, for constant a and b, the solution is given by [34] (59) and (61) with δ 0 = δ(t w = 0). Indeed, for e −b δ 0 1 and a b t w 1, one obtains logarithmic behavior, as desired. In contrast, for sufficiently long waiting times, a b t w 1, one obtains δ 0. In view of the system visiting ever deeper energy states in the course of physical aging, it seems more reasonable to have a limiting value for δ for long waiting times, rather than a logarithmic behavior on all time scales. 4.6 Discussion of the above comparison, and a further possibility for Y αβ and Z αβ In Sect. 4, we have compared our model with other models known in the literature [3,6,10,14,15], in order to identify some choices made in the literature for the coupling tensors Y αβ , Z αβ , and the driving force for viscoplastic deformation, d αβ . The results of that comparison are listed in Table 1. Notably, the cases of Buckley et al. [10] and Kamrin and Bouchbinder [14], discussed in Sects 4.2 and 4.3, respectively, contain a yet undetermined positive quantity T with the units of temperature. While the general expressions for Y αβ , Z αβ , and d αβ are given in the respective sections, in Table 1 only the special case T = T is listed, for illustration purposes. This choice is motivated by the relation implied by (59) and (61), namely T d αβ = σ d αβ . As shown in Table 1, only in the work of Kamrin and Bouchbinder [14] there is a nonzero entropy-coupling term Y αβ . Particularly, according to Table 1 and with the help of the stress tensor contributions (47) and (48), one can write Since according to [14] the kinetic and configurational subsystems are only weakly coupled, bothê andŝ can be assumed to be independent of θ , and thus the term in parenthesis in (84) can be written in the formê ,C e γ δ |ŝ. Using that Y αβ is of quasi-potential form as discussed in Sect. 3.2, one obtains from (84) Therefore, a possible choice is y 1 = T /2 and y 2 =ê, implying that both y 1 and y 2 depend only onŝ and C e . According to the temperature definition (15), one can thus writeê ,η |ŝ t ,C e = −ê ,ŝ | C e = −T for the case studied in [14], which implies ∂ y 2 ∂η ŝ t ,C e = −2y 1 .
Inspection of the Jacobi condition (27) immediately leads to the conclusion that the Jacobi identity is automatically satisfied, without any restriction on the thermodynamic potential for the kinetic subsystem. Therefore, the choices made in Kamrin and Bouchbinder [14] are fully compatible also with the Hamiltonian structure of the reversible dynamics, which is noteworthy as this criterion, namely the Jacobi identity, has not been used The tensor-valued functions, Y αβ and Z αβ , that couple the two thermal subsystems in the reversible and irreversible dynamics, respectively, are still quite general. Specifying their form, one can represent many other models known in the literature, as elaborated in detail in Sect. 4. An essential part of this paper is devoted to devise appropriate coupling of the two thermal subsystems in the reversible dynamics. In this respect, using GENERIC proves beneficial, since this framework respects the Hamiltonian structure, and therefore imposes severe constraints for the model formulation. Particularly, the importance of the Jacobi identity is highlighted, which represents the time-structure invariance of the reversible dynamics [27]. In this paper, the implications of the Jacobi identity on the reversible coupling (quantity Y αβ orŶ αβ , respectively) of the kinetic and configurational entropies have been examined. It seems that this is a research direction that to date has received only limited attention. To foster further research activities in this direction, some comments in relation to the developments in Sect. 3.2 are added. While the condition (25) for Y αβ is fully general, we believe that a large fraction of physically reasonable models are of quasi-potential typeŶ αβ = (1/y 1 )(∂ y 2 /∂C e αβ )|ŝ t ,η as discussed in Sect. 3.2, in which case the relevant condition on the two scalar-valued functions y 1 and y 2 is given by (27), or (28), respectively. It is noteworthy that a possible solution to comply with the Jacobi identity, (27) or (28), is given by ∂ y 2 ∂η ŝ t ,C e + 2y 1 = 0 .
As a matter of fact, it can be shown that (91) results in a vanishing material time-derivative of y 2 , completely analogous to the evolution ofŝ t , D tŝt = 0. Therefore, instead of working with the set of variables X = (m, s, η, F e ) with reversible entropy-coupling described by (91), one could perform a transformation of variables to X = (m,ŝ t , y 2 , F e ). The latter set is technically more convenient, since the Jacobi identity is known to be fulfilled for (m,ŝ t , F e ) [28]. Viewing this in a Lagrangian setting [35], it is immediately clear that including an additional variable (namely y 2 ) with dynamics given by (92) does not disturb the compatibility with the Jacobi identity. Since the Jacobi identity is invariant with respect to a transformation of variables X ↔ X , there is thus a clear understanding why the class of models described by (91) respects the Jacobi identity, namely by viewing y 2 as a fundamental dynamic variable. However, it is emphasized that from a physical, modeling perspective, it may be beneficial to keep the configurational entropy (η orη, respectively) as a fundamental variable instead. In any case, when formulating the thermomechanics of amorphous solids with two thermal subsystems, the condition (91) makes it possible to include a reversible entropy-coupling of quasi-potential form in a straightforward way. While the case of Kamrin and Bouchbinder [14] discussed above is an example thereof, the class of models covered by (91) is even richer.