Introduction

As a rapid developing discipline, polymer science has attracted much attention in the past few years, since it provides a powerful tool to study the macromolecules with various structures1,2. Currently, an important issue within the field of polymer science is to unveil how the underlying architecture of macromolecules affects their dynamic behaviors. Among various models, the Generalized Gaussian Scheme (GGS)3, which extends the well-known Rouse model4, has provided useful insights into many static and dynamic properties of the polymers. In particular, it has been successfully applied to a large variety of flexible polymer structures, such as dendrimers5, mesh-like polymers6,7, fractals8,9,10, dendritic11,12, regular hyperbranched structures13,14, scale-free15,16 and small-world17,18 networks and so on. All these examples are important representatives of different classes of hyperbranched macromolecules, which possess macroscopically distinguishable behavior. Experimentally, topological features of the materials are evident, despite the averaging, smoothing out character inherent due to the structural disorder or polydispersity in the relaxation measurements19. While the GGS approach provides a basic means for the relation of the structure to the physical properties of polymers, it neglects some realistic features such as excluded volume and polymer stiffness, which are very important especially for biological macromolecules (proteins, DNA, etc.)20,21,22.

While for analytic theories the inclusion of the excluded volume remains to be a hard task, recently, there have been a considerable number of theoretical investigations on the dynamic properties of semiflexible polymers23,24,25,26,27,28,29. Extending pioneering works for linear chains30,31 and stars32, it was put forward as a framework for the arbitrary semiflexible treelike polymers (STP)24, for the complex semiflexible loop polymers25 and so on. In the theoretical approaches the semiflexiblility was introduced by restricting the orientations of the bonds, which can be monitored through the related stiffness parameters23,24,25,26,27,28,29,30,31,32. Under such additional restrictions, it turns out that these approaches provide a more realistic description of semiflexible polymers23,29. However, the price one has to pay for it is an increase of the complexity for an analytical solution of the problem. Indeed, the dynamical matrix of the semiflexible polymers, which couples the set of equations of motion, contains more non-vanishing elements than that of the fully-flexible polymers. Nevertheless, also in this case, based on the numerical and analytical methods, one can determine the spectra of the dynamical matrix, which allows one to study the dynamic properties of semiflexible polymers33,34.

The main purpose of this article is to uncover the dynamic properties of the regular treelike polymer networks in the framework of the STP model24. For this, we investigate in the STP framework the mechanical relaxation properties for a family of recursive treelike polymer networks which display a small-world behavior35. We first analyze the structure of the corresponding dynamical matrix and derive and analyze its spectra. Based on them we investigate for large polymer structures the mechanical loss moduli [G″(ω)], which allows one to observe the dynamical behavior of polymeric systems at different scales3,36,37. By comparing the [G″(ω)] of the flexible polymer networks with that of the corresponding semiflexible ones, we observe that the stiffness parameter has dramatic influence on the dynamic properties of the polymers, especially in the high-frequency domain. Here, the STP framework allows us to carry on an in-depth analysis of the interplay between the structure and the dynamics of the polymer system. In particular, following the scheme of Refs. 33, 34, we show that the eigenmodes of the fully-flexible recursive treelike polymers keep their structures also in case when one introduces semiflexibility. This allows us to obtain analytic expressions for at least half of the corresponding eigenvalue spectra as a function of the stiffness parameter. It turns out that introducing semiflexibility leads to a reduction of the degeneracy of the spectra.

The report is structured as follows: Section Results introduces the construction of the recursive networks, whose dynamics is studied in the STP-framework. Here we also analyze the corresponding dynamical matrices, their eigenvectors and eigenvalue spectra, which allow us to calculate and to discuss the mechanical relaxation loss moduli G″(ω) for the polymer networks. Section Discussion summarizes our main results and conclusions. Finally, section Methods recalls briefly the tools which use in the report, namely, the STP-model24 and the eigenmodes of the fully-flexible recursive treelike polymers35.

Results

Network construction

We start by introducing a model for a class of treelike networks with exponential growth35, which are constructed in a deterministic iterative way. We denote by Ug (g ≥ 0) the deterministic treelike networks after g iterations. Then the networks can be generated as follows. Initially (g = 0), U0 consists of an isolated node, called the central node. For g = 1, f (f is a natural integer) new nodes are generated connecting the central node to form U1. For g ≥ 1, Ug is obtained from Ug−1 by adding f new nodes to each existing node in Ug−1. Fig. 1 illustrates schematically the construction process of a network for the particular case of f = 3 for the first several iterations.

Figure 1
figure 1

Construction of the network Ug corresponding to f = 3 and g = 3.

Full size image

According to the construction algorithm of Ug, it is easy to see that at each iteration gi (gi ≥ 1) the number of newly generated nodes is . Thus, the network order (i.e., the total number of nodes), Ng, at iteration g leads to

The corresponding number of edges (bonds) at iteration g is Ng − 1, which holds for all treelike networks38.

This class of networks displays typical features of macromolecules in the polymer science5,6,7,8,9,10,11,12,13,14,15,16,17,18. Their cumulative degree distribution Pcum(k), defined as the probability that the degree of a uniformly chosen node is greater than or equal to k, decays exponentially with k as 35. The average path length (APL), which represents the mean of the shortest distance between two nodes over all node pairs, increases logarithmically with the network size35. The diameter, defined as the maximum length of the shortest path between two nodes over all node pairs, also grows logarithmically with the network order35. The features of the APL and of the diameter indicate that this class of polymer networks shows small-world behavior39. In addition, the degree correlations of the networks depend on the functionality f. The network is uncorrelated for f = 3. And the network is assortative and disassortative when the functionality f is in the interval [1,2] and [4, ∞), respectively35.

After introducing the topological characteristics of the polymer networks, next we will investigate the dynamical properties on them, which is the primary topic of the present report.

Dynamics of semiflexible Ug

We start with a summary of the main formulas concerning the dynamics of semiflexible Ug; the details of the STP-model used here are provided in section Methods.

In the STP-framework24, the semiflexibility is modeled through the complementary interactions between the next-nearest neighboring beads. In particular, one introduces the so-called stiffness parameters, which are related to the pairs of adjacent bonds. Exemplarily, the stiffness parameter qi of bead i which connects bonds and is defined as , where denotes an average, l2 is the mean-square length of each bond and the sign depends on the orientation of bonds. We note also that the beads of functionality fi = 1 do not connect any pair of bonds, so that no stiffness parameters are associated with them, see section Methods for details.

The dynamics of the polymer networks is described by a set of Langevin equations3,37, i.e., say, for the x-component of the position vector ri = {x, y, z} one has:

Here τ0 = ζ/K, ζ is the friction constant, K is the spring constant (see also equation (45) in Methods). Moreover, is the x-component of the usual fluctuating Gaussian force acting on the ith bead, for which and hold.

The interaction between the beads are described through the dynamical matrix , whose entries are known in closed form24. Here it is worthwhile to introduce a notation to present the situation of the bead sites, see Fig. 2. Starting from an arbitrary bead i in the structure, we denote its neighbors by ik and the neighbors of ik by iks. The corresponding functionalities of beads i and ik are fi and EquationSource math mrow msub mif mrow msub mii mik respectively, while the corresponding stiffness parameter of junctions i and ik are qi and qk respectively. For these all the non-vanishing elements of the matrix ASTP are given by:

and

whereas all other elements of ASTP vanish24.

Figure 2
figure 2

Schematic drawing of the nearest and next-nearest neighbors of a bead i in a treelike network.

Here we denote one of the nearest neighbors of bead i by ik and one of the next nearest neighbors of bead i by iks.

Full size image

Based on the structure of ASTP, equation (3) ~ (5), it is a simple matter to determine all elements of the matrix for Ug. For simplicity, we assume that all beads of the same functionality fi in the semiflexible polymer Ug have the same stiffness parameter qi. Note that the peripheral beads of the polymer Ug have functionality fi = 1 and hence do not connect any bonds, thus we only have to consider the stiffness parameter of inner beads. In the following, we assume that the stiffness parameter qi of the inner bead i which has functionality fi(fi > 1) is taken to be , where q is a real number between 0 and 1.

Now, we can categorize the diagonal elements into two different situations:

1. If i is a peripheral bead then . It is connected by a single neighbor ik of functionality to the rest of the Ug. Thus, according to equation (3) the value of follows:

2. Otherwise, the bead i has functionality , hence the value of turns out to be:

where the set contains only the neighbor beads {ik} of bead i, which have functionality in the gth generation.

As a second step we consider the non-diagonal nearest-neighboring (NN) elements of . There are also two distinct cases:

1. If either i or ik is a peripheral bead. From equation (40) the value of leads to:

2. Otherwise, both two beads have functionality and , so that turns out to be:

Finally, we have to determine the next nearest-neighboring (NNN) elements of . These elements depend only on the bead ik of functionality lying between the beads i and iks. From equation (41), the can take only one single value, which is

All other non-diagonal elements of the matrix vanish.

According to the construction of Ug and the elements of discussed above, it is easy to observe that the matrix has the following block form:

where each block is a (f + 1)g−1 × (f + 1)g−1 matrix. Here Lg represents the (f + 1)g−1 inner beads. Cg and Dg correspond to the f(f + 1)g−1 peripheral beads. Bg describes the interaction between them.

Furthermore, to bead i, at each subsequent iteration, f new beads will be linked. Hence, the degree of bead i in the generation g evolves as:

Thus, in the generation g, all the inner beads have functionalities larger than one. Therefore, from equation (7) all the diagonal elements of matrix Lg are:

and the NN elements of matrix Lg are:

Moreover, the NNN elements of Lg are determined as:

All other non-diagonal elements of the matrix Lg vanish.

Now, as can be inferred from the Ug construction, the structure of the matrix Bg takes the form:

where each block is a (f + 1)g−2 × (f + 1)g−2 matrix. The component matrices , and are diagonal matrices. The matrix obeys similar relations as matrix Bg.

Furthermore, in equation (11) the block Cg represents the NNN interaction among the peripheral beads of the network. Therefore, its elements are given by equation (10). Moreover, Cg is a diagonal matrix. The main diagonal entries of Cg are given by

and by

where (1 + f)j−2 < i ≤ (1 + f)j−1 and 2 ≤ j ≤ g.

Finally, let us consider the blocks Dg, which are also diagonal matrices. According to equation (6) and the bead degree, the main diagonal entries of Dg can be expressed as:

and

where (1 + f)j−2 < i ≤ (1 + f)j−1 and 2 ≤ j ≤ g.

Based on the structure of dynamical matrix discussed above, in the next subsection, we will determine the eigenvalue spectra of this matrix.

Eigenvalue spectra of

In order to study the dynamical properties of the semiflexible polymer Ug, we will use the results of the previous subsections. In particular, the solution of the set of equations (2) requires diagonalization of . For several highly symmetric structures, such as dendrimers and Vicsek fractals, the set of eigenvalues can be obtained with the help of the similarity between the eigenvectors' structure of the fully flexible polymer and that of the corresponding semiflexible case33,34, fact which also is of judicious use in case of Ug-networks. We develop our study of the eigenvalues of the semi-flexible polymers in three steps: First, in Methods we recall the spectra distribution of the fully flexible case; here we discuss the features of the corresponding eigenvectors and show that these features are also applicable to the semiflexible Ug polymers – we use a combined analytical and numerical approach to derive their eigenvalues.

Let v denote the eigenvector, whose corresponding eigenvalue is λ. v can be expressed as

where all vi of the vector v have same sizes. For fully-flexible Ug (see Methods), we distinguish between the following two cases:

In the first case, the vector v satisfies following relations:

Based on the above two expressions of the eigenvectors, one can readily observe that, in this case, only beads of highest generation can move, while their ascendants(inner beads) are immobile, see Fig. 3. Moreover, the sum of the amplitudes associated with the f mobile descendants, which are generated by arbitrary bead in the (g − 1)th generation, vanishes.

Figure 3
figure 3

Schematic drawing of the beads movement corresponding to the case when the inner beads are immobile.

Here Ug has functionality f = 3 and generation g = 2. The sum of the amplitudes of all the f mobile beads generated by the same ascendant in the second generation vanishes.

Full size image

In the second case, the vector v satisfies following relation:

In this case, all of the f mobile descendants generated by the same bead in the (g − 1)th generation have the same amplitude, see Fig. 4.

Figure 4
figure 4

Schematic drawing of the beads movement corresponding to the case when the inner beads are mobile.

Here Ug has functionality f = 3 and generation g = 2. All the f mobile beads generated by the same ascendant in the second generation have same amplitude.

Full size image

We make a brief summary by noticing that for Ug the eigenmotions of beads can be categorized into two groups: (i) Motions involving immobile inner beads. (ii) Motions involving mobile inner beads. As we proceed to show, this finding allows us to readily study the spectrum of the dynamical matrix in the corresponding semiflexible Ug polymers.

First we show that Ug have similar structure of eigenvectors in the fully-flexible and in semiflexible cases, i.e. we consider the following eigenvalue problem:

where is the dynamical matrix of Ug. Based on equation (11), equation (25) can be expressed as

where vectors vi(1 ≤ i ≤ f + 1) are components of v. Equation (26) leads to the following equations:

When the inner beads are immobile, the eigenvector v satisfies equations (22) and (23). Then the system (27) can be reduced to

Let Eg = Dg − Cg be the matrix that have the same size of Dg and Cg, then Eg is also a diagonal matrix. The system (28) defines an eigenvalue problem which is corresponding to the following matrix:

Obviously, each eigenvalue of matrix Eg is also an eigenvalue of matrix ANg with multiplicity f−1. As a real and diagonal matrix, Eg has exactly (f + 1)g−1 eigenvalues and the eigenvalues of matrix Eg is equivalent to the main diagonal entries of matrix Eg.

Based on equations (17)- (20), the main diagonal entries of matrix Eg can be easily evaluated:

and

where (1 + f)j−2 < i ≤ (1 + f)j−1 and 2 ≤ j ≤ g. From equation (30) and (31), there are g distinct eigenvalues of matrix Eg, namely,

whose multiplicity is 1 in the matrix Eg and f − 1 in the matrix and

where 2 ≤ j ≤ g. The multiplicity of this eigenvalue is (f + 1)j−2f in the matrix Eg, while in the matrix its multiplicity is

We now calculate the total number Nn of eigenvalues for the group that the inner beads are immobile

A special case of equation (32) and (33) is the fully flexible case, for which q vanish. We observe that in this case λ1 → 1 and λj → 1, which corresponds to the eigenvalue set FE1g of the fully flexible case.

When the inner beads are mobile, the eigenvector v satisfies equation (24). Then equation (27) can be reduced to

which is equivalent to determining the eigenvalues of the following matrix:

in which each block is a (f + 1)g−1 × (f + 1)g−1 matrix and ASg is a 2(f +1)g−1 × 2(f + 1)g−1 matrix. The diagonalization of matrices ASg can be performed numerically. As discussed in previous subsection, Lg and Bg are sparse matrices, while Dg and Cg are diagonal matrices, that is to say, the matrix ASg is considerably sparse, which is of great help to the diagonalization procedure.

In summary, the total number of eigenvalues corresponding to the group that inner beads are mobile is:

By summing those from group (i), equation (35), and, from group (ii), equation (38), we have

which indicates that we have found the all eigenvalues of Ug.

In order to demonstrate the influence of stiffness, we plot in Fig. 5 and Fig. 6 the distribution of the eigenvalue spectra for Ug. The left part and the right part of Fig. 5 display the eigenvalues in ascending order for Ug of f = 3 and of f = 4, respectively. From Fig. 5, it is easy to notice that with the increase of the stiffness parameter q, the large eigenvalues of the dynamical matrix grow, while the small ones decrease. By comparing the left part and the right part of Fig. 6, we observe that the number of distinct eigenvalues in the semiflexible Ug is higher than for its flexible counterparts. In case of Ug, the reason for this phenomenon lies in the the dynamical matrix, which has more non-vanishing elements for semiflexible Ug. However, one should note, that semiflexibilty not always leads to a different structure of the eigenvalue spectra. For example, the multiplicities of the eigenvalues for fully-flexible and for semiflexible dendrimers are the same33.

Figure 5
figure 5

Spectrum curves of polymer network Ug of generation g = 6 and functionality f = 3 (a) or f = 4 (b) for different degrees of stiffness q.

Full size image
Figure 6
figure 6

Number of distinct eigenvalues for polymer network Ug of generation g = 6 and functionality f = 3 for different degrees of stiffness q, with q = 0 (a) and q = 0.666 (b).

Full size image

Another feature which becomes apparent in both parts of Fig. 5 is that the eigenvalue , which is located in the in-between part of the spectra, has the highest multiplicity. As the stiffness parameter q decreases from 1 to 0, λg increases gradually and approaches 1 (Note that in the corresponding fully flexible network, this eigenvalue equals to 1 for any functionality). From equation (34) we know that the multiplicity of λg is (f + 1)(g−2)f(f−1), which is increasing with functionality f and generation g. For Ug of any g, λg occupies exactly [f(f − 1)/(f + 1)2] part of the whole spectrum. For example, λg takes 3/8 part for f = 3 and 12/25 part for f = 4, respectively. As we will show, these differences of spectra between the two types of polymer networks Ug will lead to different dynamic behaviors.

To conclude this part, it turns out that the determination of eigenvalues discussed here is based on a judicious use of the typical bead motions in the highly symmetric polymers. This approach can provide analytical expressions of a large part of the spectra and their multiplicities in the semiflexible polymer network with arbitrary functionality and genaration. Instead of the traditional brute-force numerical diagonalization of the complete dynamical matrix, this approach considerably reduces the computational complexity of the eigenvalue problem and may pave the way for exploring other semiflexible polymers with complex topology, especially for symmetric architectures.

Mechanical relaxation

The eigenvalue spectrum {λk} of plays a fundamental role in the static and dynamic properties of polymers3,23. In this report we focus on the mechanical relaxation, which is represented by the complex shear modulus37

where its real part G′ (ω) and imaginary part G″ (ω) are the storage modulus and the loss modulus respectively. The dimensionless storage [G′(ω)] and loss [G″(ω)] moduli are given by3

and

where τ0 is as in in equation (2) and the {λk} are the nonvanishing eigenvalues of the matrix .

Based on the eigenvalues discussed in the previous subsection, we calculate the reduced loss moduli [G″(ω)] of Ug. In the left part and the right part of Fig. 7, we display in double logarithmic scales the [G″(ω)] of Ug of f = 3 and f = 4, respectively, for different values of the stiffness parameter q. Here we keep the generation fixed by taking g = 6 and vary the stiffness parameter by increasing q from the pure Rouse (fully-flexible) case, q = 0, to the semiflexible case q = 0.888. From Fig. 7, it is easy to observe that in both Rouse case and semiflexible case, for very low frequencies ω, we have [G″(ω)] ~ ω1; and that for very high frequencies ω, we have [G″(ω)] ~ ω−1. Note that these well-known universal scaling laws hold for nearly all finite polymer networks3. Hence, the particular structure of a polymer leaves its significant fingerprints in the intermediate region. Here, in the double-logarithmic scale, nearly no straight lines are observable in the intermediate region of [G″(ω)], which means that it appears a nonscaling behavior in this region. In Fig. 7, we can infer from the curve, that differences in the stiffness parameter q affect the intermediate behavior of [G″(ω)] dramatically.

Figure 7
figure 7

Reduced loss moduli [G″(ω)] for Ug of functionality f = 3 (a) and f = 4 (b), based on the spectra of Fig. 6.

Full size image

As it is shown in Fig. 7, the [G″(ω)]-curves show a major peak in the intermediate region. With increasing stiffness parameter q, the [G″(ω)] curve starts to bend downwards in the in-between region. The semiflexibility is reflected in the [G″(ω)] through a local minimum and a second minor peak appearing at intermediate frequencies. The reason for this fact lies in the unique spectra distribution of the polymer network. In Fig. 5, we can see that there is a clear gap in each spectrum near the value 1.0: In the left part of Fig. 5, for q = 0.333, the largest eigenvalue smaller than 1.0 is λ≈0.9808 and the smallest eigenvalue greater than 1.0 is λ≈5.1127, thus this gap is about 4.1319, while it is equal to 7.489 for q = 0.666 and to 20.947 for q = 0.888, respectively. Hence, the gap in the spectra depends strongly on q, growing with q dramatically. Moreover, with increasing stiffness parameter q, the [G″(ω)] curves get wider especially in the high frequency region, which is also observed before in a series of semiflexible polymer networks23,33,34. The explanation for this phenomenon lies in the broadening of the eigenvalue spectra with growing stiffness parameter q.

In order to have a deeper understanding of [G″(ω)] curves, in Fig. 8 we plot the derivative for the [G″(ω)] curves of Fig. 7 as a function of log10(ωτ0/2) in a simple logarithmic scale. As we can see from Fig. 8, the value of α(ω) starts to decrease from 1.0 in the low frequency domain and falls down around ωτ0/2 ≈ 1. In the intermediate region, the value again goes up rapidly and then it descends to the value −1.0 in the high frequency. From the left part of Fig. 8, we can see that there are 3 intersection points between the straight line α(ω) = 0 and the α(ω) curve for the q > 0.729, which indicates that a local minimum of the [G″(ω)] can be observed when q > 0.729. Moreover, the intersection points between these two curves shift to higher ω-region with growing stiffness parameter, which means that, in the [G″(ω)], the local minimum and the second minor peak shift to the high frequency region with increasing q. This feature is also apparent from Fig. 7.

Figure 8
figure 8

The corresponding slopes α(ω) of reduced loss moduli [G″(ω)] for Ug of functionality f = 3 (a) and f = 4 (b), plotted for different degrees of stiffness.

Full size image

We conclude that the mechanical relaxation functions show that the networks Ug belong to the class of polymeric structures whose dynamics does not scale in the intermediate frequency (or time) domain. Such a non-scaling behavior show other well-known structures, namely dendrimers and some structurally disordered SFNs (SDSFNs)12,15,27,33,41. Nevertheless, their relaxation dynamics differs from that of the Ug, especially for the semiflexible case. Namely, the [G″(ω)] of semiflexible dendrimers show a considerable broadening towards both low and high frequencies33, whereas for Ug the curves in the low-frequency domain display only little differences. Hence the relaxation at the large scales is for the Ug networks less affected by stiffness than for dendrimers. In comparison to SDSFNs41, inclusion of stiffness leads for several SDSFNs to a local minimum in [G″(ω)] which, however, considerably less pronounced than by Ug. This is also obvious from the inspection of the corresponding local slopes, which for SDSFNs do not show big jumps as in Fig. 841.

Discussion

In summary, we presented a systematic theoretical investigation of the dynamic properties of a family of growing semiflexible treelike polymer networks in the framework of STP. The main goal of this article is to explore the impact of various polymer structures and of the different degrees of stiffness q, on the dynamic behavior of the polymer networks. To achieve this goal, we analyzed the dynamical matrix for the polymer networks and characterized their spectra. We succeeded to obtain a large part of the spectra analytically.

Based on the eigenvalue spectra, we have investigated the mechanical relaxation forms for the semiflexible polymer networks. In the in-between region, [G″(ω)] curve shows a nonscaling behavior, while in the very low and very high frequency, it shows ω−1 and ω1 behaviors respectively. Moreover, it turns out that the relaxation behavior is very sensitive to the stiffness. Indeed, with the increasing stiffness parameter q, the [G″(ω)] curves get broader especially in the high frequency domain (showing that the stiffness is very important at the local scales) and they start to show a local minimum and another minor peak in the in-between region. Such a distinct qualitative behavior is observed when the stiffness parameter q gets higher than 0.729 (for f = 3), as it follows from the analysis of the derivative for the [G″(ω)]. These observations show that the considered networks Ug belong to the same class as the well-known, extensively synthesized42 dendrimers. As we have found, however, the semiflexibility allows one to distinguish the Ug networks from dendrimers and from other representatives of the common class. We hope that recent progress in the synthesis, in particular, the realization of dendrimers with hypermonomers (monomers with high number of functional groups) and with monomers consisting of both active and inactive groups42, will pave the way to the compounds with the properties of Ug. Finally, it is expected that the methods presented here can be extended to other classes of large semiflexible polymer networks.

Methods

Modelling Semiflexiblity

The polymer structure is modeled as a network, which consists of N beads located at connected by springs (bonds)

We can express the transformation (43) from the bond to the positions' variables by the incidence matrix G38,

Here the GT is the transposed matrix to G; the matrix G = (Gia) has nonzero entries only Gja = −1 and Gia = 1, where the bond a starts in bead j and ends in bead i.

For the fully flexible polymers, the potential energy between beads is purely harmonic, so that it is diagonal in the bonds,

Here the sum runs over all the bonds that build up the polymer. In equation (45), the spring constant K equals 3kBT/l2, where kB denotes the Boltzmann constant, T is the temperature and l2 is the mean-square bond length (all bonds are of the same mean-square length). From equation (45) it follows that in the GGS picture the equilibrium bond-bond correlations evaluated with respect to the Boltzmann distribution exp(−Vfl/kBT) vanish24.

However, for semiflexible polymers, the bonds are correlated, i.e. their orientations are not arbitrary. In order to account for this feature, as discussed in Ref. 24, one may extend equation (45) and use the generalized potential energy Vs:

Under the assumption that follow normal distribution, the average of with respect to the Boltzmann distribution is

In this way one has the relation between the potential energy and the mean scalar products of bonds. For the latter the following traditional23,24,25,26,27,28,29,30,31,32 conditions are taken: First, the mean-squared length is

For adjacent bonds and which are connected by the bead i one has

Here the parameter qi reflects the stiffness degree of the bead i; the plus or minus sign of qi depends on the connection of bonds and , e.g. the plus sign holds for a head to tail arrangement whereas the minus sign appears in the other cases. In a three-dimensional space qi is restricted by qi < 1/(fi − 1), where fi is the functionality of the bead i40. Another limit qi = 0 leads to a fully-flexible model. For non-adjacent bonds and , one has as in the freely rotating chain model,

Here (a, b1, b2, …, bk, c) denotes the unique, shortest path from bond to bond .

Under the conditions of equations (48) – (50), the potential energy of equation (4) has an analytic-closed form24. In order to get the Vs the frame of beads, we substitute equation (44) into equation (49), leading to

where ASTP is the so-called dynamical matrix and is given by ASTP = GWGT. Based on equation (51), one can readily write down the equation of motions (2), where the structure of ASTP is presented in equations (3)~ (5).

Spectra of flexible polymer networks Ug

According to the construction of Ug, the dynamical matrix AGGS of fully-flexible Ug is the connectivity matrix and is given by35

whose characteristic polynomial satisfies following recursive relation35

We use notation FEg to represent the eigenvalues set of matrix . Note that the polymer networks have (f + 1)g nodes, which indicted that there are (f + 1)g eigenvalues in the set FEg. Based on the recursive relation of the characteristic polynomial Q(λ), FEg can be divided into two subset FE1g and FE2g. That is to say, , where FE1g contains eigenvalue 1 with multiplicity (f − 1)(f + 1)g−1. Thus,

FE2g consists the remaining 2(f + 1)g−1 eigenvalues which are determined by the following equation35

where is an arbitrary eigenvalue in the set FEg−1. Note that equation (55) have two roots for each . Thus, each eigenvalue in FEg−1 generates two new eigenvalues in FE2g. And the FEg set can be fully determined by recursively applying above two equations.

Similarly to the eigenvalues, the eigenvectors of can be determined from those of . Let v denote the eigenvector, whose corresponding eigenvalue is λ. v can be expressed as

where all vi of the vector v have same sizes. We can solve equation to determine the vector v. We distinguish two cases: λ ∈ FE1g and λ ∈ FE2g, which will be addressed as follows.

For the first case of λ ∈ FE1g, where all λ = 0, the equation leads to the following equations:

For the second case of λ ∈ FE2g, we can have following relations.

Resolving equations (59) and (60) to find

which shows that all the vi(2 ≤ i ≤ f + 1) have same values and are uniquely determined by the v1. Equation (61) together with equation (55) indicates that vi is an eigenvector of matrix associated with the eigenvalue determined by . Thus, eigenvector v can be expressed as

which implies that v satisfies the following relation: