Boson peak, elasticity, and glass transition temperature in polymer glasses: Effects of the rigidity of chain bending

The excess low-frequency vibrational spectrum, called boson peak, and non-affine elastic response are the most important particularities of glasses. Herein, the vibrational and mechanical properties of polymeric glasses are examined by using coarse-grained molecular dynamics simulations, with particular attention to the effects of the bending rigidity of the polymer chains. As the rigidity increases, the system undergoes a glass transition at a higher temperature (under a constant pressure), which decreases the density of the glass phase. The elastic moduli, which are controlled by the decrease of the density and the increase of the rigidity, show a non-monotonic dependence on the rigidity of the polymer chain that arises from the non-affine component. Moreover, a clear boson peak is observed in the vibrational density of states, which depends on the macroscopic shear modulus G. In particular, the boson peak frequency ωBP is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{G}$$\end{document}G. These results provide a positive correlation between the boson peak, shear elasticity, and the glass transition temperature.

1 Division of Chemical Engineering, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka, 560-8531, Japan. 2 Graduate School of Arts and Sciences, The University of Tokyo, Tokyo, 153-8902, Japan. 3 Division of Materials Science, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki, 305-8573, Japan. 4 Elements Strategy Initiative for Catalysts and Batteries, Kyoto University, Katsura, Kyoto, 615-8520, Japan. *email: hideyuki.mizuno@ phys.c.u-tokyo.ac.jp; kk@cheng.es.osaka-u.ac.jp; nobuyuki@cheng.es.osaka-u.ac.jp The vibrational properties and the BP of polymeric glasses have been studied by both of experiments [18][19][20][21][22][23] and MD simulations [58][59][60][61][62] . The effects caused by non-covalent bonds including bending forces and chain length represent an important feature of polymer glasses. Previous experiments 18,19 have investigated the effects of the pressure or densification on the frequency and intensity of the BP in polymeric glasses. It was demonstrated that the evolution of the BP with pressure cannot be scaled by the Debye values (i.e., the Debye frequency and the Debye level). Therefore, the pressure effects cannot be explained only by the variation of macroscopic elasticity. In contrast, another experiment 20 has shown that the polymerization effects on the BP is explained by the change in macroscopic elasticity as the frequency and intensity variations of the BP are both scaled by the Debye values.
In addition, Zaccone et al. have recently performed MD simulations to calculate the vibrational density of states (vDOS) in polymeric glasses by changing the chain length and the rigidity of the chain bending 61 . This work studied the vibrational eigenstates in a wide range of frequencies and the effects of the chain length and bending rigidity on the high-frequency spectra. Furthermore, Giuntoli and Leporini studied the BP of polymeric glasses having chains with highly rigid bonds 62 . It was demonstrated that the BP decouples with macroscopic elasticity and arises from non-bonding interactions only. Although these studies 61,62 have helped understand polymeric glass properties, the effects of bending rigidity and chain length on the low-frequency spectra and BP need to be further studied.
Herein, the vDOS and the elastic moduli of polymeric glasses are analyzed through coarse-grained MD simulations (see Methods). In particular, the connection between the BP and elasticity as well as the glass transition temperature is explored by systematically changing the bending stiffness of short and long polymer chains. The contributions of the present study are given as follows. We demonstrate that polymeric glasses can exhibit extremely-large non-affine elastic response (compared to atomic glasses), whereas the BP is simply scaled by the behavior of macroscopic shear modulus. This behavior of the BP can be explained by the theory of elastic heterogeneities [24][25][26] . Our results indicate that effects of the bending rigidity on the BP are encompassed in change of macroscopic elasticity, which is in contrast to effects of pressure 18,19 , but instead is similar to effects of polymerization 20 . Furthermore, we show the positive correlation among the BP, elasticity, and the glass transition temperature. Finally, we will discuss the relaxation dynamics in the liquid state, in relation to our results of low-frequency vibrational spectra.

Results
Glass transition temperature. When the polymeric system is cooled down from the liquid state under a constant pressure, the volume of the system monotonically decreases with decreasing the temperature. Figure 1a shows the specific volume v as a function of the temperature T for several different bending rigidities ε bend and the chain length = L 50. For each value of ε bend , the slope of the v-T curve clearly presents a discontinuous change at a certain temperature, which is defined as the glass transition temperature T g . Figure 1b (triangles) presents the value of T g as a function of ε bend . As the rigidity increases from ε = 1 bend to 10 3 , T g progressively increases from .  T 0 45 g to 0.75. Below ε = 1 and above ε = 10 3 , the variation of T g is low or even negligible. In addition, Fig. 1b (circles) plots the density ρ = v ( 1/ ) of the system that is quenched down to = T 0. The density decreases from ρ .  1 09 to 0.97 as the rigidity increases from ε = 1 bend to 10 3 . As the chain bending becomes rigid, the glass transition occurs at a higher temperature, and as a result, the density in the glass state becomes lower. The similar observation was obtained by Milkus et al. 61 .
These behaviors of T g and ρ can be understood by studying the microscopic conformation of the polymeric chains. Figure 2a presents the probability distribution of the angle formed by three consecutive beads along the chain, θ P( ), when changing the rigidity ε bend . Two peaks are observed at approximately θ 70 and θ 120 for a low rigidity (ε ≤ 1 bend ). A similar distribution θ P( ) was also reported in ref. 61 . As the rigidity increases, the peak position in θ P( ) shifts towards θ = .  109 5 0 . It is noted that the bending potential θ U ( ) bend in Eq. 5 tends to stabilize the angle θ at θ = .  109 5 0 . In addition, Fig. 2b presents the radius of gyration R g as a function of ε bend . It can be observed that R g increases from .  R 11 5 g to 16.5 with an increasing ε bend . Importantly, these variations of conformation are induced intensively when the rigidity increases from ε = 1 bend to 10 3 , which exactly matches the region where variations of T g and ρ are observed in Fig. 1b. Therefore, it can be concluded that the conformation changes of the polymeric chains control the glass transition temperature and the density. In fact, as the rigidity of the chain bending increases, the angle θ of the polymer chains tends to be stabilized at θ = .  109 5 0 and the radius of inertia increases. As a result, the glass transition occurs at a higher temperature and the lower density (larger volume). At  ε 1 bend , the effect of the bending interaction of Eq. 5 is weak compared to those of the LJ and finitely extensible nonlinear elastic (FENE) components of Eqs. 3 and 4 (see Methods). However, at  ε 10 bend 3 , the opposite phenomenon occurs.
It is noted that the glass transition occurs at a lower temperature for = L 3 than for = L 50, which is consistent with a previous report 64 . Correspondingly, the values of ρ for = L 3 becomes larger than that of = L 50. However, common results were observed between = L 3 and 50 with respect to the dependences on the rigidity ε bend . Specifically, T g and ρ, as well as the conformation of the polymeric chains progressively change when the rigidity increases from ε = 1 bend to 10 3 , which also occurs for = L 50.
Elastic properties. The elastic properties of polymer glasses are studied by changing the strength of bending rigidity. An external strain is applied to the system at = T 0, which enables to measure the corresponding elastic moduli. Specifically, the volume-changing bulk deformation and the volume-conserving shear deformation are applied, which provide the bulk modulus K and the shear modulus G, respectively 63 . Figure 3 presents the values of K and G as functions of ε bend . Disordered systems exhibit large non-affine elastic responses 2 . The elastic moduli, First, the bulk modulus K is analyzed for = L 50 and presented Fig. 3a. The affine component K A decreases from  K 155 A to 130 as ε bend changes from 1 to 10 3 . The reduction of K A is caused by the decrease of the density ρ with the increasing ε bend (see Fig. 1b). In contrast, the non-affine component K NA shows a non-monotonic dependence on the ε bend . In particular, K NA slightly increases from ε = 1 bend to 30, which is induced by the decrease of the density ρ. As ε bend is further increased above ε = 30 bend , K NA decreases. This is because the non-affine relaxation process is constrained due to the large rigidity of ε bend . As a result, the total modulus of = − K K K A N A also presents a non-monotonic behavior, which is demonstrated in Fig. 3a. From ε = 1 bend to 10 2 , K decreases from  K 80 to 60, which is caused by the reduction of K A . Moreover, K increases from  K 60 to 65 above ε = 100 bend , which is caused by the reduction of K NA . Therefore, the ε bend dependence of the bulk modulus K is determined by the competition between the density reduction and the increase in the bending rigidity.
Further, the shear modulus G is analyzed for = L 50 and presented Fig. 3c. We note that a non-monotonic behavior of the shear modulus as a function of bending stiffness is observed in the previous study 60 . Figure 3c demonstrates that the bending rigidity strongly affects the shear modulus compared to the bulk modulus. Particularly, above ε = 10 bend 2 , both of the affine G A and non-affine G NA components considerably increase. As the shear deformation is anisotropic and causes deformations of the angles θ of polymeric chains, its response is expected to be highly affected by the bending rigidity. Interestingly, contrary to the important increases of G A and G NA , the total shear modulus = − G G G A N A shows a low variation (by comparing   G G 900 ). The bending rigidity increases the affine shear modulus but, at the same time, the non-affine component also increases to cancel the increase in G A , and as a result, the total shear modulus presents a low www.nature.com/scientificreports www.nature.com/scientificreports/ increase. The elasticity of the shear deformation is therefore different from that of the bulk deformation, which is obvious when the elastic moduli are decomposed into affine and non-affine components. Figure 3 also shows K in (b) and G in (d) for = L 3. The values of K and G of = L 3 are smaller than those of = L 50, due to the bonding energy, ε FENE , connecting the monomers along the polymeric chains. The responses of K and G to the variation of ε bend are also weaker for = L 3. However, K and G, as well as affine K A and G A and non-affine K NA and G NA , exhibit overall common dependences on ε bend between = L 3 and 50. Therefore, the decrease in ρ and increase in ε bend engenders similar effects on the elasticity for = L 3 and 50. Finally, it is remarked that the polymer glasses present larger non-affine elastic components than the atomic (LJ) glasses 63,68 . Even under an isotropic bulk deformation, the non-affine K NA (80 for = L 50 and 50 for = L 3, at ε ≤ 1 bend ) is approximately half of the magnitude of the affine K A (155 for = L 50 and 120 for = L 3, at ε ≤ 1 bend ). This result is different from that of the LJ glasses, where a negligible value of .
) was obtained 63 . Larger non-affine moduli reflect various elastic responses due to the multiple degrees of conformations in polymeric chains. Therefore, the non-affine deformation process must be considered to characterize the elastic property of polymeric systems. Interestingly, ref. 69 has reported that non-affine displacements also play an important role on melting of amorphous polymers.
Low-frequency vibrational spectra. Reduced vDOS. Finally, the spectra of vibrational eigenmodes in polymer glasses are studied. The vibrational mode analysis is performed on the configuration of the polymeric system at = T 0, which corresponds to the inherent structure 71,72 . The Hessian matrix is diagonalized to obtain the eigenfrequencies ω k that corresponds to the square root of the eigenvalues λ k , i.e., ω λ = k k ( = . . . k N 1, 2, , 3 p ). The specific expression of the Hessian matrix is given in Supplementary Material.
The statistics of the eigenfrequency provide the vDOS, ω g( ). Figure 4 presents the reduced version of the vDOS, ω ω g( )/ 2 , when changing the rigidity ε bend and for = L 50 in (a) and = L 3 in (b). The reduced vDOS, www.nature.com/scientificreports www.nature.com/scientificreports/ ω ω g( )/ 2 , of the Debye theory is the so-called Debye level A D 71,72 . A D is calculated from the elastic moduli, K and G, as follows: are the longitudinal and transverse sound speeds, respectively. Figure 5 presents the values of ω D and A D as functions of ε bend . As the bulk modulus is approximately four times larger than the shear modulus, ω D and A D are mostly determined with the shear modulus, i.e, ω π ρ ≈ c (9 ) Fig. 4, the polymer glasses present clear excess peaks over the Debye level, i.e., the BP. The BP frequency, ω BP , is defined as the frequency at which ω ω g( )/ 2 is maximal. As ε bend increases, ω BP shifts to a higher frequency. In addition, the height of the reduced vDOS, ω ω g( )/ BP BP 2 , becomes lower. These behaviors of BP are consistent with the previous observation 61 . These shifts are observed in the region from ε = 10 bend to 10 3 for = L 50 and 3. Importantly, this region corresponds to the shear modulus G variations, as shown in Fig. 3c,d. As the bulk modulus is much larger than the shear modulus, the bulk modulus should only have minor effects on the low-frequency spectra. Therefore, the BP of the proposed system should only be controlled by the shear elasticity.

. As shown in
To confirm this hypothesis, the scaled vDOS ω ω g A ( )/( ) 2 D is plotted as a function of the scaled frequency ω ω / D and presented in Fig. 6. As discussed above, A D and ω D are determined mostly by the shear modulus G. Although deviations are observed, the scaled vDOSs collapse for different values of ε bend . This result indicates that the effects engendered by the bending rigidity on the low-frequency spectra are comprised of the the shear modulus changes. A same collapse was observed in effects of pressure on the BP in the covalent-bonding network glass (Na 2 FeSi 3 O 8 ) 6 . In addition, a previous experiment 20 demonstrated that the effects of the polymerization are also comprised by the macroscopic elasticity changes. The collapsed results for (a) = L 50 and (b) = L 3 are consistent with the experimental observation.
The collapses observed in Fig. 6 indicates that ω ω / BP D does not depend on ε bend . As stated previously, when . As ρ varies in a range of 15%, as shown in Fig. 1, the effect of ρ on ω D is weak. Thus, ω D is approximately proportional to G, which leads to ω ∝ G BP in the variation of ε bend . The ε bend dependence of the BP frequency is determined by the shear modulus, which is a macroscopic quantity describing the entire system in an averaged manner. According to the heterogeneous elasticity theory 24-26 , the spatial www.nature.com/scientificreports www.nature.com/scientificreports/ fluctuations of the local shear modulus δG control nature of the BP. The value of δG is quantified by the standard deviation of probability distribution function of the local shear modulus 63 . The collapse of ω ω g A ( )/( ) 2 D as a function of ω ω / D indicates that the shear modulus fluctuations relative to the macroscopic value, δG G / , are constant for all the cases of different bending rigidities. Therefore, the results of this study can be explained as follows. The increase in bending rigidity does not affect the shear modulus fluctuations (relative to the macroscopic moduli) but only affects the macroscopic shear modulus, which leads to the collapse of the scaled vDOS. Furthermore, it is also noted that the recent theoretical study explains the origin of BP in both crystals and glasses in terms of the competition between elastic phonon propagation and diffusive damping 33 . The model based on the phonon Green's function leads to the result ω ∝ G BP . This necessitates further investigations regarding the effects of ε bend on phonon transport and the phonon's Green function.
Participation ratio. To further study the vibrational eigenstates, the participation ratio P k that measures the extent of localization of the eigenmodes k is calculated as follows 34,35 :   Figure 7 presents the value of P k as a function of the scaled frequency ω ω / D , for different ε bend . It is noted that the presented data are the binned average values. Below the BP frequency ω BP , P k progressively decreases when ω decreases due to the spatially localized vibrations. The low-frequency localization below ω BP has also been observed in multiple glasses 34,35,48,49 . Importantly, P k below ω BP collapses between different values of ε bend . This result indicates that the variations of not only the vDOS and the vibrational states due to ε bend can be characterized by the macroscopic shear modulus changes. However, P k does not collapse above ω BP , as also shown in Fig. 7. This result is attributed to the fact that the high-frequency modes above ω BP reflect microscopic vibrations that cannot be captured by the macroscopic elasticity.
Comparison with LJ glasses. The low-frequency spectra are comparable to that of atomic LJ glasses reported in ref. 70 . As observed in Fig. 4, the height of ω ω g( )/ BP BP 2 of LJ glasses is higher than that of polymer glasses, and ω BP is lower than that of polymer glasses. These observations are different from the study reported in ref. 62 , which demonstrated that the low-frequency spectra of polymer glasses correspond to those atomic LJ glasses. In ref. 62 , the bonded monomers interact via a harmonic potential with a large bonding energy scale of = k 2500. This value is two orders of magnitude larger than ε = 30 FENE , investigated in this study. With respect to the large bonding energy, the rigidity of the polymeric chains has a smaller effect on the low-frequency spectra. Therefore, the low-frequency spectra are mainly determined by the non-bonding LJ interactions, whereas the elasticity is mainly determined mainly by the bonding rigidity. As a results, the BP decouples with the macroscopic elasticity, as demonstrated in the previous study 62 .
In contrast to the the results presented in ref. 62 , the rigidity of the polymeric chains is necessary to determine the elasticity and the low-frequency spectra with respect to the bonding energy scale of ε = 30 FENE . In fact, the In this case, the BP couples with the macroscopic elasticity. However, the plot of the scaled ω ω g A ( )/( ) 2 D as a function of ω ω / D does not collapse between the polymer glasses and LJ glass, as shown in Fig. 6. The height of ω ω g A ( )/( ) 2 D is consistent between the polymer glasses and LJ glass, but ω ω / D of the LJ glass is lower than that of the polymer glasses. This result indicates that vibrational states differences between polymer glasses and LJ glasses cannot be described only by changes in macroscopic elasticity, changes in the local elastic properties should be considered as well 6,18,19,38 .
Here we make a note on the finite system size effects. We consider that the present system size of  N 5000 p is not enough large to sample the low frequency modes below the BP. The previous work on the LJ glass 70 employed the very large system size up to = N 1000000 to study the low frequency regime without any system size effect. In contrast, our results of polymer glass for the region below the BP are affected by finite system size effects. However, the system size of  N 5000 p is enough large to study the vibrational modes in the BP. Indeed, we confirm that our results in the BP regime are not contaminated by the size effects. Thus, as long as we discuss on the BP, we do not need to care for the system size effects.
In addition, the length scale of collective vibrational modes in the BP region is discussed. For atomic LJ glasses, the length scale was evaluated as ξ π ω = c 2 / BP T BP , which corresponds to the size of approximately 23 particle 68 . This length scale diverges near the isostatic point or the marginally stable point, theoretically [27][28][29][30] as well as numerically 41,46,47,73,74 . The present study evaluates the length scale of collective vibrational modes in polymeric glasses as ξ π ω = ≈ c 2 / 12 BP T BP , which corresponds to half of that for LJ glasses. The vibrational modes in the BP region are more localized nature due to the polymerization. Moreover, the value of ξ BP is independent of the bending rigidity ε bend because of ω ω ∝ ∝ c BP D T . In other words, the bending rigidity does not affect the length scale of the collective vibrational motions in the BP region. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The glass transition temperature, elastic properties, and the low-frequency vibrational spectra were studied in polymeric glasses. In particular, the bending energy scale was highly varied for long chains ( = L 50) and short chains ( = L 3). As the system becomes rigid by increasing the bending rigidity, the glass transition occurs at a higher temperature, leading to a lower density in the glass phase. The lowering density directly affects the isotropic bulk deformation, but does not affect the shear elasticity. The shear elasticity is mainly controlled by the bending rigidity. The non-affinity of polymeric glasses is much larger than that of atomic LJ glasses. This is due to the more complex conformational relaxations of the polymeric chains during non-affine deformation. Even under an isotropic elastic deformation, the non-affine relaxation process should be considered to describe the elastic response.
In addition, it is demonstrated that the BP frequency ω BP and its intensity are simply scaled by the Debye frequency ω D and the Debye level A D which are mainly determined by the macroscopic shear modulus G. This result indicates that the BP is controlled by macroscopic shear modulus and that the bending rigidity has a small impact on heterogeneities of local elasticity properties. The effects of the bending rigidity on the BP is similar to that of the polymerization, which has also been explained by macroscopic elasticity changes 20 .
The presented results provide a simple relationship between the BP and the elasticity as well as the glass transition temperature. As the system becomes more rigid by increasing the bending rigidity, the glass transition temperature T g and the shear modulus G are increased. On the contrary, the bulk modulus K decreases due to the decrease in the density ρ caused by the increase in the glass transition temperature T g . However, the BP is mainly determined by the shear modulus G: . Therefore, the glass transition temperature, the shear elasticity, and the BP frequency are positively correlated. A similar relationship between T g and ω BP was observed experimentally in ionic liquids systems 75 and also numerically in LJ glasses 76 . It is noted that the studies of refs. 75  www.nature.com/scientificreports www.nature.com/scientificreports/ provided the relationship of ω ∝ T g BP 2 , but a clear power-law like relationship between T g and ω BP was not observed in polymeric glasses.
Finally, it is worthwhile to discuss the relationship between structural relaxation above the glass transition temperature and the elastic properties. In fact, there have been proposed the shoving model, which characterizes the activation energy of the structural relaxation time τ α in terms of the shear modulus G 77,78 . In the literature, a recent study 79 has demonstrated the scaling relationship between the structural relaxation time τ α and the Debye-Waller factor u 2 as τ ∝ + (where a b , are constants) for multiple glass-forming liquids including polymeric glasses. Here, the Debye-Waller factor in the harmonic approximation 80 is estimated as . It is naturally expected that the relaxation dynamics become drastically slow by increasing the bending rigidity because of the following relationship: where α β α β ′ ′ , , , are constants. This simple relationship demonstrates that the BP below T g and the structural relaxation above T g are well correlated in the polymeric glasses with varying the bending rigidity. Further work is necessary to evaluate its validity by calculating τ α , despite the equation analogous to Eq. (2) has empirically been proposed from MD simulations of polymeric glasses 81 .

Methods
Coarse-grained MD simulations are performed by using the Kremer-Grest model 82 , which treats polymer chains as linear series of monomer beads (particles) of mass m. Each polymer chain is composed of L monomer beads, and two cases are considered in this study: long chain length with = L 50 and short chain length with = L 3. In a three-dimensional cubic simulation box under periodic boundary conditions, = N 5000 p and 4998 is defined as the total number of monomers for = L 50 and = L 3 respectively, which means that the number of polymeric chains is = N L / 100 p for = L 50 and 1666 for = L 3. The polymer chain is modeled by three types of inter-particle potentials as follows. Firstly, all the monomer particles interact via the LJ potential: LJ LJ 12 6 where r is the distance between two monomers, σ is the diameter of monomer, and ε LJ is the energy scale of the LJ potential. The LJ potential is truncated at the cut-off distance of σ = . r 2 5 c , where the potential and the force (first derivative of the potential) are shifted to zero continuously 70 . Throughout this study, the mass, length, and energy scales are measured in units of m, σ, ε LJ , respectively. The temperature is measured by ε k / LJ B (k B is the Boltzmann constant). Secondly, sequential monomer-beads along the polymeric chain are connected by the FENE potential: , according to ref. 61 . Finally, three consecutive monomer beads along the chain interact via the bending potential defined as follows: bend bend 0 where θ is the angle formed by three consecutive beads, and ε bend is the associated energy scale. This potential intends to stabilize the angle θ at θ 0 that we set as θ = .  109 5 0 . Here, the value of ε bend in a wide range from ε = − 10 bend 3 to 10 4 , and the effects of the bending rigidity on the vibrational and mechanical properties of the polymeric system are studied.
MD simulations are performed by using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 83 . The polymeric system is first equilibrated in the melted, liquid state at a temperature = .
T 1 0. Further, the system is cooled down under a fixed pressure condition of = p 0 and with a cooling rate of = − dT dt / 10 4 . During the cooling process, the glass transition occurs at a particular temperature, i.e., the glass transition temperature. After the glass transition, the system is quenched down towards the zero temperature, i.e., = T 0 state.

Data availability
The data supporting the findings of this study are available from the corresponding authors upon reasonable request.