How interface size, density, and viscosity affect creep and relaxation functions of matrix-interface composites: a micromechanical study

Matrix-inclusion composites are known to exhibit interaction among the inclusions. When it comes to the special case of inclusions in form of flat interfaces, interaction among interfaces would be clearly expected, but the two-dimensional nature of interfaces implies quite surprising interaction properties. This is the motivation to analyze how interaction among two different classes of microscopic interfaces manifests itself in macroscopic creep and relaxation functions of matrix-interface composites. To this end, we analyze composites made of a linear elastic solid matrix hosting parallel interfaces, and we consider that creep and relaxation of such composites result from micro-sliding within adsorbed fluid layers filling the interfaces. The latter idea was recently elaborated in the framework of continuum micromechanics, exploiting eigenstress homogenization schemes, see Shahidi et al. (Eur J Mech A Solids 45:41–58, 2014). After a rather simple mathematical exercise, it becomes obvious that creep functions do not reflect any interface interaction. Mathematical derivation of relaxation functions, however, turns out to be much more challenging because of pronounced interface interaction. Based on a careful selection of solution methods, including Laplace transforms and the method of non-dimensionalization, we analytically derive a closed-form expression of the relaxation functions, which provides the sought insight into interface interaction. The seeming paradox that no interface interaction can be identified from creep functions, while interface interaction manifests itself very clearly in the relaxation functions of matrix-interface materials, is finally resolved based on stress and strain average rules for interfaced composites. They clarify that uniform stress boundary conditions lead to a direct external control of average stress and strain states in the solid matrix, and this prevents interaction among interfaces. Under uniform strain boundary conditions, in turn, interfacial dislocations do influence the average stress and strain states in the solid matrix, and this results in pronounced interface interaction.

List of symbols a Coefficient of ordinary differential equation in γ a 1 Radius of first interface family a 2 Radius of second interface family A I Relaxation capacity associated with characteristic time τ relax,I A II Relaxation capacity associated with characteristic time τ relax,II A j Third-order strain-to-dislocation downscaling tensor, quantifying the concentration of macroscopic strain into the average dislocation encountered across the interfaces; j = 1, 2 refers to the two interface families A Σ j Third-order influence tensor describing the influence of macroscopic stress on the dislocation of interfaces; j = 1, 2 refers to the two interface families b Coefficient of ordinary differential equation in γ B j Third-order influence tensor describing the influence of interface traction on the macroscopic stress; j = 1, 2 refers to the two interface families B Σ j Third-order influence tensor describing the influence of interface traction on the macroscopic strain; j = 1, 2 refers to the two interface families c j Constant used in Laplace transformation method, related to dislocation of interface family j = 1, 2 C j Integration constant used in the elimination-non-dimensionalization method, related to the dislocation of interface family j = 1, 2 C 1 Constant appearing in relaxation function C 2 Constant appearing in relaxation function C j Domain of interface belonging to interface family j = 1, 2 C s Fourth-order elastic stiffness tensor of solid Inverse of C s C hom Fourth-order homogenized stiffness tensor of matrix-interface composite Inverse of C hom d j Interface density parameter of interface family j = 1, 2 D j Integration constant used in the elimination-non-dimensionalization method, related to dislocation of interface family j = 1, 2 D 1 Constant appearing in relaxation function D 2 Constant appearing in relaxation function e j Constant used in Laplace transformation method, related to the dislocation of interface family j = 1, 2 e 1 , e 3  Number of interfaces per unit volume of a matrix-interface composite; j = 1, 2 refers to the two interface families n Outward unit normal at any point of the boundary of the volume element representing the composite material n j Unit vector oriented orthogonal to the interface family j = 1, 2 R xzxz Relaxation function s As subscript: index for solid phase s As variable: Laplace space parameter T j Interface traction vector acting on interface phase j = 1, 2 T i Fourth-order morphology tensor for 2D interface inclusion ("sharp crack" morphology) Second-order tensor contractioṅ • Partial derivative with respect to time ("rate"), of quantity "•" ⊗ Dyadic product

Introduction
Matrix-inclusion composites are known to exhibit interaction among the inclusions [1][2][3][4][5][6][7]. When it comes to the special case of inclusions in form of flat interfaces, i.e., to that of matrix-interface composites, interaction among interfaces would be clearly expected as well; however, the two-dimensional nature of interfaces is responsible for particularly surprising interaction properties [8,9], reminiscent of the situation encountered with micro-cracked materials [10,11]. The present contribution tackles the question of how interaction among microscopic viscous interfaces, believed to be at the origin of creep of hydrated biomaterials and geomaterials [12][13][14], affects the overall macroscopic creep and relaxation functions of matrix-interface composites. To this end, we analyze matrix-interface composites consisting of a linear elastic solid matrix and parallel viscous interfaces, and we consider that the creep and relaxation behavior of such composites results from micro-sliding within adsorbed fluid layers filling the interfaces. The latter idea was recently elaborated in the framework of continuum micromechanics [15], based on eigenstress homogenization schemes [16].
Extending the analysis of [15]-where we mainly focused on identical interfaces-we here consider interaction among two classes of interfaces (referred to as "interface families"), differing in interface size, viscosity, and density. The theoretical basis for our analysis is the topic of Sect. 2, where we briefly recall fundamental state equations governing the time-dependent behavior of matrix-interface composites comprising two-dimensional interfaces filled with viscous fluids [15]. This review comprises both types of so-called Hashin boundary conditions [17], i.e., uniform strain boundary conditions and uniform stress boundary conditions, as means for the study of relaxation and creep properties, respectively. The latter properties are comparatively simple to be derived from the microscopic interface behavior, and the corresponding result can be directly recalled from our previous study [15]. The derivation of relaxation functions turns out as formidable mathematical task, which Sect. 3 of the present paper is fully devoted to: In order to (i) derive a compact closed-form solution of the sought relaxation functions and to (ii) provide detailed insight into interface interaction, we carefully select solution methods for coupled systems of linear differential equations. They are: Laplace transformation, a decoupling strategy in time domain based on an elimination scheme, and the method of non-dimensionalization. The analytically derived relaxation functions are evaluated numerically, in order to (i) study their sensitivity with respect to interface size, density, and viscosity of two interface families, and to (ii) provide insight into interface interaction. Comparing creep and relaxation functions with the aim to identify the reason for interface interaction is the focus of Sect. 4. There, we discuss the seemingly paradox situation that no interface interaction can be identified from the mathematical structure of the creep functions, while interface interaction is clearly manifested in the relaxation functions. In fact, this situation can be understood when recalling the stress and strain average rules for materials hosting interfaces, relating to loading in terms of uniform stress boundary conditions and uniform strain boundary conditions, respectively. Section 4 closes the paper with final conclusions.

Matrix-interface micromechanics for different interface families and review of creep functions
We consider a matrix-interface composite consisting of a linear elastic solid matrix and of two families of embedded 2D, circular, viscous interfaces which are parallel to the x, y-plane (Fig. 1). The material behavior of the solid matrix follows the generalized Hooke's law, In Eq. (1), σ (x) and ε(x) denote Cauchy micro-stresses and linearized micro-strains at any position x in the volume Ω s of the solid phase, and C s denotes an isotropic elastic stiffness tensor, In Eq. (2), E s and ν s , respectively, denote Young's modulus and Poisson's ratio of the isotropic matrix. In addition, I vol and I dev stand for the volumetric and the deviatoric part of the symmetric fourth-order identity tensor I (see the "Appendix" for components of these tensors).
The two interface families are denoted by indexes 1 and 2, and they comprise interfaces with radius a 1 and a 2 , respectively ( Fig. 1). Denoting the number of interfaces per unit volume of the composite as N 1 and as N 2 , respectively, the interface density parameters d 1 and d 2 read, according to Budiansky and O'Connell [18], as d 1 = N 1 a 3 1 and d 2 = N 2 a 3 2 . The interfaces are considered to host fluids adsorbed to electrically charged solid surfaces [15]. Accordingly, the polar fluid molecules (such as water) build up layers, i.e., such fluids are in a "glassy" or "liquid crystal" state, right in between the long-range positional and orientational order found in solids and the long-range disorder found in liquids. The lubricant effect of such fluid layers promotes gliding of the solid surfaces along the fluid sheets [12,13,19,20]. This fluid behavior is mathematically described by relations between components of (i) interfacial dislocation vectors [[ξ ]] 1 and [[ξ ]] 2 quantifying the displacement jumps across the interfaces, and of (ii) interfacial traction vectors acting on the interface planes, T 1 and T 2 . Molecular ordering-related joining forces prevent the interfaces from opening, i.e., displacement jumps in the interface normal direction e z (see Fig. 1) vanish in both interface families, As for the in-plane behavior, the adsorbed fluids are considered to result in linear viscous relations between the components of the traction vectors (T 1,x and T 1,y as well as T 2,x and T 2,y ), and the corresponding dislocation rates, with interface viscosities η i,1 and η i,2 as proportionality factors, In Eq. (4), the dot stands for the time derivative, At this point, it appears suitable to shortly discuss the physical meaning of the interface viscosities as introduced in Eq. (4): as seen from this equation, the interface viscosities relate the (average) tractions acting on the interface plane, to the rate of the (average) displacement jumps encountered when passing through the interface, orthogonal to the interface plane. Accordingly, the dimension of the interface viscosities reads as ML −2 T , with M, L, and T being abstract positive numbers which are related to different choices of units of measurements for mass, length, and time [21]. The aforementioned dimension differs by a factor of L −1 from the standardly defined bulk viscosity of fluids, with dimension ML −1 T . This results from the fact that the observation scale relevant for Eq. (4) is so large that the interface thickness does not appear any more as finite measure, but instead, as infinitely small. Hence, the thickness of the interface does not explicitly appear any more, and the interface viscosity follows to be of the form "bulk viscosity over interface thickness." Accordingly, when zooming into an even smaller scale, where the interfaces would again appear as a three-dimensional domain, then the bulk viscosity of the latter would result as the product of the aforementioned interface viscosity and the thickness of the interface. Typical viscosity magnitudes for the case of (water-filled) interfaces between hydroxyapatite crystals in bone have been recently given in [22]: an interface viscosity of η i = 1.83 × 10 12 GPa s m −1 ; with a corresponding bulk viscosity of the fluid filling the interface, of η b = 1.83 × 10 12 Pa s; and an interface thickness of 1 nm.
In the sequel, we recall fundamental relations of interfacial micromechanics [15], which were derived in the framework of eigenstress homogenization [16]. Thereby, we consider that the matrix-interface composite of Fig. 1 is subjected to a specific type of so-called Hashin boundary conditions [17], i.e., either to uniform strain boundary conditions or to uniform stress boundary conditions, respectively.

State equations for uniform strain boundary conditions
Consider that the boundary ∂Ω of the studied matrix-interface composite ( Fig. 1) is subjected to linear displacement functions which are related to a uniform macrostrain state E by In addition, consider that traction forces T 1 and T 2 prevail in the interface phase. The corresponding macroscopic state equation expresses the macrostress Σ as a function of the macrostrain E and of the interface traction vectors T 1 and T 2 , reading as [15] In Eq. (7), C hom denotes the homogenized stiffness tensor, while B 1 and B 2 , respectively, denote influence tensors, (see the "Appendix" for components of these tensors). The formulae of the "Appendix" are based on so-called Mori-Tanaka homogenization schemes [23], where the mechanical interactions between inhomogeneities (here viscous interfaces) in a homogeneous matrix are considered not up to complete detail, but rather in an average sense, see, e.g., [1,2,24] for more detailed discussions: each interface (or more generally, each inhomogeneity) is formally embedded in a fictitious matrix undergoing the mean strains of the actual material matrix, these strains being different from the macroscopic strains acting on the entire RVE; the matrix and inhomogeneity strains are enforced to fulfill the strain average rule for the entire RVE. The Mori-Tanaka scheme has been repeatedly used for modeling the behavior of cracked materials [25][26][27]. In particular, it has been shown that in the case of sharp cracks (exhibiting the same morphology as the interfaces), the aforementioned average consideration of inhomogeneity interaction leads to remarkably precise results for the overall homogenized stiffness properties, when compared to solutions accounting for fully explicit interactions, as has been shown by Kachanov et al. [25]. Hence, we expect the results of the present contribution to be only marginally affected when replacing the Mori-Tanaka estimates of the "Appendix" by more explicit (and hence also much more complex) formulations. This expectation is supported by the theoretical derivation of Ponte Castañeda and Willis [28] that the Mori-Tanaka scheme actually considers spatial distributions of inhomogeneities following the same ellipsoidal shape as the inhomogeneities themselves. For material systems targeted at by our approach, including calcium silicate hydrates in concrete [29], extrafibrillar spaces in bone [30], or montmorillonite interlayers in clay [31], these spatial distributions appear as natural and appropriate. The remaining two state equations establish relations between dislocation vectors of the two interface families, [[ξ ]] 1 and [[ξ ]] 2 , on the one hand, and the macrostrain E as well as the interface traction vectors T 1 and T 2 , on the other hand. These concentration-influence relations read as [15] [ In Eqs. (8) and (9), A 1 and A 2 denote third-order strain-to-dislocation downscaling tensors of dimension L (length)-and in this regard, they are different from the dimensionless fourth-order strain-to-strain downscaling tensors called strain concentration tensors in classical continuum micromechanics [2], while D 11 , D 12 , D 21 , and D 22 represent influence tensors (see the "Appendix" for components of these tensors). Notably, D 12 and D 21 account for the interaction of the two interface families, i.e., the components of these two tensors quantify the influence of interface tractions of one interface family, on the dislocation evolution of the other interface family, and vice versa. This interaction will be the central topic of Sect. 3, where we derive relaxation functions of the matrix-interface composite of Fig. 1.

State equations for uniform stress boundary conditions
Consider that the boundary ∂Ω of the studied matrix-interface composite ( Fig. 1) is subjected to a field of traction vectors T which are related via Cauchy's formula to a uniform macrostress state Σ, where n(x) denotes the outward unit normal at any point x of the boundary of the composite. In addition, consider that traction forces T 1 and T 2 prevail in the interface phase. The corresponding macroscopic state equation expresses the macrostrain E as a function of the macrostress Σ and of the interface traction vectors T 1 and T 2 . It follows from solving Eq. (7) for the macrostrain E [15], and it reads as In Eq. (11) B Σ 1 and B Σ 2 , respectively, denote influence tensors (see the "Appendix" for components of these tensors).
The remaining two state equations establish relations between dislocation vectors of the two interface families, [[ξ ]] 1 and [[ξ ]] 2 , on the one hand, and the macrostress Σ as well as the interface traction vectors T 1 and T 2 , on the other hand. These concentration-influence relations follow from specifying Eqs. (8) and (9) for (11), and they read as In Eqs. (12) and (13), A Σ 1 , A Σ 2 D Σ 11 , and D Σ 22 represent influence tensors (see the "Appendix" for components of these tensors). Very remarkably, Eqs. (12) and (13) do not contain influence tensors D Σ 12 and D Σ 21 , i.e. interface tractions of one interface family do not influence the evolution of the dislocations of the other interface family, and vice versa. This renders the situation with uniform stress boundary conditions (10) as fundamentally different from the situation encountered with uniform strain boundary conditions (6), compare the concentration influence relations (12) and (13), with those of (8) and (9). Clarification of the reason for this interesting difference is one of the central aims of this paper. To this end, we derive and compare macroscopic creep and relaxation functions defined on the homogenized matrix-interface composite. While relaxation functions are reserved for Sect. 3, we already here quickly recall the derivation of creep functions as anticipated in [15].

Review of creep functions
In order to derive the creep functions related to the matrix-interface composite of Fig. 1, we consider a stress state which is suddenly imposed at time t = 0 and kept constant thereafter (t > 0), see the uniform stress boundary conditions (10). Since time-dependent behavior of the composite is related to in-plane dislocations of the interfaces, see the viscous interface law (4), and given the structure of state equation (7) together with the non-zero components according to (123-128), the time-dependent macroscopic stresses obey the form of pure shear, Specification of the concentration-influence relations (12) and (13) for (14), as well as consideration of the components of the involved influence tensors according to the "Appendix," and of the viscous interface behavior according to (4), deliver the following two differential equations for the time evolutions of the in-plane dislocation functions Since Eqs. (15) and (16) are uncoupled, the solution for the two sought dislocation histories is straightforward, and it reads as The creep function is obtained in two steps. First, the interface traction histories T 1,x and T 2,x are calculated by specifying the viscous interface law (4) for the time derivatives of the dislocation histories (15) and (16). Thereafter, the macroscopic state equation (11) is specified for the imposed stress state (14) and for the interface traction histories T 1,x and T 2,x obtained in the first step. Considering the components of the homogenized stiffness tensor and the involved influence tensors according to the "Appendix," and comparing the resulting expressions with deliver the sought creep functions as [15] The creep functions (20) contain two exponentials underlining the existence of two different characteristic creep times. They read as Since the concentration-influence relations (12) and (13) are free of interaction terms, it was mathematically quite simple to derive the creep functions of the matrix-inclusion composite of Fig. 1. This will be different when it comes to the determination of relaxation functions, as discussed next.

Derivation of relaxation functions for interacting interfaces of different size, viscosity, and density
In order to derive relaxation functions of the matrix-interface composite of Fig. 1, we consider a deformation state which is suddenly imposed at time t = 0 and kept constant thereafter (t > 0), see the uniform strain boundary conditions (6). Since time-dependent behavior of the composite is related to in-plane dislocations of the interfaces, see the viscous interface law (4), and given the structure of the state equation (11) together with the non-zero components according to (139-142), the time-dependent macroscopic strains obey the form of pure shear [note the analogy to Eq. (14) considered for the derivation of the creep functions], Specification of the concentration-influence relations (8) and (9) for (22), and consideration of the viscous interface behavior according to (4), deliver the following set of coupled, linear, inhomogeneous, first-order, ordinary differential equations (with constant coefficients) for the dislocation histories of the two interface families:  (23) and (24). This underlines the coupled nature of Eqs. (23) and (24). In order to simplify this situation, the differential equation (23)  from the differential equation (24). Vice versa, the differential equation (24) is solved for [[ξ ]] 2,x , and the resulting expression is used for elimination of [[ξ ]] 2,x from the differential equation (23). This yields This set of differential equations can be brought into the more comfortable forṁ based on the following definitions of constants: The most popular method for solving a system of differential equations like the one given by (27) and (28) is based on Laplace transformation. This is described next.

Determination of dislocation histories based on Laplace transformation
Let us recall that Laplace transformation L of functions γ i (t) yields functions Γ i (s), according to the definition where s is a complex number. Equation (32) implies the following transformation rules for first-order and second-order time derivatives of functions γ i (t) [32]: Accordingly, determination of dislocation histories by means of Laplace transformation becomes more straightforward when considering the time derivatives of Eqs. (27) and (28), resulting in the following set of coupled, linear, homogeneous, second-order, ordinary differential equations (with constant coefficients): Relations (33) underline (i) that initial conditions are considered already during Laplace transformation and (ii) that these initial conditions need to be formulated in the dimensionless dislocations and their rates (according to (30.3) and (31.3)). At time t = 0, i.e., at the time instant of sudden loading, the dislocations vanish in all interfaces The (dimensionless) dislocation rates right after sudden loading at time t = 0 follow simply from specification of (27) and (28) for (36), and from solving the resulting expressions forγ 1 (0) andγ 1 (0), yieldinġ The Laplace transformations of Eqs. (34) and (35) are carried out under consideration of transformation rules (33), as well as of initial conditions (36) and (37). This yields Rearranging terms in (38) and (39) delivers the following system of two algebraic equations: ⎡ The solution of (40) yields the functions Γ 1 (s) and Γ 2 (s) as Notably, the two solutions (41) and (42) are structurally identical. In other words, the solution for Γ 1 (s) according to (41) can be converted into the solution for Γ 2 (s) according to (42), and vice versa, simply by permutation of indexes, i.e., by setting 1 → 2 and, at the same time, 2 → 1. Consequently, it is sufficient to back transform only Γ 1 (s) in order to obtain the solution for γ 1 (t) in physical time domain, because γ 2 (t) follows simply from subsequent index permutation. For the back transformation of Γ 1 (s) from Laplace space to physical time domain, we will use the following basic rules [32]: We are left with modifying Γ 1 (s) according to (41), such that the back transformation rules (43)-(45) can be applied. This is done in two steps. At first, we transform (41) into the format by the definition of new coefficients a, b, c 1 , and e 1 , reading as In the second step, we expand the expression for Γ 1 (s) according to (46), such that we arrive at the following expression: Back transformation of (48) into physical time domain yields, based on transformation rules (43)-(45), In order to replace the trigonometric functions in (49) by exponential functions, we use the following two Euler's formulas [32]: Specifying (49) for (50) and (51), together with consideration of Finally, we consider that the product rule for exponential functions exp(x) exp(y) = exp(x + y), yielding γ 1 (t) as γ 2 (t) is derived from (54) by means of index permutation, i.e., by replacing index 1 by index 2 and at the same time, by replacing index 2 by index 1, which leads to The coefficients c 2 and e 2 follow from c 1 and d 1 according to (47), by means of the same index permutation, so that we arrive at The two solutions (54) and (55) clearly indicate interaction among the two different interface families, but the mechanical reason remains somewhat obscure. In more detail, we observe that: -The square-root expressions in Eqs. (54) and (55) originate from preparation of Γ 1 (s) according to (51), for back transformation into physical time space, see (49). This square-root expression combines information on both interface families and, therefore, represents a coupling term. The mechanical meaning of it, however, is not directly obvious, because of the markedly mathematical character of the derivation. -The obtained solutions for the dislocation histories are still quite mathematically expanded, and more compact expressions appear to be out of direct reach.
Both aspects provide the motivation to more deeply investigate the solution of Eqs. (27) and (28), by (i) solving the governing system of coupled differential equations (27) and (28) in physical time space by an uncoupling strategy based on an elimination scheme, and by (ii) combining this approach with the method of non-dimensionalization, in order to arrive at a compact closed-form solution of the relaxation function. This is described next.

Revisiting dislocation histories based on an elimination scheme combined with the method of non-dimensionalization
Our first aim is to combine the two governing differential equations (27) and (28), such that an uncoupled differential equation, exclusively in γ 1 (t), is obtained. To this end, we solve (27) for γ 2 (t), and we calculate its first-order time derivative, Insertion of (57) and (58) into (28) eliminates γ 2 (t) as well asγ 2 (t), and delivers a linear, inhomogeneous, second-order, ordinary differential equation exclusively in γ 1 (t), which reads, after rearranging and collecting terms, as It is interesting to interpret Eq. (59) from a mechanical viewpoint. To this end, we recall that the original problem is coupled, see (27) and (28). The used elimination scheme delivered an uncoupled differential equation for γ 1 (t), see (59), but this mathematical modification, of course, does not change the underlying physics. The coupled nature of the problem manifests itself in the coefficients of the uncoupled differential equation, see Eq. (59). Notably, the same two coefficients were found during the Laplace solution, compare the abbreviation a and b in (59) with the definitions (47). The solution of the differential equation (59) contains two parts: the particulate integral γ 1 p (t) and the solution of the homogeneous differential equation, called complementary function γ 1h (t), We start with the particulate integral γ 1 p (t). Since the right-hand side of (59) is time independent, also the sought particulate integral is a constant. It follows from specifying (59) for γ 1 (t) = γ 1 p , forγ 1 p (t) = 0, and forγ 1 p (t) = 0, as well as from solving the resulting expression for γ 1 p , as Notably,γ 1 p (t) = 0 andγ 1 p (t) = 0 underline that γ 1 p is a stationary (time independent) solution, such that γ 1 p can be interpreted as the dislocation, which is asymptotically reached after infinite time, i.e., the dislocation that is finally reached once the relaxation process has come to an end. In addition, we note that γ 1 p is identical to constant c 1 , which we have introduced during the derivation of the Laplace solution, compare (61) with c 1 in (47). We are left with the determination of the complementary function γ 1h (t). To this end, we rewrite the differential equation (59) under consideration of a and b according to (47), and we set the right-hand side equal to zero, in order to obtain a homogeneous differential equation, As for the solution of (62), we make the following ansatz involving an exponential function and two constants β and λ: Specifying the differential equation (62) for ansatz (63) yields Dividing Eq. (64) by β exp(λ t) results in the following quadratic equation for λ: The two roots of the quadratic equation (65) are two independent solutions for λ, i. e.
Equation (66) implies that the sought complementary function γ 1h (t) consists of two terms, containing two exponentials which are multiplied by two integration constants C 1 and D 1 , The complete solution for the dislocation history γ 1 (t) follows from specification of (60) for the particulate integral from Eq. (61) and the complementary function from Eq. (67) as While noting that integration constants C 1 and D 1 are conceptually to be identified from the initial conditions (36) and (37), we emphasize that the solution (68) exhibits exactly the same structure as the solution derived with the Laplace transformation method, compare (68) with (54). Still, the used uncoupling strategy based on an elimination scheme provides valuable insight into the question why a square root shows up in the exponential function. It stems from a quadratic function which needs to be solved during the derivation of the complementary function. Still, the square-root terms render the derivation of a compact version of a closed-form solution for the relaxation function difficult. In order to further improve the situation, we revisit the derivation of the complementary function γ 1h (t) based on the method of non-dimensionalization. This approach allows for reducing the mathematical complexity of the solution to a possible minimum. Given that time t is the parameter of the studied relaxation problem, the idea of non-dimensionalization is to set the dimensional time t equal to dimensionless time τ , multiplied with an arbitrary time constant t c , t = τ t c .
Appropriate choice of the time constant t c will allow us to reduce the mathematical complexity of the solution to the sought minimum. Consideration of (69) in the form dt = dτ t c underlines that every derivative with respect to dimensional time t can be replaced by a derivative with respect to dimensionless time τ , multiplied with 1/t c . Applying this strategy to the homogeneous differential equation (62) yields Since Eq. (71) is completely analogous to Eq. (62), and because (62) has led to the two solutions for λ according to (66), we conclude that avoiding a square-root expression in the solution function γ 1h (t) can be achieved by choosing the time constant t c according to the following condition: Specifying condition (72) for a and b according to (70), respectively, and solving the resulting expression for the time constant t c delivers t c = 2 With help of condition (72), and the corresponding special choice (73) for the time constant t c , the solution of the differential equation (71) can now be retrieved in a much simpler form, which reads as whereby we followed the same line of reasoning that has led us from Eq. (59) all the way through Eq. (67). Notably, γ 1h (τ ) in Eq. (74) is a function of dimensionless time τ . In order to obtain an expression for γ 1h (t), i. e. a function where dimensional time t appears as the argument, we specify (74) for τ = t/t c , The complete solution for the dislocation history γ 1 (t) follows from specification of (60) for the particulate integral from Eq. (61) and the complementary function from Eq. (75), yielding The two integration constants C 1 and D 1 are identified from the initial conditions (36) and (37). To this end, we calculate the time derivative of γ 1 (t) given in Eq. (76) aṡ Specifying the initial conditions (36) and (37) for γ 1 (t) given in (76) and forγ 1 (t) given in (77) yields the following two equations for C 1 and D 1 as Solving the initial conditions (78) and (79) for the integration constants C 1 and D 1 delivers The solution for γ 2 (t) is found by analogy to (76), i.e., by permutation of indexes: The integration constants C 2 and D 2 follow from initial conditions (36) and (37) and by analogy to C 1 and D 1 given in (80) and (81), respectively, as The dislocation histories (76) and (82) exhibit a sufficiently compact mathematical form, ready for the derivation of a closed-form representation of the relaxation function. This is described next.

Relaxation functions
Relaxation functions are obtained in two steps. At first, interface traction histories T 1,x and T 2,x are calculated. To this end, the relations between γ 1 (t) and [[ξ 1 ]], as well as between γ 2 (t) and [[ξ 2 ]], as defined in (31.1-3), are solved for the dislocation histories: The sought interface traction components T 1,x and T 2,x follow from specifying the viscous interface law (4) for the time derivatives of the dislocation histories (85): In the second step, the macroscopic state equation (11) is specified for the imposed strain state (22), which yields, under consideration of the components of the homogenized stiffness tensor and the involved influence tensors according to the "Appendix," the following scalar equation for the time evolution of the macroscopic shear stress Σ xz : Specifying (87) for the tractions (86) yields, under consideration of the definitions (31.1-3), the following compact result: Comparing Eq. (88) with Σ xz = R xzxz (t) 2E xz (89) delivers the sought relaxation functions as The relaxation functions (90) contain two exponentials which are underlining the existence of the following two different characteristic relaxation times: Equations (95) and (96) together with the characteristic evolution of a relaxation test (discussed in the sequel) provide the motivation to re-formulate the relaxation functions (90) finally as with Noting that A I + A II = 1 allows for the following mechanical interpretation of Eq. (97). Specifying (97) for Equation (101) states that the effective stiffness of the composite is equal to the solid stiffness at the beginning of a relaxation test. Specifying (97) for t = ∞ yields Comparison of Eqs. (101) and (102) shows that Δμ denotes the loss of effective stiffness of the composite during the relaxation test. Finally, we note that A I and A II stand for relaxation capacities associated with the two characteristic relaxation times.
For the transition to micromechanical quantities, i. e. to Young's modulus E s and Poisson's ratio ν s of the solid matrix, to interface sizes a 1 and a 2 , to interface viscosities η i,1 and η i,2 , as well as to interface densities d 1 , and d 2 , consider the definitions of a according to (47), of t c according to (73) and of μ s , μ 1 , μ 2 , η 1 , η 2 according to (31). Exemplarily, this transition is carried out for the two characteristic relaxation times (95) and (96). This delivers

Study of interface interaction in a relaxation test
Because the concentration-influence relations (8) and (9) contain interaction terms, it was mathematically quite challenging to derive the relaxation functions of the matrix-inclusion composite of Fig. 1. In the following, we evaluate our analytical results , so as to gain insight into the nature of interface interactions. To this end, we consider the example of creep of bone stemming from viscous interfaces between hydroxyapatite crystals, slightly extending the discussion given in [14]. In this reference, needle-type bulk viscoelastic mineral phases were introduced as part of the porous polycrystal forming the extrafibrillar space of bone [30]-for further details on the mechanics of such porous polycrystals, we refer to [33]. Herein, we regard the viscoelasticity of these aforementioned mineral phases stemming from viscous interfaces being embedded, at a yet smaller scale, within the bulk crystal phases. Hence, the solid elastic matrix between the viscous interfaces would exhibit the elastic properties of hydroxyapatite, namely [34] E s = 114 GPa, ν s = 0.27.
In order to study the sensitivity of the relaxation function (97) with respect to changes in size, viscosity, and density of two interface families, we consider that the sum of the two interface densities amounts to 0.3, We will investigate all cases between d 1 = 0.     which is consistent with the interface viscosity of η i = 1.83 × 10 12 GPa s m −1 as given in [22], and an interface radius of 50 nm, consistent with electron microscopic images [30]. In order to study interaction among two interface families exhibiting different interface sizes and viscosities, the corresponding product a 2 η i,2 , referring to the second interface family, is considered to be a multiple of the amount given in (108), where n will be set equal to 2, to 5, and to 10, respectively, see Figs. 2, 3, and 4. (109) implies the following interpretation of the multiplication factor n: n is equal to the ratio of the two interface radii, n = a 2 /a 1 , provided that the two interface viscosities are the same, η i,1 = η i,2 . n is equal to the ratio of the two interface viscosities, n = η i 2 /η i 1 , provided that the two interface radii are the same, a 1 = a 2 .
As for discussing the characteristic relaxation times (95) and (96), see also (103), (104), and (105), and the corresponding relaxation capacities (99) and (100), it is useful to study limit cases d 1 → 0 and d 2 → 0. Considering the limit case that all interfaces belong to the first family (d 2 → 0) yields lim d 2 →0 Equations (110) and (111) indicate that the relaxation function (97) degenerates, for d 2 → 0, to an expression involving only one exponential with τ relax,I as argument, and this characteristic time is equal to the one which was derived in [15] for interfaces of identical size and viscosity, see the white square in Figs. 2a, b, 3a, b, as well as 4a, b. The second characteristic time, τ relax,II , has, in the limit d 2 → 0, no physical meaning, because the related relaxation capacity A II vanishes, see (111) and the black circles in Figs. 2a, b, 3a, b, as well as 4a, b. Considering the other limit case that all interfaces belong to the second family (d 1 → 0) yields by analogy to (110) and (111) lim Again, Eqs. (110) and (111) indicate that the relaxation function (97) degenerates, for d 1 → 0, to an expression involving only one exponential with τ relax,II as argument, and this characteristic time is equal to the one which was derived in [15] for interfaces of identical size and viscosity, see the white circle in Figs. 2a, b, 3a, b, as well as 4a, b. The first characteristic time, τ relax,I , has, in the limit d 1 → 0, no physical meaning, because the corresponding relaxation capacity A I vanishes, see (112) and the black squares in Figs. 2a, b, 3a, b, as well as 4a, b. Between the two discussed limits, a continuous transition of characteristic relaxation times and corresponding relaxation capacities is obtained, see Figs. 2a, b, 3a, b, as well as 4a, b. Also the relaxation functions exhibit a continuous transition from the limit case d 1 = 0.0 and d 2 = 0.3, i. e. from one interface family with a larger interface radius and/or larger interface viscosity, to the other limit case d 1 = 0.3 and d 2 = 0.0, i. e. to one interface family with a smaller interface radius and/or smaller interface viscosity, see Figs. 2c, 3c, and 4c. Since dislocation processes evolve faster at smaller scales and slower at larger scales, the asymptotic ("relaxed") state is reached the faster, the smaller or less viscous interfaces are present at the microstructure of the composite.
In the mentioned transition regime from d 1 = 0.0 and d 2 = 0.3 to d 1 = 0.3 and d 2 = 0.0, our analytical results strongly underline significant interaction between the two interface families. The reason for this behavior is discussed next.

Identification of the mechanism responsible for interaction among interfaces and concluding remarks
In order to study interaction between interfaces of different size, viscosity, and density, we have recalled the derivation of creep functions, see (20), and we have provided new derivations of the relaxation functions, see (97), both for a matrix-interface composite containing two interface families. The creep functions (20) contain two exponentials underlining the existence of two different characteristic creep times, see Eqs. (21). The latter equations highlight that each of the two interface families is associated with one of the two characteristic creep times. In other words, the characteristic creep times suggest that the two interface families do not interact. The relaxation functions (90) also contain two exponentials which are-again-underlining the existence of two different characteristic relaxation times, see Eqs. (103). However, contrary to the case discussed before, these equations now highlight that properties of both interface families influence both characteristic relaxation times. In other words, the characteristic relaxation times suggest that the two interface families are interacting, see also Figs. 2, 3, and 4. These results suggest the question why a matrix-interface composite exhibits interaction of interfaces, provided that it is subjected to uniform strain boundary conditions, while interfaces are not interacting under uniform stress boundary conditions. In order to clarify the interaction properties of different interfaces, it is useful to consider the average rules for stresses and strains. For any multiphase composite with volume Ω, which is subjected either to uniform stress or strain boundary conditions, they read as [35] Next, we specify Eqs. (114) and (115) for a matrix-interface composite. In this context, it is important to consider that the interfaces are two dimensional, i.e., interfaces occupy a vanishing volume such that the solid matrix fills the entire volume of the composite: Ω = Ω s . Because finite interface tractions are acting in a vanishing interface volume, the stress average rule (114) simplifies to [36], Note that the integral in (116) refers to the volume of the solid. Equation (116) underlines that the macrostress Σ is equal to the average stress in the solid matrix. The situation is different when it comes to the strain average rule, where the dislocations contribute to the macrostrain [36], In Eq. (117), summation index j refers to the interface family number, j = 1, 2, each one of the interface families comprising N j |Ω| interfaces defined on surface domains C j , and n j is standing for the unit normal to all the interfaces of the j-th family. Equation (117) underlines that the macrostrain is decomposed into (i) the average strain in the solid matrix and (ii) a contribution related to the dislocations of the interfaces. Notably, the average stress of the solid, see the stress average rule (116), and the average strain of the solid, see the strain average rule (117), are related via generalized Hooke's law reading as ⎛ The stress average rule (116) and the strain average rule (117) allow for the sought interpretation of interface interaction, described next. Given the matrix-interface morphology illustrated in Fig. 1, each interface family "feels" the stress and strain states of the surrounding solid matrix, and they drive the evolution of interfacial dislocations and tractions. Accordingly, uniform stress and uniform strain boundary conditions have markedly different effects on the interaction between the interfacial dislocations, on the one hand, and the average stress and strain states of the solid matrix, on the other hand: -Under uniform stress boundary conditions, average stresses and strains in the solid are controlled from outside of the composite; in more detail, the average stress of the solid matrix is equal to the externally controlled loading, see the stress average rule (116). Generalized Hooke's law (118), in turn, underlines that also the average strain in the solid matrix is directly proportional to the externally prescribed loading. This implies that the dislocation histories of both interface families evolve under the externally controlled average stress and strain states of the solid matrix. Hence, the interfaces do not interact, because increasing dislocations do not change the average stress and strain states of the solid matrix, but result only in an increase of the macrostrain, see (117). With increasing time, the dislocations increase, and they take over part of the imposed macrostrain, see (117), such that the average strain of the solid matrix decreases. According to (118), also the average stress level in the solid matrix decreases ("relaxes") proportionally. In other words, evolving dislocations of one specific interface result in a reduction in the average stress and strain states in the solid matrix, and this is "felt" by all other interfaces. This explains the interaction underlined by the coupled relaxation times (103).
We here investigated the interaction between two different interface families, but the conclusions can be extended to consideration of n interface families. The creep functions contain n exponentials with n different characteristic creep times, whereby each interface family influences just one characteristic creep time. Also the relaxation functions contain n exponentials with n different characteristic relaxation times, but all interface families influence each and every characteristic relaxation time. In the herein studied special case of two interface families, the characteristic times involve square-root expressions, stemming from the solution of a second-order polynomial. In the general case of n interface families, the characteristic relaxation times involve nth-order root expressions, stemming from the solution of an nth-order polynomial. This underlines that macroscopic creep tests are conceptually clearly preferable for top-down identification of interfacial properties. , Equation (11) contains influence tensors B Σ 1 and B Σ 2 . Their nonvanishing components read as and B Σ 2,zx x = B Σ 2,xzx = B Σ 2,zyy = B Σ 2,yzy = Equations (12) and (13) and Equations (12) and (13) contain influence tensors D Σ 11 and D Σ 22 . Their nonvanishing components read as