Eigenvector Centrality is a Metric of Elastomer Modulus, Heterogeneity, and Damage

We present an application of eigenvector centrality to encode the connectivity of polymer networks resolved at the micro- and meso-scopic length scales. This method captures the relative importance of different nodes within the network structure and provides a route toward the development of a statistical mechanics model that correlates connectivity with mechanical response. This scheme may be informed by analytical and semi-analytical models for the network structure, or through direct experimental examination. It may be used to predict the reduction in mechanical performance for heterogeneous materials subjected to specific modes of damage. Here, we develop the method and demonstrate that it leads to the prediction of established trends in elastomers. We also apply the model to the case of a self-healing polymer network reported in the literature, extracting insight about the fraction of bonds broken and re-formed during strain and recovery.

Predicting the macroscopic mechanical response of polymeric solids based on the underlying material structure remains an elusive problem, even after many decades have passed since the development of the earliest molecularly-aware theories. Rather than applying a detailed evaluation of structure at either the micro-or meso-scopic length scales, in practice, one often resorts to using toy mechanical models that require many empirical coefficients and are disconnected from the underlying molecular picture. Though much effort has been invested in capturing the rich phenomena at the molecular scale 1-3 , including chain conformational statistics, chemical and physical cross-linking, and entanglements, we still lack a satisfactory means to pass that information up to the macroscopic length scales 4 . This problem becomes more challenging when damage occurs during the use of the material and is especially relevant to the emerging field of self-healing materials [5][6][7] , in which the mechanical characteristics evolve with time as the result of molecular processes.
Here, we propose that developments in graph theory that are finding broad application in the arena of information science 8-10 may provide a systematic means to encode specific structural details and to do so in a manner cognizant of the relative importance of different portions of the material. This latter trait should prove especially useful in evaluating the robustness of a given material in a specific application. Moreover, these developments also provide a means to capture the network details in a reduced order model.
Eigenvector centrality 11 and its variants find use in such diverse information science applications as contemporary search engines 8 and community identification 9,10 . It is especially useful in ranking the importance of each node in a network based not only on its own functionality or degree, but also on the connectivity of its neighbors. To gather this information, one need only find the principal eigenvector of the adjacency matrix for the network. That is, one finds the eigenvector with the largest eigenvalue for the matrix that has non-zero element entries for pairs of nodes that are connected (when represented in the natural basis of node-identity). The relative value of the elements of the eigenvector indicates the relative importance that the corresponding node holds in the network. In the context of academic papers, this may highlight those articles that hold the most impact in a given scientific community. In the context of polymers, we postulate that network importance ranking can inform models of mechanical response. While the application of graph theory ideas to polymer mechanics has a long history 12 , the application of centrality appears not to have been previously explored.
The adjacency matrix Â provides a table denoting which nodes in the network are linked. For example, Â for a linear chain whose nodes are numbered consecutively with node identification number n ranging from 1 to N is given by eigenvector  e 1 holds the importance ranking of each node within the network. It has unit magnitude, but its elements do not sum to one. The eigenvalue λ is the maximum of either the importance-weighted average degree of the network or the square root of the largest value of degree found within the network 13 . Thus, it is a measure of the functionality of the most important nodes and encodes the relative quantities of different kinds of bonds in the network.
We construct a simple statistical mechanics model to exploit these characteristics via the following ansatz. Let us define an importance operator ιˆ that is represented in the n basis as a diagonal matrix with elements that are simply those of  e 1 . Since the elements of  e 1 provide the relative importance of each node in the network, we associate  e 1 with a distribution function that partitions importance within the graph. Note that importance means the relative contribution that a node makes to the network in the sense of its connectedness to other nodes; in short, high importance means a node is connected to many other nodes that are also highly connected. Thus, the trace of ιˆ may be associated with a partition sum ι ≡Ẑ Tr[ ]. We also identify an underlying Hamiltonian operator Ĥ such that = − ΓẐ Tr H [exp( )]. Here, Γ is the inverse of the Lagrange multiplier that enforces average degree in the network and is the analog of the inverse thermal energy. We propose that this multiplier is proportional to the probability for any two nodes to form a link, P. Within this framework, we define a network free energy ≡ − F P Z log( ). Linear chains provide insight into how to link the Young's tensile modulus E to this free energy. Numerical interrogation of a number of linear graphs reveals that the partition sum is well approximated by We also find that the value of the scale parameter may be approximated by ξ ≈ 1.48 × 10 −3 + 1.39 N −1/2 over the range tested herein. If the exponent of −1/2 varies at higher values of N, the prefactor merely changes in the following analysis leading to Eq. 1. Recall that P is the ratio of the number of bonds formed in the network to the number possible. In the linear case, we have P = 2(N − 1)/(N 2 − N). Thus, in the large N limit, Z ∝ N 0.5 , P ∝ 1/N, and F = −P log(Z) decays as log(N)/N. See Fig. 1 (top), which contains a plot of Z versus N 1/2 for the linear chain case. Now, for linear chains the modulus E drops as 1/N 14 . We therefore predict that the modulus is related to the network free energy by Eq. 1, the guiding hypothesis for our analysis.
Random networks (Erdös -Renyi graphs 15 ) provide a more realistic model for real polymeric solids and have been applied in the past to describe both the onset of a percolated structure 16,17 , as well as to capture the thermodynamics of elastomers 12 . Consider a set of N junctions that form links with any other junction with probability P. Here, one does not expect that the node identification number plays a role. Rather,  We assume that the junction density is fixed and that the volume of the solid grows as N. Applying the relationship between E and F above, Eq. 1, yields: is always small, one may further simplify the prediction for E by expanding the logarithm of Z to find that . Thus, the model predicts an N (or size) dependence for E only for very small networks, as expected.
The connection to classical rubber elasticity theory may be made clear by making explicit the connection to cross-linking density. In the classical theory, the modulus varies directly with the number of cross-links per unit volume, X. This may be approximated as X ∝ PN/V ∝ PN/N, where V is the sample volume. Thus, we find that E ∝ P ∝ X. That is, for a fixed network density, the application of eigenvector centrality and the hypothesis that Eq. 1 holds predicts that the modulus varies directly with cross-link density, in accord with classical models 3 .
To illustrate how heterogeneity impacts the modulus and how to encode it in the proposed scheme, consider a modification to the random network constructed as above. Let there be a ninety percent chance that a given link is "weak" in the sense that it contributes less to the network structure, while the remaining ten percent are "strong". This may mimic, for example, a network composed of physical cross-links (weak) and covalent cross-links (strong). We encode this weakness in Â by assigning a value of 0.1 rather than 1.0 to the entries for the weaker edges. Numerically estimating the modulus E confirms that this quantity is linear in P for both cases. However, both the homogeneous and heterogeneous networks produce trends that fall nearly upon one another; F only provides insight into the global organization of the network (i.e. random versus ordered), but does not probe details at the micro-structural level. The eigenvalue λ, presented in Fig. 2 (top), does capture the variations due to heterogeneity; it captures microstructure. Now, λ is a measure of the importance-averaged degree in the network and is, thus, proportional to PN. Therefore, this formalism provides a direct means to pass specific mesoscale structure up to the macro-scale.
Materials often experience changes in their mechanical properties during their lifecycle 18,19 . Changes such as mechanical chain scission 6 , chemical degradation 20,21 , and morphological rearrangement 2 that result in the loss of macroscopic material properties are often grouped under the rubric of "damage". Quantifying and predicting these changes based on the micro-and meso-scopic pictures presents a serious challenge when modeling specific materials. Though an active area of research reported in the literature 22 , we still lack a robust mathematical definition for damage that goes beyond noting changes at the macroscopic length scale. Here, we propose such a metric and illustrate its application.
Defining a new measure of damage follows from a straightforward extension of the two observations above: i) F captures information about how the network was prepared (random, ordered, etc.) and provides a link to E and ii) λ captures the micro-structural details and plays a role in the prefactor, discussed below. For a homogenous network, one may simply define a rate of damage as the change in F with change in P. This will be characteristic of the particular network type and will predict variations in E in the homogeneous case. This, however, fails to capture the micro-structural details that play a role in damage in real, heterogenous materials. Rather, we define two damage rates, one that captures network class and one that captures microstructure. Though F is insensitive to the particular details of how P is changed, λ does vary with P in a manner that depends upon the specific mode of link deletion employed. For example, consider a heterogeneous network like the one described above, beginning with P = 0.5. Figure 2 (bottom) presents two scenarios, the case in which the deletions are carried out randomly and the case in which deletions are performed only on the "strong" links. The latter scenario may correspond to, for example, the case of chemical degredation of the covalent bonds. Clearly, λ falls much more quickly as a function of P in the case where only the strongly contributing bonds are removed.
We therefore propose a new damage model. First, we introduce a function Ψ to complete the equality between E and F. It is dependent on λ/(PN) − Λ in order to preserve the scaling relied upon to arrive at Eq. 1 and a numerical prefactor that absorbs the log(N) term (discussed further below in the context of the model's application). We further assume that Ψ is a power law, giving Eq. 2 as our final model form.
The new parameter Λ is the value of λ/PN upon the formation of a percolated network, while the prefactor exp(b) and exponent m are free fitting parameters that capture the chemical details. Next, we define a global damage rate defined as δ ≡ F dF dP and a microstructure damage rate δλ ≡ λ d dP M . Here, the suffix "M" implies that the links are deleted in a specific fashion or mode (e. g. randomly chosen, only those of a certain class, etc.). Combining these two terms yields a total damage rate measure suitable for predicting the trend in E under a given mode of link removal: Clearly, if Eq. 3 does not vary with λ under the mode of link deletion, then the variation in E resides entirely in the second term. One expects that this is true when bond deletion occurs randomly rather than targeting specific link categories. For example, the slope for the heterogenous case in Fig. 2 (top) is identical to that for the random case in Fig. 2 (bottom) in which all bonds are randomly deleted; one does not expect that the prefactor changes as a function of cross-linking if the relative proportion of the junctions, weak to strong, remains fixed. That is, there is a variation in Ψ only if the ratio of λ/(PN) changes. The random deletion case, in which λ/(PN) remains fixed is illustrated by the "Random" line for one specific ratio of "strong" to "weak" in Fig. 2. The "Strong" line corresponds to how λ varies with P under a non-fixed ratio of linkage types. These slopes correspond to variations of two different functions; that found in the prefactor Ψ and that in δλ, respectively.
The recent study by Bao and coworkers 23 provides an excellent test application for this model, since their polymer structure is well defined, the cross-linking density is controlled, damage is induced, and the nature of the entanglements are well understood for their chemistry. Those investigators prepared a self-healing elastomer composed of chain-extended poly(dimethyl siloxane) (PDMS) linked by moieties (2,6-pyridinedicarboxamide, or H2pdca) that act as ligands in iron complexes. Each of their chain-extended threads had roughly nineteen iron binding sites and eighteen runs of 6,000 g/mol PDMS. To these threads, they added FeCl 3 , which undergoes ligand exchange to form bonds with the polymer threads (via the H2pdca sites) to form a cross-linked network. They demonstrated that the elastomer thus prepared suffered damage upon elongation followed by subsequent healing, supporting the idea that the polymer-iron-polymer bonds were broken and re-formed. The cross-linking density was varied by changing the ratio of the iron-binding sites to the amount of FeCl 3 added to the material.
The essential step in applying the proposed model to this material lies in generating an ensemble of random adjacency matrices that conform to the known facts about the material. We build this ensemble based on a coarse-grained model of the polymer threads: each thread is composed of nineteen iron-binding sites and eighteen potential PDMS entanglement sites, with the two site types alternating along the chain backbone. If we assume a mono-disperse collection of threads, then we know that there are well-defined static covalent bonds in the system. We give those bonds a weighting of 1 in every realization of Â . We also know that there should be on average nine entanglement sites per thread, since the entanglement molecular weight of PDMS is roughly 12,000 g/mol 24 and the total molecular weight of the threads is about 108,000 g/mol. Thus, each entanglement site forms a bond with another such site with probability of a half. Those bonds are assigned a weight of 0.1 in Â , reflecting the physical intuition that entanglement plays a lesser role in determining the modulus of the system. More on this choice follows below. Finally, we also have some knowledge of the probability for forming a polymer-iron-polymer bond for a given stoichiometric ratio of ligand to iron. Roughly speaking, each thread must have about two polymer-iron-polymer bonding sites at the onset of a percolated network. This is rough since some chains may be part of the network but form dangling ends, while other chains that satisfy this criteria may be bound only intramolecularly. Given this, examining Bao's data, we estimate that the FeCl 3 is only 50% efficient in forming cross-links, since one would need ≈ 5% 1 2

19
FeCl 3 upon network formation, while extrapolating the data to zero modulus indicates that the threshold value is closer to 10%. The factor of a half reflects the fact that each polymer-iron-polymer bond has one iron per every two bound polymer sites. That the efficiency is less than one comports with the physical expectation that not every iron finds two polymer ligand sites due to kinetic frustration. Thus, each iron-binding site has a probability P f ≈ 17% of forming a polymer-iron-polymer bond for the scenario in which one FeCl 3 molecule is added for every six iron-binding sites (i.e., 1:6 Fe:Hpdca-PDMS), the lowest ratio reported by Bao. To those sites, we assigned a weight of 0.99 in Â . We choose a representative sample size of sixty chains to generate an ensemble of one thousand independent realizations of Â within this picture, generating average values for F and λ/(PN). This choice, as well as the choices for the weights, set the scale for the two fitting parameters reported below; nevertheless, once these choices are made, quantitative predictions regarding the changes in the numbers of bonds may be calculated. The assignment of the weights for more general problems merits further study. However, for the particular problem and material that we have studied herein, the weights can be arbitrarily set. This is true because we are breaking only one type of bond (polymer-iron-polymer), and those bonds also set the percolation threshold. A straightforward calculation provided in the Appendix demonstrates that changing the weights only impacts the magnitude of the prefactor fitting term in Ψ.
We note that we make no assumptions regarding affine deformation or spatial heterogeneity. Heterogeneity is captured with the numerical evaluation of the eigenvector and eigenvalues. That is, statistical variations across the network are accounted for if they are encoded within the adjacency matrix. To consider something that is strongly heterogeneous spatially, such as a gradient in cross-linking density across a sample, requires modification of the actual numerical simulations of the networks. This can be done by building correlations into the connectivity of the graph. The basic theoretical model, however, is capable of capturing strong heterogeneity.
The details of the algorithm employed herein to both fit the model and apply it to estimate the unknown values of P f for damaged samples are given below. Note that three functions are used to invert the model so that P f may be found for a given value of the modulus and one function is used to extrapolate the value of Λ. These functional forms are motivated by examination of the data rather than any theoretical considerations. The former three are simple functions that fit the numerical data and whose estimates are ultimately checked for consistency in step 6. The latter function is a straightforward linear interpolation. The fit coefficients (a, b, c, and m) are given subscripts indicating in which step they are calculated. Table 1 contains a summary of the dependencies and roles of the major parameters of the theory. The algorithm follows. l og( /( ) ) 3 3 , where 〈F〉 is calculated from 〈Z〉. 4) For a sample with unknown P f , construct an estimate of the sample's PN/λ (note that this quantity rather than its inverse is modeled as a matter of numerical convenience): a) fit 4 4 4 for the full range of values calculated in step 1 and b) find the value of F (and, thus, PN/λ) that minimizes the squared difference between the target modulus E* and the model, using the fit in step 4.a -that is,  Working within the model described above by Eq. 2 and the enumerated algorithm, we find that λ ≈ . − Λ . E P N F 4522 8( /( ) ) 1 3 for this system when the five lowest concentrations of iron are considered. The value of Λ ≈ 1.12 is found by estimating λ/(PN) for the case of the probability of forming polymer-iron-polymer bonds P f = 0.05. Figure 3 (top) plots Bao's data in units of MPa in black circles as a function of this fit. The line is the ideal relationship. Note that the value excluded from the fit, that is, the black circle data point with the highest modulus value, falls well off the trend line. This data point corresponds to Bao's 1:1 Fe:Hpdca-PDMS sample, a stoichiometry that should match one iron to two ligands (if the assumed efficiency described above is correct) and could, in principal, result in all ligands participating in polymer-iron-polymer bonds. Clearly, however, this is physically improbable, since kinetic effects will always prevent the complete formation of all of the possible bonds in the system. We may calculate the true fraction of bonds formed at this saturated concentration by inverting the relationship between the modulus and bonding probability P f . This leads to an estimate of P f ≈ 79% for this iron concentration. Calculating an ensemble of adjacency matrices with that value produces the value along the horizontal axis for the red circle in Fig. 3 (top). Following the same procedure, we estimate the values of P f for the case of the second highest iron concentration considered, the 1:2 Fe:Hpdca-PDMS sample, after damage induced by elongation to 1500% strain. The stress-strain curves immediately following the strain and after a one-hour recovery period were presented by Bao. From those curves, we estimate the corresponding values of the modulus. Given that data, we estimate that P f ≈ 14% immediately after the 1500% strain and P f ≈ 39% after one hour of recovery. Remarkably, given an original P f = 50%, that translates into an 72% loss of iron cross-links after the damaging strain measurement with a subsequent recovery to a 22% loss after one hour. Experimental validation of our analysis and tracking of the recovered bond fractions as a function of time could provide interesting insight into the self-healing process. Figure 3 (bottom) presents all values of the modulus as a function of bonding probability, illustrating the overall dependence of the tensile properties on that parameter as captured by the model presented herein. Finally, we note that we employed an estimate of the efficiency of cross-linking by the FeCl 3 (50%) in the analysis above. Experimental determination of this efficiency would improve the fidelity of our predictions for P f . We varied this value from 45% to 55% and found that both the exponent m and percolation threshold value Λ were insensitive to this choice. However, the prefactor exp(b) obtained values of 3816.6, 4522.8, and 3675.5 for efficiencies of 45%, 50%, and 55%, respectively. Similarly, P f for the 1:1 Fe:Hpdca-PDMS sample varied as 74%, 79%, and 91% for those efficiencies. The same quantity for the 1:2 Fe:Hpdca-PDMS sample immediately following damage presented values of 12%, 14%, and 14%, while the recovered sample displayed values of P f equal