Dynamic behavior of dual cross-linked nanoparticle networks under oscillatory shear

Via computer simulations, we investigate the linear and nonlinear viscoelastic response of polymer grafted nanoparticle networks subject to oscillatory shear at different amplitudes and frequencies. The individual nanoparticles are composed of a rigid spherical core and a corona of grafted polymers that encompass reactive end groups. With the overlap of the coronas on adjacent particles, the reactive end groups form permanent or labile bonds, and thus form a ‘dual cross-linked’ network. The existing labile bonds between particles can break and reform depending on the bond rupture rate, extent of deformation and the frequency of oscillation. We study how the viscoelastic behavior of the material depends on the energy of the labile bonds and identify the network characteristics that give rise to the observed viscoelastic response. We observe that with an increase in labile bond energy, the storage modulus increases while the loss modulus shows a more complex response depending on the labile bond energy. Specifically, in the case of the samples with the weaker labile bonds, the loss modulus increases monotonically, while for the samples with the stronger labile bonds, the loss modulus exhibits a minimum with an increase in frequency. We show that an increase in the storage modulus corresponds to an enhancement in the average number of bonds in the samples and the characteristics of the loss modulus depend on both the bond kinetics and the mobility of the particles in the network. Furthermore, we determine that the effective contribution of the bonds to the storage modulus decreases with increase in strain amplitude. In particular, while bond formation at small amplitude drives an increase in storage modulus, at large amplitudes it promotes clustering and formation of voids leading to strain softening. Our simulations provide a mesoscopic picture of how the nature of labile bonds affects the performance of cross-linked polymer-grafted nanoparticle networks.


Introduction
Polymer-grafted nanoparticle (PGN) networks are a relatively new class of materials [1][2][3][4] composed of nanoparticles that are decorated with long, end-grafted polymers. The free ends of the grafted polymer 'arms' encompass reactive groups, which bind to reactive groups on neighboring PGNs and thereby interconnect the nanoparticles into an extensive network. These nanoparticle networks have the dual advantage of being inorganic-rich and due to the interconnecting polymer chains, being relatively 'fluidic' (encompassing a relatively low T g ). Inorganic-rich materials are vital for their superior opto-electronic and mechanical behavior, while the low T g could permit these materials to be processed at ambient temperatures through the application of external forces [5]. The presence of solvent within the network further facilitates the chains' motion, and thus makes these composites ideal candidates for lowtemperature processing via applied force.
Using computational modeling, our aim is to determine the response of such PGN networks to mechanical deformations and thus, gain insight into factors that can affect the materials' ability to be dynamically reconfigured. We anticipate that the interlinking, swollen polymers impart the composite with sufficient flexibility (fluidity) that the material could be deformed and reconfigured without undergoing catastrophic failure. Via our recently developed computational approach, we have in fact shown that the PGN networks can withstand significant tensile deformation when the particles are cross-linked by both strong bonds, which break irreversibly, and weaker labile bonds, which are readily broken and reformed. We refer to these systems as 'dual cross-linked' PGN networks (see figure 1). Notably, the stronger ('permanent') bonds formed a 'backbone' that provided structural integrity and strength [2][3][4], while the labile bonds enabled broken interconnections to be re-made and hence, permitted reshuffling of the bonds in the network. Acting in concert, these two types of bonding interactions led to hybrid materials that exhibited pronounced self-healing behavior, even after multiple cycles of mechanical deformation [3]. Furthermore, by optimizing the energies of the labile bonds and fraction of permanent bonds, we could significantly improve the material's toughness and the rejuvenation of its structural integrity [3].
This dual cross-linking might provide useful benefits when the PGN networks are subjected not only to tensile, but also more complex forms of deformation. Here, we extend our studies to determine the response of the dual cross-linked PGN networks to oscillatory shear. In particular, we now formulate a procedure for applying a strain-controlled oscillatory shear to a dual cross-linked PGN network and determine the resulting mechanical response at varying strain amplitudes and frequency of oscillation. We specifically calculate both the elasticity and dissipation of energy that characterize the PGN network under the imposed shear. Through these studies, we can relate the inherent interactions between the particles and shear response of the material in both the linear and nonlinear regimes. Moreover, we can predict effective reconfiguration conditions by pinpointing regions in parameter space that permit deformation of the material without disrupting its structural and mechanical integrity. Below, we first provide the details of our multi-component model. Next, we describe the behavior of the system in the linear regime and then analyze the performance of the material in the nonlinear regime.

Methodology
Our system consists of a concentrated solution of PGNs, which are cross-linked by a combination of strong or 'permanent' bonds and labile bonds to form an extended network (see figure 1). The rigid core of each nanoparticle has a radius r 0 . The grafted polymers form a corona of stretchable chains around the nanoparticles; the thickness of this corona is given by = H qr 0 . In the ensuing studies, all the length scales in the system are given with respect to the core radius r 0 and hence, we consider a PGN of core radius unity and a corona thickness q.
The interaction between two PGNs is modeled through a sum of interaction potentials and is given by link . The first term characterizes the repulsive interactions between the coated nanoparticles and is given by [2,6,7]: Figure 1. Hierarchy of interactions in a network of dual cross-linked polymer grafted nanoparticles (PGNs). The grafted polymer chains contain reactive end groups, which interact to form labile or permanent bonds. At the bond-scale a force sensitive energy barrier separates the bound and the free states and thus, the model captures phenomena on the scale of bond formation (and rupture). At the nanoscopic scale, each spherical nanoparticle is represented by a rigid core that encompasses a corona of grafted chains. The overlapping of the coronas on neighboring particles enables the formation of multiple bonds between the nanoparticles. Finally, via the interactions between the PGNs and formation of bonds, the coated nanoparticles form an extended mesoscopic network. Here, f is the number of arms and σ is the range of the potential. The range of this repulsive potential is related to the diameter of the last blob in the Daoud−Cotton model [7,8] and is given by The second term in the potential, U coh , describes the attractive cohesive interaction between the coated particles. This term is chosen to be a constant for small values of R, the separation between the particles, but balances the repulsion at the edges of the corona in order to allow for the overlap between neighboring coronas [2,9]. Hence, this term is written as: where C is an energy scale, and A and B are length scales that determine the respective location and width of the attractive well in the potential. The final term in the potential, U link , describes the attractive interaction between particles that are linked by the bonded polymer arms. The attractive force acting between the two bonded particles is given by the following equation: where N b is the number of bonds formed between the given pair of particles and κ r ( ) is the spring constant, which increases progressively with the chain end-to-end distance = − r R 2 [2]. The stiffening of the chain is described by the following equation obtained for a worm-like chain [10]: In the above equation, L 2 is the contour length of the chain formed by bonding two corona arms of length L.
The value of N b in equation (3) depends on the extent of overlap between the coronas of the nanoparticles and on the kinetics of bond formation and rupture = − r R 2 [2]. The evolution equation for the number of bonds can be written as = − r R 2 [2]: is the probability of contact of two chain ends, N R ( ) max corresponds to the maximum number of bonds that can be formed between two nanoparticles, and k R ( ) r and k R ( ) f are the respective rupture and formation rates at the individual bond level. The probability of contact, P R ( ) c , is determined by integrating the product of the individual free-end distributions in each corona over the volume of overlap of the PGNs = − r R 2 [2]. The parameter N R ( ) max is estimated from purely geometric considerations by making the simplifying assumption that there is negligible distortion of the coronas.
We use the Bell model [11][12][13][14] to describe the rupture and re-formation of bonds between the polymer arms due to thermal fluctuations. In accordance with the Bell model, the rupture rate, k r p l ( , ) , is an exponential function of the force applied to the bond and is given by is the rupture rate in the zero force limit, which depends on the energy barrier U p l 0 ( , ) separating the bound and the free states. For reversible bond formation in the zero force limit, the ratio of the formation to rupture rate is given by [12,13]. In general, the rate constant of bond formation decreases under an applied force [14]. For the sake of simplicity, however, we neglect the dependence of the rate constant of bond formation on force and consider it to be constant, k f l 0 ( ) [2]. Only the labile bonds can reform after breaking; the permanent bonds can only break. In the following section, we refer to the set of bonded arms as a link; in other words, each link in the network encompasses multiple bonds.
The dynamics of the system is assumed to be in the over-damped regime, with the motion of each particle given by the equation where μ is the mobility and F tot is the total force on the polymer grafted particle. The simulations are performed off-lattice. The total force acting on a particle can be written as where F ext is the external force acting on the edge particles of the particle array (see figure 1). This equation is solved numerically in two steps since the polymer spring force (within the expression for F tot ) in the dynamic equation depends on the number of bonds between particles and consequently, on the evolution of the chemical kinetics given by equation (5). In the first step, we determine the number of bonds at any given time, should take discrete integer values. In order to  , with a random number ξ distributed uniformly between 0 and 1. If , then we truncate the result; otherwise, we increment the integer part of the result by 1. In the second step, we use this value for the number of bonds to calculate the spring force (see equation (3)) in the dynamic equation and numerically integrate the resulting equation using a fourth-order Runge-Kutta algorithm with a time step of Δ = − t T 10 2 0 . In the studies described below, we consider a particle of core radius = particles, which are initially arranged in 13 rows, with 14 particles in each row (see figure 1). The initial state of the system is generated using the following five-step procedure. In the first step, the nanoparticles are placed in an array such that their centers are separated by a horizontal spacing of + q 1.8 (1 ) and a vertical spacing of + q 1.62 (1 ), with a horizontal offset position of + q 0.855 (1 ) between adjacent rows. In the second step, we hold the sample in the initial configuration for 4000 T 0 units of time to allow for the formation of the labile bonds. In the third step, we equilibrate the sample for a time 6000 T 0 , during which the stressed labile bonds are broken and new bonds are formed between the adjacent PGNs according to equation (5). Note that the particles rearrange their positions in continuous space because the simulations are performed off-lattice. In the fourth step, we introduce the permanent bonds into the system. In this step, a certain number of already established labile bonds are converted randomly to permanent ones; as a result, each pair of neighboring particles has an average of = P 1 permanent bonds. Finally, in the fifth step, the resultant sample composed of a dual cross-linked network of permanent and labile bonds is equilibrated for 10 4 time steps, with the breakage and formation of the labile bonds being allowed.
In order to impose an oscillatory shear deformation, the coordinates of the particles at the lower edge of the sample are fixed in space, and each particle at the top edge is moved in the X direction at small increments so that the time- . Here, h 0 is the height of the sample, which is kept constant in the course of simulations, γ 0 is the amplitude of the shear and ω is the frequency of the oscillations. To minimize the finite size effects, the shear displacements are also imposed on the particles at the left and right edges of the sample. It is worth recalling that the simulations are performed offlattice, so that the particles in the bulk of the network move in continuous space in the course of the shear deformation. The X component of the force, F x , that acts on the particles at the top layer (shear force) is recorded as a function of time t at various values of ω and γ 0 . We plot Lissajous curves showing the behavior of F d ( ) x x , the force in the X direction as a function of the time-dependent displacement, d x . These curves are generated and analyzed to characterize the elasticity and dissipation in the system as functions of ω. At small shear amplitudes (the linear regime), the elasticity and dissipation are described by the storage modulus, ω ′ G ( ), and the loss modulus ω ″ G ( ), respectively. Namely, ω ′ G ( ) is the Fourier component of the shear force that is in-phase with the shear strain, and is determined from the tilt of the Lissajous curve. In turn, ω ″ G ( ) is the Fourier component of the shear stress that exhibits a phase difference of π /2 from the shear strain, and is determined from the area bounded by the Lissajous curve. At large shear amplitudes (the nonlinear regime), the Lissajous curves provide the measures of elasticity and dissipation that are generalizations of ′ G and ″ G , respectively. Using the approach described above, we study the viscoelastic response of the material to an applied oscillatory shear deformation characterized by γ γ ω = t sin( ) 0 . The shear strain amplitude γ 0 is calculated as the ratio of the shear deformation of the sample to its height in the undeformed state, i.e., in the state achieved in the course of the sample preparation described above. The shear amplitude and angular frequency of oscillations are varied within the respective ranges of

Linear viscoelastic behavior at small shear strain
We first consider the dynamic response of the system to oscillatory shear deformations for the case of a small amplitude of the applied shear, setting γ = 0.09 0 . In figure 2, we plot the shear force and average number of bonds in the network as functions of the shear strain, and show that the response of the network depends on both the energy of the labile bonds and frequency of oscillation. In particular, depending on the values of U l 0 ( ) and ω, the force-strain Lissajous curves form a straight line (see figures 2(a) and (b)) or an ellipse (see figures 2(b)-(d)). The Lissajous curves indicate that the PGN network exhibits linear viscoelastic behavior at γ = 0.09 0 . The tilt of an ellipse (slope of a straight line) is proportional to the storage modulus ′ G , which characterizes the elasticity of the system. The width of an ellipse is proportional to the loss modulus ″ G , which is a measure of the dissipation in the system due to a lag between the applied strain and the resulting force.
At a given frequency, ω, the samples with the stronger labile bonds exhibit a larger tilt of the force-strain Lissajous curves, indicating that they have higher storage moduli (see figures 2(a)-(d)). The latter behavior can be attributed to the increase in the number of bonds with the increase in bond strength (see figures 2(e)-(h)). Furthermore, figure 2 shows that for both the weaker and stronger labile bonds, the number of labile bonds and the storage moduli (tilt of the Lissajous curves) increase with an increase in the frequency of oscillation, ω. Finally, the dissipation in the system is seen to depend on ω, as indicated by the variations in the width of the Lissajous curves in figures 2(a)-(d).
The storage and loss moduli as functions of the oscillation frequency (at the shear amplitude of γ = 0.09 0 ) are shown in figures 3(a) and (b), respectively. To calculate the moduli ′ G and ″ G , we first determine the amplitude of the force, F 0 , and the phase difference, δ, by fitting the shear force F, which is obtained from the strain-controlled simulations at . Then, the moduli are found as: The storage and loss moduli could also be determined directly from the Lissajous curves [15].
The storage modulus of the network shows two distinct features. First, as noted above, ′ G increases with an increase in the energy of the labile bonds U l 0 ( ) , so that it is consistently greater figure 3(a)). Second, ′ G exhibits a transition from a lower plateau to a higher plateau value with an increase in frequency for both the U l 0 ( ) considered here. The increase in U l 0 ( ) shifts the transition to the lower frequencies ( figure 3(a)). On the other hand, the loss modulus ″ G as a function of ω exhibits distinctly different behavior at the weaker and stronger labile bonds (see figure 3(b)). Namely, ″ G increases monotonically at , the loss modulus exhibits a peak (see curve 2 figure 3(b)). The latter peak is located within the frequency range where the storage modulus of It is instructive to plot the moduli ′ G and ″ G as functions of the normalized frequency ω ω / 0 (see figures 3(c) and (d), respectively), where ω π = k 2 r l 0 0 ( ) is the characteristic frequency, which depends on the energy of labile bonds U l 0 ( ) through the rupture rate constant, k r l 0 ( ) . The storage modulus ′ G as a function of ω ω / 0 shows the same behavior for the networks containing the weaker and stronger labile bonds. Namely, ′ G exhibits a lower plateau regime at ω ω < / 1 0 , a higher plateau regime at ω ω > / 100 0 , and a smooth transition between the latter two regimes at Figure 3(c) indicates, therefore, that at low strain amplitudes, the network elasticity (storage modulus) as a function of frequency is controlled by the rupture rate k r l 0 ( ) of the labile bonds. The latter conclusion cannot be made, however, for the frequency-dependent dissipation in the network. Figure 3(d) shows that the loss modulus ″ G as a function of ω ω / 0 behaves quite differently at The reason for the observed behavior of ″ G is associated with dissipation due to the particle-solvent friction, which does not depend on the energy of the labile bonds and is linearly proportional to ω.
The effect of the particle-solvent friction is demonstrated in figure 4, where we plot the storage and loss moduli as functions of ω ω / 0 at three values of the mobility, which is inversely proportional to friction. Figure 4(a) shows that the storage modulus ′ G is unaffected by the variations in the value of the mobility. In contrast, the loss modulus ″ G exhibits notable changes as the mobility is varied. Within the high frequency region (see figure 4(b)), ″ G decreases with an increase in the mobility, i.e., with a decrease in the particle-solvent friction. Namely, the effect of friction on ″ G becomes noticeable at ω > 0.1 for both the weaker and stronger labile bonds. It is worth noting that at 0 at the mobility of μ 2.5 indicates the onset of a maximum (at ω ω ≈ / 6 0 ) that is not observed in the samples with lower mobility.
To gain greater insight into the simulation results, we consider a simple model, which captures all the salient features of the viscoelastic behavior arising from the dynamics of the labile bonds. In our simulations, the nanoparticles in the PGN network form a hexagonal lattice at equilibrium (as observed in recent experimental studies on a comparable system [16]). We note that the hexagonal lattice is a stable configuration and alternate initial configurations of square and face centered square lattices spontaneously transform to the hexagonal lattice during equilibration. Hence, we consider the system of three particles shown in figure 5 as a representative cell of the lattice shown in figure 1, and determine the response of this cell to a strain-controlled oscillatory shear deformation. Specifically, particles 1 and 2 are fixed in space, whereas particle 3 exhibits oscillations in the horizontal direction, as indicated in figure 5. The position of particle 3 in the vertical direction is held constant at Δ = Y 3 /2 0 1/2 0 and the horizontal coordinate changes in time as , where Δ 0 and γ 0 are the equilibrium center-to-center distance between the particles and shear amplitude, respectively. The time-dependent distance between the particles in pairs 1-3 and 2-3 is calculated as , where the signs '+' and '−'correspond to the pairs 1-3 and 2-3, respectively. Within each pair, the dynamics of the rupture and formation of the labile bonds is described by the following master equation [17] for the probability of finding ⩾ n 1 labile bonds p n : Here, k R ( ) r is the rate coefficient for the rupture of labile bonds defined in section 2 and ( ) is the rate coefficient for the bond formation (cf equation (5)). The probability of having no labile bonds between a pair of particles, p 0 , is found from the normalization condition to be = − ∑ = ∞ p p 1 n n 0 1 . It is worth recalling that the inter-particle distance R in equation (8) is a function of time t (as described above).
The system of equations generated from equation (8) was truncated at = n 17 and solved numerically for the two pairs of particles using the Mathematica™ software. The model parameters used in these calculations were the same as used in the simulations of the array of particles. Upon solving equation (8), the value of attractive force between two particles due to bonded polymer arms was calculated according to equation (3)  . The latter equation means that the link consists of one permanent bond and an average number of labile bonds given by = ∑ n np n n . Finally, the values of ′ G and ″ G were determined using the numerical Fourier transform applied to the X component of the total force acting on particle 3.
In figures 6(a) and (b), we compare the behavior of the storage and loss moduli as functions of ω ω / 0 for the three particle model ( figure 5(a)) and for the array of particles . Only a qualitative comparison is possible since, in our simulations, the magnitudes of the calculated modulus is based on the force required to shear the array of particles and, hence, depends on the size of the array. To make these comparisons, the moduli ′ G and ″ G were normalized by G 0 , which corresponds to the plateau value of the storage modulus at the lower frequencies, i.e., . The normalized storage moduli obtained from the three-particle model (solid lines) and from simulations of the array of particles (symbols) are in quantitative agreement at both values of the labile bond energies considered here (see figure 6(a)). Note that the three-particle model takes into account only the processes of rupture and formation of the labile bonds between the individual pairs of PGNs. The latter result (figure 6(a)) indicates, therefore, that at low shear strains, the individual PGN pairs contribute independently to the elastic properties of the material, i.e., there are no collective effects. Figure 6(b) shows the normalized loss moduli for the three-particle model (solid lines) and PGN array (symbols). For the three-particle model, ″ G exhibits a peak at the value of ω ω / 0 where ′ G exhibits an inflection point within the transition zone between the lower and higher plateaus ( figure 6(a)). Such behavior is typical for cross-linked viscoelastic polymeric systems [18]. In networks of cross-linked PGNs, the viscoelastic behavior is caused by the rupture and formation of the labile bonds, and the dissipation is maximal when the characteristic times of shear deformation and of bond rupture are comparable ω ω( / 1) 0 . Notably, for = U k T 37 l B 0 ( ) (the stronger bond), the peak in the loss modulus for the three-particle model coincides with the peak for the PGN array (figure 6(b)), indicating the predictive ability of the simplistic three-particle model. The peak in ″ G at = U k T 33 is not observed for the cross-linked array of PGN because it is masked by dissipation due to the particle-solvent friction. (We will return to this figure to discuss the behavior in the nonlinear regime, which is captured by figures 6(c)-(f).) The three-particle model served to highlight the importance of the breaking and making of labile bonds in the viscoelastic response of the system to the oscillatory shear. As noted with respect to the discussion of figure 6(b), the dissipation in the system (value of the loss modulus) reaches a maximum value when ω ω/ 1 0 . The rupture and formation of the labile bonds do not contribute to the loss modulus at both the lower and higher frequencies (see the solid lines in figure 6(b)), although the variations Δn are drastically different at these  7(a) and (c)). At the lower frequency, the labile bonds break and reform faster than the progression of the shear deformation. Therefore, no hysteresis is observed in the number of bonds in the system, as evident from figure 2(e) (see curve 2), and hence, there is no contribution to dissipation. In contrast, at the higher shear frequency, the rate of rupture of the labile bonds is much lower than the shear rate. The number of labile bonds does not change noticeably in the course of an oscillation cycle (see curve 2 in figure 2(h)), so again there is no contribution to dissipation. The rupture and reformation of the labile bonds contribute to the loss modulus when the bond rupture rate and shear rate are comparable, so that the number of labile bonds exhibits hysteretic behavior (see loops in curve 2 in figure 2(f)).
As demonstrated above, the viscoelastic properties of the PGN networks that are caused by rupture and formation of the labile bonds are governed by the bond rupture rate constant, k r l 0 ( ) , which depends on the bond energy. The viscoelastic behavior of the PGN networks cannot, however, be considered as a result of a single relaxation process characterized by the relaxation time τ~k 1/ r l 0 ( ) . This is evident from figure 8, which shows Cole−Cole plots for the PGN . In a Cole-Cole plot, the properly normalized frequency-dependent loss and storage moduli are plotted against each other. The normalization is performed as indicated in figure 8, where G 0 is the low frequency plateau value of the storage modulus, and ΔG is the difference between the high and low frequency plateau values of ′ G . It is known that if viscoelastic relaxation is characterized by a single relaxation time, then the Cole-Cole plot forms a semi-circle indicated by curve 3 (dashed line) in figure 8 [19,20]. In contrast, figure 8 shows that the Cole-Cole plots for the viscoelastic moduli obtained from both the results of simulations (symbols) and the threeparticle model (solid curves) lie well below the semicircle. The only exceptions are the symbols on the right that correspond to the simulations at the higher frequency of shear oscillations where the particle-solvent friction dominates (see above).
The features of viscoelastic behavior revealed through the Cole-Cole plots (figure 8) are characteristic of rheologically complex materials, for which the stress relaxation exhibits a stretched exponential dependence on time and hence, the spectrum of relaxation times is continuous [19,20]. The viscoelasticity of rheologically complex materials could be described in terms of the customary linear springs and dashpots only if an infinite, hierarchical system of springs and dashpots is considered [19,20]. It is rather remarkable that the simple three-particle model, which takes into account only the breakage and formation of bonds (see figure 5 and equation (8)), is capable of reproducing the characteristic features of the rheologically complex behavior of the PGN network as seen in figure 8. Hence, we can conclude that the relaxation processes in this system are controlled mostly by the kinetics of bond rupture and formation. In order to demonstrate how the bond kinetics results in a spectrum of relaxation times, in figures 9(a)-(c) we plot the evolution of the average number of bonds, n , obtained for the three-particle model through solving equation (8)   Fourier harmonics. As the attractive force between two particles is proportional to the average number of bonds, the stress relaxation exhibits a spectrum of relaxation times, which is a characteristic feature of the rheologically complex behavior.
It is important to note that the average numbers of bonds n 13 and n 23 shown in figures 9(a)-(c) are obtained by solving the master equation, equation (8), and in this way, we account for stochastic effects. The contribution of the stochastic effects can be estimated by calculating the value of relative deviation ε = 〈 〉 〈 〉 − n n ( 1 / 2 ; the stochastic effects can be  9(e)). It is worth noting that the frequency ω ω ≈ / 6.6 0 , at which the stochastic effects are especially important, corresponds to the peak in the loss modulus (see figure 6(b)).

Nonlinear behavior under large amplitude oscillatory shear
The response of the PGN networks to oscillatory shear deformations is nonlinear under sufficiently high shear amplitudes. In the nonlinear regime, the force-strain Lissajous curves deviate from the elliptical shape. Figure 10 shows Lissajous plots obtained from the simulations of the PGN network at shear frequencies of ω ω ≈ . The shape of the Lissajous curves changes dramatically as the shear amplitude γ 0 is increased, and that the shape distortion is more pronounced at the lower relative frequency of oscillation ω ω / 0 . Furthermore, the shear force (marked in gray in figure 10) exhibits noticeable fluctuations, which increase with an increase in the amplitude of shear γ 0 .
The nonlinear response of the material to the large amplitude oscillatory shear (LAOS), , can be analyzed by representing the time dependent shear force F t ( ) as a Fourier series, i.e., Here, ′ G n and ″ G n are the Fourier coefficients, which depend on the shear frequency ω and amplitude γ 0 . It is more useful, however, to consider the force F as a function of the instantaneous values of shear strain γ and shear strain rate γ, and then decompose F into the elastic part γ F ( ) e , which depends only on the strain, and the strain rate-dependent viscous contribution γḞ e v The above decomposition can be performed in a straightforward manner if the elastic and viscous contributions are represented by Chebyshev polynomials of the first kind [21,22]: e n n n 0 : odd  The first harmonic coefficients in the Fourier series in equation (9) provide a measure of the average (over one cycle of oscillation) values of the storage and loss moduli ′ G and ′ G , respectively, which are related to the first Chebyshev coefficients in equations (11) and (12) as follows: To determine the average storage and loss moduli from the results of the simulations, we employ equations (10)- (12), in which the first four terms of the Chebyshev expansion (n = 1, 3, 5 and 7) are used to fit the shear force over 20 consecutive cycles of oscillation. Fitting was performed using the Mathematica™ software. The Lissajous curves obtained by the fitting procedure are shown in figure 10 by the thick colored lines. Finally, the values of ′ G and ″ G are calculated according to equations (14) and (15)  G obtained from the three-particle model is in a good agreement with the data from the simulations of the network (figure 6(c)), and the agreement is fairly good at the larger strain amplitude γ 0 = 0.45 (figure 6(e)). Furthermore, similar to the linear regime, the three-particle calculations reproduce the peak observed in the loss modulus ″ G for the PGN network containing the stronger labile bonds (compare curves 2 in figure 6(b) and in figures 6(d), (f)). Therefore, we can conclude that the viscoelastic relaxation in the LAOS regime is governed by the same processes as in the linear case. We now utilize the three-particle model (figure 5) to discuss the effect of strain softening observed in the LAOS regime (see figure 11(a)). Figure 13 indicates that the strain softening   shear strain in the vicinity of a particle is calculated by averaging the shear strains determined in the two triangular elements to which the particle belongs in the undeformed state. In figure 15(a), particle 1 belongs to the elements (123) and (145)  , where x y { , } n n are the coordinates of the particle n = … n ( 1, 2, , 5). Then, the shear strain localized at the particle 1 is calculated as γ γ . Figure 15 figure 15(b)). In contrast, at the lower frequency of ω ω ≈ / 0.3 0 , the distribution of strain is rather featureless but remarkably wide, revealing the presence of areas where the local shear strain is lower γ ≈ < ( 0.28 0.45) and higher γ ≈ > ( 0.52 0.45) than the imposed strain ( figure 15(b)). The observed non-uniformity of the local strains is evidence of local rearrangements of particles that occur because the labile bonds between the neighboring PGNs break in the course of deformation if the deformation is sufficiently slow (low frequency). If the deformation is fast (high frequency), the breakage of labile bonds during deformation is less prominent so the shear strain is distributed more uniformly within the sample. Note that this is the relative frequency of oscillations ω ω / 0 that controls the effect. It is

Conclusions
To gain further insight into the mechanical behavior of the dual cross-linked PGN networks, we used our computational model to analyze the response of the material to oscillatory shear. Through these simulations, we could capture the elasticity and viscoelastic properties of the system. In particular, the simulations involved the imposition of a strain controlled oscillatory shear deformation on the samples over a range of frequencies that encompassed the characteristic rupture frequency of the labile bonds. The resulting viscoelastic response was examined both in the small amplitude (linear) and large amplitude (nonlinear) limit.
Our findings showed that the ratio between the rate of rupture of labile bonds and the frequency of the imposed shear deformations determined the in-phase (storage modulus) response of the PGN networks. In particular, we observed a smooth transition of the storage modulus from a lower to higher plateau value with an increase in the oscillation frequency relative to the characteristic rupture frequency of the labile bonds. Furthermore, we found that this increase in storage modulus with frequency of oscillations was accompanied by an enhancement in the average number of bonds in the samples. We demonstrated that unlike the storage modulus, the out-of-phase (loss modulus) response depends on both the bond kinetics and the mobility of the particles in the network. Specifically, we identified that the dominant contribution to loss modulus shifted from the hysteresis in the rupture and formation of labile bonds at lower frequencies to dissipation through particle-solvent friction at high frequencies.
The dissimilarity in the nature of these contributions resulted in distinct differences in the loss modulus response with change in labile bond energies. Notably, in the case of the weaker labile bond samples, the loss modulus increased monotonically, while for the stronger labile bond samples, it exhibited a maxima and minima with an increase in frequency.
We found that the qualitative features of the frequency dependence do not change in the nonlinear limit, indicating that the viscoelastic relaxation in the LAOS regime was governed by the same processes as in the linear regime. We did, however, observe that the local heterogeneity in the shear strain and bond distribution with increase in strain amplitude led to a decrease in effective contribution of the bonds to the storage modulus. In particular, while bond formation at small amplitude drove an increase in storage modulus, at large amplitudes it promoted clustering and formation of voids that lead to strain softening.
Interestingly, although the network response was quite complex, its key qualitative features could be captured by a simple model for a representative cell of three particles. It is noteworthy that the storage and loss modulus response from even this simple model exhibited the features characteristic of rheologically complex materials. Specifically, in the regime where the particle-solvent friction does not dominate, we demonstrated through the Cole-Cole plot that the PGN network displayed the existence of a spectrum of relaxation times that is usually attributed to some cooperative processes [19,20]. Our results indicated that in the PGN networks, the spectrum of relaxation times could arise solely due to the stochastic kinetics of bond rupture and formation.
In summary, we have demonstrated that the viscoelastic response of the dual cross-linked PGN networks depends on two distinct aspects of the network: the labile bond kinetics and particle-solvent friction. We found that there was similarity between the linear and nonlinear regimes in the in-phase response of the networks that was determined by the competition between labile bond rupture kinetics and the imposed frequency of oscillation. We also showed that there was a distinction in the out-of-phase response of different labile bond energy networks that was determined by the contribution due to the particle-solvent friction. In the simulations, these observed features of the dynamic response evolved in response to local rearrangements in the network of strongly interacting particles. Hence, even a simple model for three particles could capture the key qualitative features of this evolution due to local rearrangements. Our simulation results provide a mesoscopic picture of how the nature of labile bonds and the interaction between particles affect the dynamic response of cross-linked PGN networks. Thus, our simulation approach can be used to establish design criteria for PGN networks with improved mechanical response to shear.