Understanding the scaling of boson peak through insensitivity of elastic heterogeneity to bending rigidity in polymer glasses

Amorphous materials exhibit peculiar mechanical and vibrational properties, including non-affine elastic responses and excess vibrational states, i.e., the so-called boson peak. For polymer glasses, these properties are considered to be affected by the bending rigidity of the constituent polymer chains. In our recent work [Tomoshige, et al., Sci. Rep. 9 19514 (2019)], we have revealed simple relationships between the variations of vibrational properties and the global elastic properties: the response of the boson peak scales only with that of the global shear modulus. This observation suggests that the spatial heterogeneity of the local shear modulus distribution is insensitive to changes in the bending rigidity. Here, we demonstrate the insensitivity of elastic heterogeneity by directly measuring the local shear modulus distribution. We also study transverse sound wave propagation, which is also shown to scale only with the global shear modulus. Through these analyses, we conclude that the bending rigidity does not alter the spatial heterogeneity of the local shear modulus distribution, which yields vibrational and acoustic properties that are controlled solely by the global shear modulus of a polymer glass.


I. INTRODUCTION
Amorphous materials exhibit anomalous mechanical and vibrational properties that have been studied for many years by experimental, numerically, and theoretical methods. The vibrational and acoustical properties of such materials have been investigated in many experiments using neutron, light, and X-ray scattering, e.g., Refs. [1][2][3][4][5][6][7][8][9]. Using these methods, anomalies in vibrational and acoustic excitations have been detected, including excess vibrational states, the so-called boson peak (BP), and strong damping of sound wave propagation.
To explain these anomalous properties, the heterogeneous elasticity theory was proposed and developed by Schirmacher and co-workers [10][11][12][13] (see also Refs. [14,15] for the theory in the context of the jamming transition and Refs. [16][17][18] for very recent developments). It is now well-established that amorphous materials exhibit spatial heterogeneity in their local elastic modulus distributions, as supported by numerical simulations [19][20][21] and experiments [22,23]. In the theory, elastic moduli heterogeneities are critical in describing anomalies in the vibrational, acoustic, and thermal properties. The theory notably predicts that the BP and the attenuation rate of sound are more significant when moduli distributions are more heterogeneous. This prediction has been tested and justified by numerical simulations [24][25][26][27][28][29][30].
Anomalous behaviours in polymer glasses have also been reported through both experiments [31][32][33][34][35][36][37] and numerical simulations [38][39][40][41][42][43]. In polymer glasses, the bending rigidity of the constituent polymer chains is an important parameter. In our recent work [44], we studied the effects of the bending rigidity on the global elastic moduli (shear modulus G and bulk modulus K) and the vibrational density of states (vDOS) g(ω) using coarse-grained molecular dynamics (MD) simulations. We demonstrated that the variation of the BP simply follows that of global shear modulus G through the Debye frequency ω D . If this simple scaling behaviour is considered in terms of the heterogeneous elasticity theory, we obtain an important implication that the spatial heterogeneity in local modulus distributions is insensitive to changes in the bending rigidity.
In this study, we examine this correlation by directly measuring the degree of elastic heterogeneity with changes in the bending rigidity. We also study transverse acoustic excitations in the polymer glasses by calculating the dynamic structure factor and examine the connection among the sound velocity, attenuation rate, and the simple scaling behaviour of the BP. Thus, we comprehensively discuss that the effects of bending rigidity in polymer glasses on vibrational and acoustic excitations from the perspective of elastic heterogeneities.
The remainder of this paper is organized as follows. Section II describes the MD simulation details used to characterise the elastic heterogeneity and the acoustic excitation. In Section III, the numerical results and discussions are presented. Finally, our conclusions are drawn in Section IV.

A. Simulation model
We performed MD simulations using the Kremer-Grest model [45], which is a coarse-grained bead-spring model of the polymer chain. Each polymer chain comprises L monomer beads of mass m and diameter σ. We studied the case arXiv:2011.04858v4 [cond-mat.soft] 31 May 2021 of 200 chains of L = 50, such that the system contained N = 200 × 50 = 10000 monomer beads in total, in a threedimensional cubic box of volume V under periodic boundary conditions.
In the Kremer-Grest model, three types of inter-particle potentials are utilised. First, the Lennard-Jones (LJ) potential acts between all pairs of monomer beads, where r and ε LJ denote the distance between two monomers and the energy scale of the LJ potential, respectively. The LJ potential is truncated at the cut-off distance of r c = 2.5σ, where the potential and the force (the first derivative of the potential) are shifted to zero continuously [46]. Second, sequential monomer beads along the polymer chain are connected by a finitely extensible nonlinear elastic (FENE) potential: where ε FENE is the energy scale of the FENE potential, and R 0 is the maximum length of the FENE bond. Following Ref. [42], we employ the values of ε FENE = 30ε LJ and R 0 = 1.5σ. Finally, the bending angle θ formed by three consecutive monomer beads along the polymer chain is controlled by where ε bend is the associated bending energy. We set the stabilised angle as θ 0 = 109.5 • [42]. In the present work, we utilise a wide range of ε bend values: ε bend /ε LJ = 10 −1 , 1, 3, 10, 30, 10 2 , 3 × 10 2 , 10 3 , and 3 × 10 3 . We performed the MD simulations using the Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [47]. Hereafter, the length, energy, and time are measured in units of σ, ε LJ , and σ(m/ε LJ ) 1/2 , respectively. The temperature is presented in units of ε LJ /k B , where k B is the Boltzmann constant. We first equilibrated the polymer melt system at a temperature T = 1.0 and polymer bead number densityρ = N/V = 0.85. We then cooled the system down towards T = 0.05 with a constant cooling rate of dT/dt = 10 −4 , under a fixed pressure of p = 0. Finally, the inherent structure at T = 0 is generated using the steepest descent method. In our recent work [44], we reported the dependence of the glass transition temperature T g and the number densityρ at zero temperature on ε bend .

B. Vibrational density of state and boson peak
The vDOS analysis was performed for the configuration at T = 0. By diagonalizing the Hessian matrix, we obtained the eigenvalues λ k (k = 1, 2, · · · , 3N), which provide the eigenfrequencies as ω k = √ λ k . The vDOS is defined as where three zero-frequency modes are omitted. The expression of the Hessian matrix of the polymeric system was given in Ref. [44]. The Debye law predicts the vDOS as Here, the transverse and longitudinal sound velocities, c T and c L , are given by the bulk modulus K, shear modulus G, and the mass density ρ = mρ as c T = G/ρ and c L = (K + 4G/3)/ρ, respectively. The reduced vDOS g(ω)/ω 2 thus characterises the excess vibrational modes exceeding the Debye prediction, i.e., the BP.

C. Global and local shear modulus
The global shear modulus G and bulk modulus K were evaluated from the stress-tensor response to the shear and volume deformations in the "quasi-static" way, respectively, applied to the inherent structure. For perfect crystalline solids, the mechanical equilibrium is maintained during affine deformation. However, the force balance is generally broken down for amorphous solids under applied affine deformations. Thus, further energy minimization causes additional non-affine deformation (relaxation) towards mechanical equilibrium. In other words, G and K are decomposed into G = G A − G NA and K = K A − K NA . Here, M A and M NA denote the affine and non-affine components of elastic moduli, with M = G and K, respectively. Our recent work [44] also reported the ε bend dependence of G and K. In particular, we demonstrated that the bulk modulus K is much larger than the shear modulus G, and thus the shear modulus has important effects on the low-frequency vibrational properties of the polymeric system.
In this study, we further study the local shear modulus. Specifically, we measure the spatial distribution of the local shear modulus G m , by using the numerical procedure of "affine strain approach", given in Ref. [21]. Note that the analysis completely neglects anharmonic effects and provide zerotemperature limit values of elastic heterogeneities. Briefly, we divided the system into 7×7×7 cubic cells and monitored the local shear stress as a function of the applied shear strain in each local cell. The linear dimension of the cell is approximately W ≈ 3σ. Here, the local strain of the small cell is assumed to be given by the global strain applied to the system. The local shear modulus G m of cell m was measured as the slope of the local shear stress versus the shear strain. The expression of the local modulus was also given in Ref. [21]. Finally, we collected the G m values for all the cells to calculate the probability distribution of the local shear modulus P(G m ). Remark that the average and standard deviation of the local shear modulus distribution is insensitive to the cell size W [21].
As in the LJ glass [19,21], we found that P(G m ) is well fitted to the Gaussian where the relative standard deviation δG m /G provides a measure of the spatial heterogeneity in the local shear modulus distribution.

D. Transverse acoustic excitation
The transverse acoustic excitation can be characterised by the (transverse) dynamic structure factor as a function of the wave vector q and frequency ω [24,27,48,49]: where q = |q|, ' * ' indicates complex conjugation, and · · · denotes the ensemble average over the initial time and angular components of q. Here, the transverse current is expressed by: whereq = q/q, and r i and v i (= dr i /dt) represent the position and velocity, respectively, of the monomer bead i. In general, the dynamic structure factor S (q, ω) exhibits two kinds of peaks: the Rayleigh (elastic) peak and the Brillouin (inelastic) peak. The Rayleigh peak is located at ω → 0 and is related to the thermal diffusion, while the Brillouin peak is related to the (transverse) sound-wave propagation. The Brillouin peak in S T (q, ω) can be fitted by the damped harmonic oscillator function [24,27,48,49], which provides information about the propagation frequency Ω T (q) and the attenuation rate Γ T (q) as functions of the wave number q. The sound velocity is then given by c T (q) = Ω T (q)/q. Note that the sound velocity c T (q) converges to the macroscopic value c T = G/ρ in the long-wavelength limit of q → 0. We numerically calculated the dynamic structure factor S T (q, ω) [Eq. (6)] of the inherent structure for each bending energy ε bend from ε bend = 10 −1 to 3 × 10 3 . Note that the thermal fluctuations are imposed at very low temperature T = 0.05, which is small enough that the derived values are consistent with the zero-temperature limit values. The values of Ω T (q) and Γ T (q) were then extracted using Eq. (8).

III. RESULTS AND DISCUSSION
A. Scaling of boson peak by the Debye frequency and Debye level Figure 1(a) plots the reduced vDOS g(ω)/ω 2 , showing the BP beyond the Debey level A D for each ε bend . The BP frequency ω BP is located at ω BP ≈ 2, but it slightly shifts to the higher frequency with increasing the bending rigidity. In addition, the peak height of g(ω)/ω 2 gradually decreases when ε bend is increased. Figure 1(b) shows the reduced vDOS g(ω)/(ω 2 A D ) scaled by the Debye level A D as a function of the frequency ω/ω D scaled by the Debye frequency ω D . This demonstrates the scaling of the BP by the Debye frequency ω D and Debye level A D for various bendig rigidities of the polymer chain. Note that the scaling property of the BP is also shown for shorter polymer chain with the length L = 3 in our previous paper [44].

B. Debye frequency and global shear modulus
We next examine the relationship between the Debye frequency ω D and the shear modulus G, which is plotted in Fig. 2. As demonstrated in Ref. [44], the bulk modulus K is approximately three to four times larger than the shear modulus G. Thus the term c L −3 becomes negligible, and the Debye frequency ω D can be approximated as which is mainly governed by the shear modulus G. Figure 2 directly demonstrates the relationship of ω D ∝ √ G with changes in ε bend . The density ρ is also changed by changing ε bend , but the effect of density on ω D is close to negligible [44].

C. Local shear modulus distribution
As demonstrated in Fig. 1, the reduced vDOS g(ω)/ω 2 in the BP frequency ω BP regime was well scaled by using the Debye frequency ω D and Debye level A D = 3/ω D 3 This suggests that the frequency and intensity of BP are controlled only by the global sear modulus G. In particular, we obtain the relationship of ω BP ∝ ω D ∝ √ G. According to the heterogeneous elasticity theory [10][11][12][13], this observation implies that the degree of the shear modulus heterogeneity δG m /G is invariant with changes in the bending energy ε bend : this implication is confirmed below.
We plot the probability distribution of the local shear modulus G m in Fig. 3(a); this plot follows the Gaussian form of Eq. (5). Figure 3(b) then plots the scaled distribution P(G m )G as a function of the scaled local shear modulus G m /G, demonstrating the data of P(G m )G versus G m /G nicely collapse for different values of ε bend . Because we can transform P(G m ) (Gaussian form) to this collapse indicates that the scaled standard deviation δG m /G remains unchanged for different ε bend values. This is verified by direct demonstration in Fig. 4, where δG m /G is plotted explicitly as a function of ε bend . Therefore, we can conclude that the bending rigidity of the polymer chain does not alter the degree of the shear modulus heterogeneity. This conclusion justifies the theoretical prediction [10][11][12][13] that vibrational excitations including the BP are controlled only by the global elastic modulus under the condition of constant heterogeneities in the moduli distributions.

D. Transverse acoustic excitation and its link with boson peak
We finally study the transverse acoustic excitation in the frequency regime including the BP. The generalised Debye model [24,50] yields the reduced vDOS g(ω)/ω 2 in terms of the propagation frequency Ω T (q) and the attenuation rate Γ T (q), as follows: with Debye wavenumber q D = (6π 2ρ ) 1/3 . This form can be scaled by ω D and A D = 3/ω 3 D as: Thus, the collapse of the reduced vDOSs g(ω)/(ω 2 A D ) for different values of ε bend indicates that c T /ω D and Γ T /ω D are both independent of the bending energy ε bend . In addition, Eq. (8), which is the damped harmonic oscillator function for the dynamic structure factor S T (q, ω), can be scaled by the Debye frequency ω D : which indicates that S T (q, ω)ω D is simply scaled by ω/ω D , when c T /ω D and Γ T /ω D are independent of ε bend . Below we show that these properties of transverse acoustic excitations are true. Figure 5(a) shows the S T (q, ω) for different values of ε bend . The wave number q is set to its lowest value q min = 2π(ρ/N) 1/3 , which ranges from q min = 0.283 (for ε bend = 0.1) to 0.295 (for ε bend = 3 × 10 3 ). The frequency of the Brillouin peak shifts to higher values with increasing ε bend . We then plot S T (q, ω)ω D versus ω/ω D in Fig. 5(b). It is evident that our calculations of S T (q, ω) are in accordance with the predicted scaling description of Eq. (13).
We also show the sound velocity c T and attenuation rate Γ T as functions of the frequency Ω T in Fig. 6. As expected from the scaling property of g(ω), the data of c T and Γ T collapse for different values of ε bend , although small deviations are detected. These collapses are also consistent with the prediction from Eq. (12) and are explained in terms of the shear modulus heterogeneity. The collapses break down in the high frequency regime above the BP frequency, Ω T /ω D 0.2 > ω BP /ω D ≈ 0.1. Because the generalized Debye model does not hold above the BP frequency [24,50], this deviation is not unexpected.

IV. CONCLUSION
In conclusion, we have numerically studied elastic heterogeneities and acoustic excitations in polymer glasses, with particular attention to the effects of the bending rigidity of the constituent polymer chains. Our main finding is that the degree of heterogeneity in the local shear modulus distribution is insensitive to changes in the bending rigidity. According to the heterogeneous elasticity theory, for unchanging elastic heterogeneities, the vibrational and acoustic properties of amorphous materials are controlled only by global elastic moduli. Consistent with this theoretical prediction, we demonstrated that the BP and properties of the transverse acoustic excitations are both simply scaled only by the global shear modulus. The present work therefore clarified remarkably simple material property relationships in polymer glasses. These originate from the invariance of the local elastic heterogeneities over an extremely wide range of bending rigidity values for polymer chains. Our results also provide good demonstrations that verify the heterogeneous elasticity theory [10][11][12][13], which is among the central theories used to describe the mechanical and vibrational properties of amorphous materials.
We note that effects of polymerization on vibrational properties can be scaled by global elastic moduli [33,37]. On the contrary, some experiments demonstrate that the pressureinduced shift of BP cannot be explained by the global elastic moduli [31,32]. From these observations, we speculate that the polymerization effect is insensitive to the elas-tic heterogeneities as is the bending rigidity, whereas the heterogeneities would be altered by the densification. Furthermore, recent MD simulations revealed antiplasticizer additives significantly modify the local elastic constant distribution in glass-forming polymer liquids [51]. It could be interesting to study how boson peak properties change with evolution of elastic heterogeneities during the antiplasticization process.
At the end of this paper, we would discuss the relationship between the structural relaxation time and the elastic properties. Remarkably, numerical work [52] has proposed and demonstrated a scaling relationship between the structural relaxation time τ α and the Debye-Waller factor u 2 for many types of glass-forming systems, including polymer glasses, as τ α ∝ exp a u 2 −1 + b u 2 −2 (where a and b are constants). Because the Debye-Waller factor in the harmonic approximation limit is estimated as u 2 = 3T ∞ 0 g(ω)/ω 2 dω ∝ T ω BP −2 ∝ TG −1 (where ω BP ∝ √ G is applied) [53], we obtain where α, β, α , and β are constants. This is the idea of the shoving model [54,55], which characterises the activation energy in terms of the global shear modulus G. Interestingly, Eq. (14) has been well demonstrated for polymer glasses by MD simulations, where the plateau modulus G p of the stress correlation function was effectively utilized as the shear modulus [56]. Our results suggest an important condition under which Eq. (14) holds. When the spatial heterogeneity in the local shear modulus distribution is unchanged, the excess vibrational excitations, i.e., the BP, are controlled only by the global shear modulus, indicating that the structural relaxation time is also controlled solely by the global shear modulus.