Polymer glass transition occurs at the marginal rigidity point with connectivity z* = 4

Abstract

We re-examine the physical origin of the polymer glass transition from the point of view of marginal rigidity, which is achieved at a certain average number of mechanically active intermolecular contacts per monomer. In the case of polymer chains in a melt/poor solvent, each monomer has two neighbors bound by covalent bonds and also a number of central-force contacts modelled by the Lennard-Jones (LJ) potential. We find that when the average number of contacts per monomer (covalent and non-covalent) exceeds the critical value z* ≈ 4, the system becomes solid and the dynamics arrested – a state that we declare the glass. Coarse-grained Brownian dynamics simulations show that at sufficient strength of LJ attraction (which effectively represents the depth of quenching, or the quality of solvent) the polymer globule indeed crosses the threshold of z*, and becomes a glass with a finite zero-frequency shear modulus, G ∝ (z − z*). We verify this by showing the distinction between the ‘liquid' polymer droplet at z < z*, which changes shape and adopts the spherical conformation in equilibrium, and the glassy ‘solid’ droplet at z > z*, which retains its shape frozen at the moment of z* crossover. These results provide a robust microscopic criterion to tell the liquid apart from the glass for the linear polymers.

---

1 Introduction

Dynamical arrest in supercooled melts

The phenomenon of a supercooled liquid transition into an amorphous solid (glass) has been studied extensively and many theories have been proposed, starting with the Gibbs–DiMarzio.^1–4 The modern mode-coupling theory (MCT)⁵ associates the dynamical arrest of ergodic liquid into the non-ergodic glass state with the emergence of a non-decaying plateau of the dynamic scattering function evaluated at the nearest-neighbour distance r_min ∼ 1/k_max.⁶ This scattering function, F(k_max,t), remains nonzero at long times for temperatures below a critical temperature. This MCT temperature T_MC is somewhat higher than the calorimetric glass transition temperature, T_g, for example, in confined polymers T_MC = 1.2–1.3 T_g.⁷ Generally, the glass transition T_g somewhat depends on the quantity measured and on the route followed in bringing the system out of equilibrium upon cooling. Nevertheless, if one focuses on the mechanical signature of the glass transition, i.e. the sharp drop of the low-frequency shear modulus G by many orders of magnitude at T ≥ T_g, the latter can be robustly determined for many different materials, and its value does not vary appreciably with the protocol.⁸

Dense supercooled liquids slightly above T_g feature a slow decay of the scattering function at short times due to local rearrangements, known as the β-relaxation, and a second more dramatic decay at longer times when F(k_max,t) finally falls to zero, known as the α-relaxation.^5,6 In glasses well below T_g there is no α-relaxation left, and the F(k_max,t) has a long-time plateau corresponding to the arrested state.^9,10 However, thermally-activated hopping processes may still lead to a further time-decay of this plateau at long-times.¹¹

In classical theories of dynamical arrest there is no way to discriminate the liquid from the solid glass by just looking at a snapshot of the atomic structure. This remains a crucial gap in our understanding of the glass transition, because we cannot explain the emergence of rigidity and the finite zero-frequency elastic shear modulus G from a purely structural point of view. This issue is related to the impossibility of relating the glass transition to any detectable change in the average number of nearest-neighbours. This number is traditionally given by the integral of the static radial distribution function, g(r), from contact up to its first minimum, which defines the first coordination shell.¹² It is well established that the total coordination number so defined is always ≃12, and remains constant from the high-temperature liquid all the way into the low-temperature glassy state.¹² The reason for this lies in the fact that the total coordination number up to the first minimum of g(r) includes many nearest-neighbour particles which are fluctuating, i.e. not in actual contact. Hence, when one integrates g(r) up to the first minimum, the average number of mechanically-active contacts is significantly overestimated. The emergence of rigidity must be associated with only those nearest-neighbours that do not fluctuate and remain in the ‘cage’ for long times. Only these permanent nearest-neighbours are able to transmit stresses and their number does of course change significantly across T_g (as one can appreciate, for example, upon looking at the first peak of the van Hove correlation function, or from MCT calculations¹²). A criterion to identify and estimate the average number z of mechanically-active, long-lived nearest neighbours is therefore much needed to understand the emergence of rigidity at the glass transition.

Vitrification from a solid-state perspective

An alternative way is to consider the glassy solid with the tools of solid state physics, in particular, using nonaffine lattice dynamics. In the presence of structural disorder, a solid lattice deforms under an applied strain very differently from well-ordered centrosymmetric crystals. Additional nonaffine atomic displacements are required to relax the unbalanced nearest-neighbour forces transmitted to each atom (particle) during the deformation. These local forces cancel to zero in centrosymmetric lattice, but are very important in glasses because they cause additional nonaffine displacements and a resulting substantial reduction in the elastic free energy.^13–15

It has been shown that nonaffinity plays a key role in the melting of model amorphous solids. In an earlier theory¹⁴ we have defined the glass transition as a point at which the equilibrium shear modulus G vanishes with increasing the temperature due to the Debye–Gruneisen thermal expansion, linking the average number of contacts per particle z to T via the monomer packing fraction ϕ(T). This criterion can be interpreted as a generalization of the Born melting criterion^16,17 from perfect centrosymmetric crystals to non-centrosymmetric and amorphous lattices. Other recent criteria of glass melting have been also proposed, along the line of the Lindemann condition,^18,19 which put an emphasis on local bonding and atomic-scale dynamics.

The nonaffine linear response theory^13,14,20,21 correctly predicts the Maxwell marginal rigidity criterion at the isostatic point where the total number of constraints zN/2 (with z the mean number of bonds per atom) is exactly equal to the total number of degrees of freedom dN, with central-force interactions in d dimensions (leading to z* = 6 at the isostatic point in 3D). The fixed value z* of the point of marginal rigidity depends crucially on the type of bonding (central-force or bond-bending, or a mixture thereof). The nonaffine theory is also able to recover the correct limit in the case when the interaction is purely covalent bonding (z* = 2.4). In this way, the dynamical arrest is understood in terms of a global rigidity transition, at which the average number of total mechanical contacts on each atom (particle, monomer) z becomes just sufficient to compensate local nonaffine relaxation and causes the emergence of rigidity.

Effects of temperature

When the volume is kept constant, as in canonical computer simulations of bulk polymers, z decreases due to thermal motions out of the cage, and the underlying physics is qualitatively described by the MCT and other localization-based approaches.¹⁹ However, in most practical systems the volume of a liquid or a glass is free to change on changing the temperature, and the basic thermal expansion becomes the leading mechanism for the decrease of z on heating. In real molecular and atomic glasses, the packing fraction is closely and directly related to temperature via thermal expansion, which is a phenomenon determined by the anharmonicity of interparticle interaction.²² The packing fraction ϕ is given typically as ϕ ∝ exp[−α_TT], with the thermal expansion coefficient α_T.

The decrease of ϕ on heating directly corresponds to the decrease in z(T). Hence the equilibrium shear modulus G must vanish at a critical value z(T_g) = z*, which gives the temperature at which the glass ceases to be an elastic solid (which we define as T_g), with a continuous critical-like dependence^14,15. Although there are only a few measurements of equilibrium shear modulus vanishing near T_g in polymer glasses,²³ this picture has been empirically verified for amorphous semiconductors²⁴ (principally alloys of Ge and Se where the stoichiometry controls the connectivity z), where atomic bonding is purely covalent and the glass transition T = T_g happens at z = 2.4 as predicted by constraint-counting.²⁵

Within this picture, the effect of the cooling rate can also be given a microscopic interpretation. If we accept that the average connectivity z is defined by counting mechanically-active, long-lived bonds per particle, it becomes clear that the relevant lifetime of intermolecular interactions between neighbours has to be defined in comparison with the characteristic cooling rate. Thus z becomes an increasing function of the cooling rate, since all contacts that are stable (unbroken) over a characteristic cooling time contribute to z, making it larger at any given T for a faster cooling process. Upon replacing this reformulation of z(T,τ_cool) in the theory for G(T) and solving for T_g by setting G(T_g) = 0, one will obtain that T_g increases upon increasing the cooling rate, in agreement with broad experimental evidence.⁴

Here we consider this problem for the glass transition of a linear polymer chain, a cornerstone problem for soft matter which has not been adequately investigated thus far, see a comprehensive review.²⁶ The dynamical arrest and the freezing of the dynamics in this case is non-standard and, as we will show below, indeed quite different from simple (monoatomic or monomeric) bulk liquids.^5,9,10 Nevertheless, the concept of glass transition as a rigidity transition at the molecular level, associated with nonaffine dynamics summarized above, is in principle applicable to this situation as well. Indeed, we present numerical simulations that confirm the existence of a critical connectivity associated with the dynamical arrest and vitrification of polymer chains.

2 Marginal rigidity condition

The number of bonds per particle requires a careful definition in amorphous systems. In a linear non-branched polymer of N units, each particle has two covalent bonds per atom: z_co = 2(1–1/N), accounting for the chain ends. In addition to these bonds, weaker physical interactions can be established between pairs of monomers in close contact, depending on the local density. All such interactions are of central-force nature, both in vacuo and in solvents (in contrast to covalent bonds which have both central and bending components), and can be modeled by the effective Lennard-Jones (LJ) potential. We shall count a contribution to z_LJ (a physical ‘contact’) when the two monomers are separated by r ≤ r_min of the LJ potential well, see Fig. 1.

Fig. 1 Schematic of the criterion to define the physical contacts: only pairs of particles that lie within the soft repulsive part of the LJ potential contribute to the z_LJ. R_c is the cutoff length and ε the depth for the attractive LJ potential used in simulations.

The Phillips–Thorpe constraint-counting analysis of marginal stability^27,28 gives the fraction of floppy modes in a purely covalent network: , where each z_co-coordinated monomer contributes 2z_co − 3 bending constraints, in addition to ½z_co stretching constraints. The contact number z* occurs when f = 0 and no more ‘floppy’ zero-frequency modes can exist, as was the case at z < z*. This point coincides with the point at which the equilibrium shear modulus becomes non-zero. The vanishing G = 0 at z = z* is a point when the affine contribution to the modulus is exactly balanced by the negative nonaffine contribution,^20,21 producing G ∝ (z − z*). This counting predicts a rigidity transition at z* = 2.4 in a purely covalently bonded network,^28,29 a result independently confirmed by the nonaffine model of linear elasticity.³⁰ In a purely central-force network with no bending constraints, , recovering the Maxwell value of isostaticity in central-force 3D structures and sphere packings: z* = 6. Clearly, the additional bond-bending constraints intrinsic to covalent bonding greatly extend the stability of lattices down to very low connectivity.

When not all bonds are covalent, and some are purely central-force, we need to add the corresponding LJ contacts, to the counting of floppy modes for stretching, but not for bending constraints:

The rigidity transition occurs at z_LJ* = 12 − 5z_co. In our case, when z_co is fixed by the linear polymer chemistry, the critical value of connectivity at which the rigidity is lost is z* = 12 − 4z_co = 4 + 8/N. For very long chains, N ≫ 1, the polymer solidifies into glass when z* ≈ 4, i.e. when each monomer acquires additional z_LJ* ≈ 2 physical bonds. Interestingly, the same condition (z* ≈ 4) has been observed in an experimental study of colloidal gel in high strain-rate flows.³¹

Only the physical central-force contacts contributing to z_LJ are changing upon increasing the packing fraction by δϕ, which is what happens on temperature change is a system without an artificially imposed fixed-volume constraint. Therefore, the critical volume fraction ϕ*, corresponding to z*, will be lower when covalent bonds are present. One may account for this decrease in the simplest meaningful way: ϕ* = ϕ_c − Λ·z_co, where ϕ_c is the packing fraction of non-covalently bonded particles (e.g. in a system of frictionless spheres ϕ_c ≃ 0.64 at random close packing). As a result, one finds the glass transition temperature T_g, dependent only on one free parameter Λ (given that the thermal expansion coefficient is experimentally measured in the glass, e.g. α_T = 2 × 10⁻⁴ K⁻¹ for polystyrene²³):

Importantly, in polymer glasses there are effectively two independent measurements that determine the fitting parameter Λ: one from the value of T_g, the other from the specific dependence of T_g on the degree of polymerisation N. This characteristic dependence of polymer T_g on the chain length N has a name of Fox–Flory equation,²⁶ and has been empirically seen since a long time ago, not only in synthetic polymers,^32,33 but also in glassy biopolymers.³⁴ Very accurate matching values are obtained in this way, e.g. for polystyrene²³ the fitting gives: Λ ≈ 0.1, and ϕ(T) = ϕ₀exp[−α_TT] with ϕ₀ ≈ 0.61.

3 MD simulation of polymer globule

We used the Brownian dynamics simulation package LAMMPS,³⁵ where we could control the strength of physical (LJ) interactions between particles on a polymer chain. We took the chain composed of N = 1000 connected monomeric units consisting of monomers – a length sufficient to not only demonstrate the stiffness-dependent dynamics of individual polymer chains during coil-globule transition, but also to demonstrate the differences between the collapse of chains at different inter-monomer attraction strengths (we have separately verified that the results were nearly identical for N = 2000, meaning that surface to volume ratio of the collapsed globule is no longer relevant at this size³⁶). The simulation is based on the classical coarse-grained polymer model of Kremer and Grest,³⁷ where each monomer has a diameter of σ and the interactions between monomers are described by:

(1) Finitely extensible non-linear elastic (FENE) potential for connected monomers on the chain:

where the maximum bond length R₀ = 1.5σ and the spring constant κ = 30w/σ², with the characteristic energy scale w = 2.5 kJ mol⁻¹ corresponding to the thermostat temperature T = 300 K, as in ref. 37.

(2) The bending elasticity described by a cosine potential of LAMMPS, chosen to produce the chain persistence length l_p = σ (i.e. the case of flexible chain).

(3) The Lennard-Jones potential for non-consecutive monomer interactions, with the strength ε measured in units of w or, equivalently, k_BT. The LJ potential reaches its minimum value of U_LJ = −ε at r_min = 2^1/6σ. For numerical simulations, it is common to use a shifted and truncated form of the potential which is set to zero past a certain separation cutoff. In poor solvent when there is an effective attraction between monomers, the cut-off was set to r_cut-off = 3σ, and the potential is given by:

As the simulation produces a specific configuration of a polymer chain at any given moment of time, we follow the earlier discussion and estimate the mean coordination number z from counting the ‘contacts’ defined as instances when two particles (monomers) have their centers separated by r ≤ r_min. It is obviously not a stringent condition, but we find it most reassuring that when we follow this rule, the z values reproduce the rigidity transition, i.e. when the polymer freezes into a solid glassy state, very close to z* = 4 (see below). However, one has to be open to a minor uncertainty in how we determine z from the simulation data.

Once the decision about the contact radius is made, it is straightforward to create what is called the ‘contact map’ (a concept widely used in protein structure analysis). This map records all instances when, for each chosen monomer, there are particles at or closer than r_min to it, in any given snapshot of fluctuating chain configuration. The total number of all such contact instances is calculated, and given a name Ω. Since the total number of particles within a sphere r_min includes the particle itself, we need to subtract the self-count from Ω. Dividing by N gives the average number of bonds each particle has in this configuration: z = (Ω − N)/N. Note that this is a slightly different way of calculating z, which still produces the correct total number of interacting pairs in the system zN/2.

Identifying the glass transition

Depending on the depth of quenching (effectively measured by the depth of LJ potential well ε), the expanded self-avoiding random walk rapidly collapses into a globule,³⁶ and we monitor the number of contacts between monomers growing as a function of time as the globule becomes increasingly dense. In the expanded coil, only two covalent contacts per particle exist due to the chain connectivity (and all the curves in the plot correctly converge to z ≈ 2). In the dense globular state, the number of physical (LJ) contacts increases. Fig. 2 shows the evolution of this average number of contacts per particle for several values of the quench depth ε, which we measure in units of k_BT, while keeping the average thermostat temperature of the simulation constant. It is clear that at the start of simulation, when the chain is a random expanded coil, z = 2 with good accuracy, and then it increases as the chain collapses into a globule. The line z* = 4 in this plot marks the zone above which the bond-counting theory predicts the polymer globule to be in a solid glassy state.

Fig. 2 The number of contacts per particle for a chain collapsing into a globule in poor solvent, increasing with time measured in simulation timesteps [ts] from an initial expanded coil at t = 0. For a ‘shallow quench’ (ε = 1) the globule must remain fluid, while for a ‘deep quench’ (ε = 10) the glass state sets early in the collapse process, as z exceeds the threshold for marginal stability z* = 4. The numbers on the right show the total number of physical (LJ) bonds in the final globule.

We will now show that the critical connectivity z* = 4 (and the associated transition temperature T_g), at which the glass is predicted by nonaffine dynamics to lose its mechanical stability (or conversely, melt to acquire rigidity), coincides with the point at which the dynamical arrest of the liquid occurs. First consider the shapes of the polymer globule for different quench depths, and compare the globules soon after the collapse (at 1 million time-steps: 1 Mts) and at a very long time after the collapse. Fig. 2 confirms that the local density, or the average number of contacts each particle has, remains constant during this period. The representative snapshots from the simulations are shown in Fig. 3. For liquid globules above the glass transition temperature (ε = 1), the system quickly acquires a spherical shape, and retains it for all subsequent times.

Fig. 3 Snapshots of a globule for different depth of quenching (labelled in the graph), at 1 Mts, and at 12 Mts of simulation. We call an ‘equilibrium’ the state where no further change in its topology has occurred for an order of magnitude longer time than the time it took to form the globule (cf.Fig. 2).

For polymer chains quenched down to T < T_g (i.e. to LJ depths ε > 3), things are very different. In this case we see from Fig. 3 that the globule is not able to equilibrate into the spherical shape (the absolute minimum of free energy). It remains instead frozen in a random shape that it happened to have when the density increased past the rigidity transition at z*, and changes little past that point, clearly remaining far away from thermodynamic equilibrium. The case of ε = 3 represents the borderline situation when the initially non-spherical globule adopts the spherical shape after a very long time.

Fig. 4 illustrates in a more detailed way the observations depicted in Fig. 2, for a polymer collapsing at different depth of quenching. The plot shows how the ‘vitrification time’ τ, measured from the beginning of simulated collapse to the point at which the globule density reaches the level z* = 4 (crossing the dashed line in Fig. 2). Clearly this time is shorter for deeper quenched systems, while the polymer globule with ε = 1 never reaches this threshold. The data in Fig. 4 can be fitted by several functions, including an exponential exp[1/(ε − ε*)], however, the plot displays the best fit with a power-law function τ = τ₀ + A/(ε − ε*)^1.3, with the asymptote τ₀ = 131 ts at very large ε (i.e. T → 0) and the coefficient A = 15 300 (or 15.3 kts). Importantly, it predicts the point of divergence of vitrification time at ε* ≈ 1.25, which corresponds to the effective LJ temperature T_g ≈ 0.8, since our ε is measured in units of k_BT. We should not be treating a precise value of this predicted T_g as accurate: too many uncertainties are present in the simulation, data analysis, and fitting, and to reduce them we would have to carry out a massive statistical averaging over many simulations. But the qualitative magnitude, and the underlying physics of vitrification time after an instant quench are clear from this plot.

Fig. 4 The vitrification time after quenching (in [kts]), at which the z = 4 compaction density is reached during the chain collapse into a globule. The colored data points correspond to the curves shown in Fig. 1. The data is fitted by the power-law curve τ ∝ 1/(ε − ε*)^1.3, with the nominal “glass transition” point ε* ≈ 1.25 labelled on the plot.

4 Analysis of cage dynamics

In order to study the dynamical arrest following the quench in a more quantitative way, let us consider the change (evolution) of LJ contacts (as defined in Fig. 1) with time. First let us define a (somewhat arbitrary) time at which the globule collapse is ‘complete’ and its density (and modulus) reaches the constant plateau value, a time t_ref = 500 kts labeled in Fig. 2. Then let us record each event of a pair of monomers changing their contact (breaking or forming new) that occur after that reference time. Since the density after collapse in all cases remains constant, we have normalized the contact-change number by the total (constant) number of contacts per particle, z_max, on this plateau. Fig. 5 shows the result of this analysis: at short times after t_ref, very few particles change their contact configuration and Δz/z_max increases from zero, apparently in a universal way irrespective of the state of the globule. At very long times, the asymptote value Δz/z_max → 1 implies that all initial contacts are broken and re-formed in a different way.

Fig. 5 The change in the internal contact topology: the difference in physical (LJ) contacts between a reference ‘contact map’ at t_ref = 500 kts and the subsequent long-time globule evolution. The data for liquid globule (ε = 1) is fitted by the t^−0.5 power-law relaxation curve, while the deeply quenched globules develop a long-time plateau after the universal fast relaxation period, indicative of the ‘cage confinement’.

The first point to notice is the behavior of the liquid globule (ε = 1), where we do see that all contacts change over time ergodically; the data is fitted to the power-law 1 − at^−1/2 with good accuracy. The crossover case around ε = 3 is evident here as well, while the glassy globules have a distinctly different evolution of Δz/z_max. At short times, the number of broken contacts initially increases with time in exactly the same way as in the liquid system, clearly reflecting the thermal motion within the cage. After a characteristic time, the change of contacts Δz/z_max is very abruptly arrested and remains frozen for a very long time. At this crossover point there is an almost singular behaviour of the first derivative of Δz(t)/z_max with respect to time, which is qualitatively different from the smooth crossover in simple liquids.⁵ Finally, at very long times, the number of broken contacts starts to slowly increase again, we assume due to isolated thermally-activated hopping events. This process is qualitatively analogous to the one seen in simulations of simple glass-formers,⁶ where it typically appears as the development of a second peak in the self-part of the van Hove correlation function (the first peak in the van Hove function at short distance corresponds to rattling motion in the cage, which we see at short times in Fig. 5). Fig. 6 graphically illustrates this difference in ‘cage’ confinement between the liquid and glassy states of the globule. All particles on the chain are marked as small dots, except the single test monomer (red) and a group of its neighbors (green), which are defined as having a center in a volume of 2 particle diameters around the test particle. These designations are set at t_ref, and then we monitor the labelled particles at later times. In the liquid droplet all particles on the chain are evidently free to diffuse and eventually evenly disperse around the allowed volume. In the glassy droplet, the cage around an arbitrary chosen particle remains essentially intact in time. Even though there are small motions and re-arrangements on a very long time scale (isolated hopping events), the diffusion is clearly arrested.

Fig. 6 An illustration of the ‘cage effect’, comparing the chosen particle and its neighbours at a reference time of t = 500 kts and the subsequent times during the globule evolution. The liquid globule (ε = 1) has the initial neighbours surrounding the particle diffusing around the volume. In contrast, the glassy globule (ε = 10) has the cage of neighbours surrounding the particle intact.

5 Analysis of the glass modulus

On the other side of the marginal rigidity threshold the shear modulus is predicted to increase in proportion to z − z*. As the collapsing polymer globule becomes increasingly dense, at the point of vitrification time τ (Fig. 4) the mean coordination number z(t) crosses the predicted threshold, as we have seen in Fig. 2. In order to estimate this modulus, we may use the expression from linear non-affine theory,²¹ which gives G = (1/30)(N/V)κR²(z − z*), where in our case the bond stiffness is taken as the curvature of the LJ potential in eqn (4) near the minimum at contact (bond) length R ≈ r_min: κ = 36 × 2^2/3ε/σ². The density ϕ = N/V can be roughly estimated by the close-packing density of spheres of size R: . As a result we can plot a predicted evolution (growth) of the emerging glass modulus on the droplet compaction, as shown in the inset of Fig. 7. At each value of effective depth of LJ attraction (which corresponds to a different quench temperature), this modulus eventually reaches a different value of final equilibrium plateau G_eq. (T), at which the compaction finally saturates.

Fig. 7 (a) The inset shows the evolution of shear modulus, predicted based on the relation G ≈ 0.4(ε/σ³)(z − z*) and the data for z(t) from Fig. 1 for a few selected values of quench depth ε. The modulus is measured in units of k_BT/σ³ ≈ 1.5 × 10⁸ Pa for polystyrene at room temperature. (b) The main plot shows the final values of G_eq. for each studied ε, plotted against “LJ temperature” T = 1/ε. Fitting line is eqn (5), producing T_g matching the ε* point in Fig. 3.

The main plot in Fig. 7 presents the values of this equilibrium glass modulus as function of the quench temperature. The temperature dependence is somewhat artificial: it is merely based on the fact that the LJ attractive depth ε is measured in units of k_BT, that is, one could say that the true depth of the van der Waals physical attraction constant varies as inverse temperature (i.e. the LJ temperature T = 1/ε). There is also a lot of uncertainty in extracted values of z, and as a result – in the plateau values G_eq.. Nevertheless, we find a very coherent dependence of the modulus on effective temperature to which the polymer is quenched.

In an earlier paper¹⁴ we have derived the temperature dependence of the glass modulus based on the ideas of thermal expansion coefficient: starting from G ∼ (z − z*), expressing z via the packing density ϕ, and then expressing ϕ(T) via the law of thermal expansion. We have predicted that at the critical point (in the vicinity of T ≤ T_g) the scaling is , which is seemingly not what one finds in the simulation data of Fig. 7. However, we attempted to fit the data with the full theoretical expression for G(T), which we reproduce here:

where ϕ_c is the density at maximum compaction and α_T the thermal expansion coefficient. In the example of polystyrene glass,²³ discussed in Section 2, the product α_TT_g ≈ 0.15. The best fit to eqn (5) in Fig. 7 is achieved with the fitting constant in the exponent equal to 0.55, which is within the right order of magnitude (we could not expect an exact match since our simulation does not use any polystyrene-specific parameters).

It is clear that, especially given the noise in the data points for G_eq., the data is consistent with the theoretical expression, while the square-root critical point behavior is confined to a close vicinity of T_g. It is remarkable and reassuring that the data for the equilibrium modulus (at saturation packing) produces the effective temperature T_g ≈ 0.8, which is the inverse of the point of divergence in Fig. 4, ε* = 1.25, even though these two values arise from very different segments of data and different analyses.

6 Discussion

Note that in the present case of linear polymer chain an additional important source of different qualitative behaviour (besides the mixed covalent/non-covalent character of bonds, which distinguishes the polymer case from the standard liquid vitrification) is the effects on the surface of the polymer globule, where monomers are locally less connected than in the interior. This effect is also present in the standard liquid glass transition; it is experimentally well known that polymer T_g changes when surface to volume ratio becomes very large (e.g. in thin films,^38,39 strings,⁴⁰ or small micelles⁴¹). We have not explored this dependence of the glass transition on the size of polymer globule: we are reasonably assured that the surface effects are not critical in our (apparently – sufficiently large) globule, since we found no significant difference in our numbers for N = 1000 or N = 2000 chains. In any case, this is an important step towards understanding dynamical arrest in polymers at the microscopic monomer level, from the point of view of packing density. This approach is currently missing, since most studies focus on larger scale and cooperative, macroscopic observables such as the viscosity.⁴²

There is a quantitative link between the number of changed contacts that we monitor in Fig. 5 and the standard correlation functions of liquid dynamics. Assuming isotropicity of monomer distribution in the dense globule, the (total) van Hove correlation function gives the probability of finding a second particle at a distance r and at time t from a test particle located at r′ at t = 0.¹² Our normalized number of contact changes is related to the integral of G(r,t) taken up to r = r_min; expressed as . The main quantity which is monitored in studies of liquid and glassy dynamics is the dynamic structure factor F(k,t) which is related to the van Hove correlation function G(r,t) via a spatial Fourier transformation. Since all the space integrations leave the qualitative behaviour as a function of time unaltered, we thus have the following relation Δz(t)/z_max ∼ [1 − F(k_max,t)]. This is almost exactly what we see in Fig. 5. This connection allows one to compare the qualitative decay in time of our contact-change parameter introduced here, with standard bulk dynamical parameters such as F(k_max,t). In future studies, our contact change parameter can be used in simulations and also in experimental systems such as colloids under confocal microscopy, to analyse the dynamical arrest transition in terms of mechanical contacts at the monomer level, thus bridging the mean-field dynamics of standard approaches like MCT with particle-based analysis of rigidity and nonaffine motions.^{20,21,43–45}

The aim of this numerical simulation study was to construct a minimal system that demonstrates the arrested dynamics of a collapsed polymer at the monomer scale, and explore the role of different type of bonding. Although one might argue that it would be useful to study the collapse process of multiple chains, the complexity of microscopic analysis would in that case increase to the point where it would be difficult to draw detailed, reliable conclusions. It would not be clear whether the frozen state is a result of interactions within a single chain or collective interactions and entanglements between several interacting chains. Another important issue to point out is that the thermostat temperature was kept constant (on average) during our Brownian dynamics simulation to avoid confusion and misinterpretation of the connectivity data. We changed the effective quenching depth by varying the LJ depth ε. It might be interesting to study the collapse dynamics at different temperatures. In this paper, our interest is mostly in the arrested dynamics dependent on the bond counting and not the glass transition temperature per se.

7 Conclusions

We have established that the linear polymer chain of length N will form a glass when each monomer, on average, gets approximately 2 + 8/N additional physical (non-covalent) attractive contacts, on top of the covalent bonds along the chain. In a dense system of attracting particles, this gives a total number of mechanical contacts z* = 4 + 8/N, different from the Maxwell's z* = 6 of purely central-force networks, or z* = 2.4 for purely covalently bonded system. The glass transition temperature can be easily derived from this criterion by using the thermal expansion relation for the packing density ϕ(T) ∝ exp[−α_TT]. We have verified the cooperative freezing of global dynamical (α-like) relaxation in the glassy state, and the retention of localized thermal motion akin to β-relaxation inside the confining cage for each particle, although the qualitative behavior differs from that of bulk simple liquids.

It needs to be emphasized that in our simulation we deliberately made an effort to model covalent bonds of the connected polymer chain, distinct from the central-force (LJ) interaction of any pair of particles in close proximity. This was achieved by the combination of FENE and bending potential for the chain-linking bonds, and we believe this is the reason for our main result of z* ≈ 4. There are many simulations (athermal and Brownian) of a colloid glass, of particles with only central forces, and those should have a different threshold.

This approach offers a different, much simpler and intuitive look at the glass properties and criteria of dynamical arrest in complex liquids. Furthermore, our view of the polymer folding into a mechanically-stable glassy state in terms of a quantitative rigidity criterion appears consistent with independent simulations studies⁴⁶ where the collapse was put in relation to the boson peak (excess of low-frequency soft modes) in the vibrational density of states. The excess of soft modes, in turn, correlates with z and with nonaffinity,⁴⁷ and our framework may lead to a unifying picture of the emergence of rigidity across the variety of glassy transitions in soft matter. It also has a direct relevance to folded protein rigidity, where certain regions (we assert: with more than 2 physical bonds between residues) are rigid, while others are flexible.

Acknowledgements

We are grateful for useful discussions and input of W. Goetze, H. Tanaka, D. Frenkel and W. Kob. This work has been funded by the Theory of Condensed Matter Critical Mass Grant from EPSRC (EP/J017639). Simulations were performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk).

References

1. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
2. E. A. DiMarzio and J. H. Gibbs, J. Polym. Sci., 1959, 40, 121–131.
3. E. A. DiMarzio and J. H. Gibbs, J. Polym. Sci., 1963, 1, 1417–1428.
4. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259–267.
5. W. Goetze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory, Oxford Univ. Press, Oxford, 2009.
6. J. L. Barrat, J. N. Roux and J. P. Hansen, Chem. Phys., 1990, 149, 197–208.
7. R. A. Riggleman, K. Yoshimoto, J. F. Douglas and J. J. de Pablo, Phys. Rev. Lett., 2006, 97, 045502.
8. K. Schmieder and K. Wolf, Kolloid Z., 1953, 134, 149–189.
9. W. van Megen, S. M. Underwood and P. N. Pusey, Phys. Rev. Lett., 1991, 67, 1586–1589.
10. W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 4134–4153.
11. L. Berthier and G. Biroli, Rev. Mod. Phys., 2011, 83, 587–645.
12. J. P. Hansen and I. R. MacDonald, Theory of Simple Liquids, Academic Press, London, 2005.
13. A. Lemaitre and C. Maloney, J. Stat. Phys., 2006, 123, 415–453.
14. A. Zaccone and E. M. Terentjev, Phys. Rev. Lett., 2013, 110, 178002.
15. J. P. Wittmer, H. Xu, P. Polinska, F. Weysser and J. Baschnagel, J. Chem. Phys., 2013, 138, 191101.
16. J. Wang, J. Li and S. Yip, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, 12627–12635.
17. J. Wang, J. Li, S. Yip, D. Wolf and S. Phillpot, Physica A, 1997, 240, 396–403.
18. H. B. Yu, R. Richert, R. Maaβ and K. Samwer, Phys. Rev. Lett., 2015, 115, 135701.
19. S. Saw and P. Harrowell, Phys. Rev. Lett., 2016, 116, 137801.
20. A. Zaccone and E. Scossa-Romano, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 184205.
21. A. Zaccone, J. R. Blundell and E. M. Terentjev, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 174119.
22. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 8th edn, 2005.
23. L. H. Sperling, Introduction to Physical Polymer Science, Wiley, New York, 2006.
24. M. Ikeda and M. Aniya, J. Non-Cryst. Solids, 2016, 431, 52–56.
25. S. R. Elliott, Physics of Amorphous Solids, Longman Scientific & Technical, London, 1990.
26. J. L. Barrat, J. Baschnagel and A. Lyulin, Soft Matter, 2010, 6, 3430–3446.
27. J. C. Phillips, J. Non-Cryst. Solids, 1979, 34, 153–181.
28. H. He and M. F. Thorpe, Phys. Rev. Lett., 1985, 54, 2107–2110.
29. G. H. Doehler, R. Dandoloff and H. Bilz, J. Non-Cryst. Solids, 1983, 42, 87–95.
30. A. Zaccone, Mod. Phys. Lett. B, 2013, 27, 1330002.
31. L. C. Hsiao, R. S. Newman, S. C. Glotzer and M. J. Solomon, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16029–16034.
32. T. G. Fox and P. J. Flory, J. Appl. Phys., 1950, 21, 581–591.
33. L.-P. Blancharde, J. Hesse and S. L. Malhotra, Can. J. Chem., 1974, 52, 3170.
34. J. Guan, Y. Wang, B. Mortimer, C. Holland, Z. Shao, D. Porter and F. Vollrath, Soft Matter, 2016, 12, 5926–5936.
35. S. J. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
36. A. Lappala and E. M. Terentjev, Macromolecules, 2013, 46, 7125–7131.
37. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086.
38. J. L. Keddie, R. A. L. Jones and R. A. Cory, Europhys. Lett., 1994, 27, 219.
39. P. G. De Gennes, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 2, 201–205.
40. M. Arutkin, E. Raphael, J. A. Forrest and T. Salez, Soft Matter, 2016 10.1039/c6sm00724d.
41. L.-M. Wang, F. He and R. Richert, Phys. Rev. Lett., 2004, 92, 095701.
42. B. A. Pazmino Betancourt, P. Z. Hanakata, F. W. Starr and J. F. Douglas, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 2966–2971.
43. M. L. Falk and J. S. Langer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 7192.
44. H. Yoshino and M. Mezard, Phys. Rev. Lett., 2010, 105, 015504.
45. H. Yoshino, J. Chem. Phys., 2012, 136, 214108.
46. C. Maffi, M. Baiesi, L. Casetti, F. Piazza and P. De Los Rios, Nat. Commun., 2012, 3, 1065.
47. R. Milkus and A. Zaccone, Phys. Rev. B, 2016, 93, 094204.

---