Microscopic carriers of plasticity in glassy polystyrene

The relevant parts of the microstructure participating in the plastic deformation of glassy polystyrene are identiﬁed. This is achieved by performing quasi-static deformation simulations of an atomistic representation of polystyrene, within a thermodynamic framework, under experimentally realistic boundary conditions (controlled compressive strain, atmospheric pressure in the lateral directions, and room temperature). In particular, it is shown that discontinuities in the free energy in the course of deformation are representative of microscopic plastic events, which are strongly correlated with the heterogeneous non-aﬃne deformation of the phenyl-rings. However, despite the induced anisotropy during deformation, the phenyl-rings do not show any preferred orientation during deformation. The rearrangements of the microstructure during the discontinuities are analyzed further in terms of a local best-ﬁt deformation gradient tensor (a modiﬁed version of the procedure introduced by Falk and Langer [1] ) as well as the stretch-and rotation-components of the local deformation of clusters of atoms that show signiﬁcant non-aﬃne deformation. This analysis shows that the local deformation is highly heterogeneous, and locally the strains and rotations can be orders of magnitude larger than the imposed deformation, giving rise to plastic deformation of the material.


Introduction
Polymer materials are increasingly used in high performance, load bearing applications where plastic deformation is generally not acceptable.It is therefore essential to understand the microscopic mechanisms underlying plastic deformation.Understanding the mechanism of plasticity can be of great value in the development of new materials with an improved resistance to plastic deformation.In contrast to crystalline polymers, where plastic deformation is known to occur due to slippage of crystallographic planes facilitated by imperfections in the crystal lattice (dislocations) [2] , relatively little is known about the microscopic mechanisms of plastic deformation in glassy (amorphous) polymers, which do not exhibit such a crystal structure due to the lack of long-range order.Unravelling the mechanisms of plasticity is a topic which has received considerable attention in the past, and has led to different views of the possible mechanisms behind plastic deformation in glassy materials.The first ideas on plasticity date back to 1930, when Eyring proposed that the rearrangement of matter is a process during which a potential energy barrier must be surpassed [3,4] .As shown by Krausz and Eyring [5] , external stress can reduce the height of an energy barrier and thereby facilitate the system to move between different potential energy minima.Despite the fact that vanishing energy barriers are in-deed related to plastic deformation [6][7][8][9] , the potential energy landscape (PEL) framework does not provide direct information about which parts of the microstructure are involved in the plastic deformation.A significant amount of computational work has been done in trying to understand the microscopic origins of plasticity in glasses, from a structural point of view.Research has shown that plastic deformation in glasses involves the local, non-affine rearrangements of atoms [1,[10][11][12][13][14] .Falk and Langer [1] observed locally irreversible shear transformation zones in a Lennard-Jones glass.By fitting a linear deformation gradient tensor to a group of atoms, they calculated the amount of deviation from a linear strain field for each atom and used it as a measure for the irreversible deformation.Papakonstantopoulos [10] observed a saw-tooth pattern in the stress-strain curve of a Lennard-Jones glass, and found that the drops in stress occur at very small strains (10 −4 ) and that these stress drops are associated with local areas of irreversibly deforming atoms, leading to local relaxation of the structure.The locality of these irreversible rearrangements suggests that glassy materials posses a heterogeneous structure composed of domains with different stiffness, where the low-stiffness domains are more prone to irreversible deformations.Despite the fact that these low-stiffness domains are microscopically small but do affect the macroscopic response of the material, continuum-based approaches have proven to be useful (provided they are used on a sufficiently large length-scale) by using effective material properties [15,16] .However, a better understanding of the behavior of glassy materials on a microscopic scale would be helpful for improving these continuum-based models.Low-stiffness domains are indeed observed by multiple researchers [10,14,[17][18][19] , and some researchers even propose that these so-called "soft spots" fulfill the same role as dislocations do in crystalline materials: They locally weaken the structure of the material and thereby facilitate local rearrangements.For instance, Papakonstantopoulos [10] found that irreversible rearrangements occur at regions with low (and sometimes even negative) stiffness.Manning and Liu [19] analyzed vibrational modes in a quasi-static, athermal sheared glass and found that regions with lowfrequency vibrational modes serve as weak spots where irreversible rearrangements occur.A significant amount of computational work described in literature focuses on relatively simple non-bonded Lennard-Jones glasses containing only a single type of interaction potential [8][9][10]17,20,21] , and therefore cannot describe the complex interactions and interplay between atoms that exist in an actual (bonded) polymer material. The studis that do focus on actual bonded polymer systems [22][23][24][25][26][27][28] often use a significant level of coarse-graining (i.e.bead-spring models) [22,[26][27][28][29] .Furthermore, most of the computational work employs molecular dynamics simulations which suffer from unrealistically high strain rates, or molecular mechanics simulations under non-physical (athermal or constant volume) boundary conditions [6,19,23,25,26,28,[30][31][32][33] .In this paper, we employ a free energy minimization framework [34] that allows us to investigate the behavior of the microstructure of a glassy polystyrene sample during plastic deformation, under physically relevant boundary conditions: applied quasi-static uniaxial deformation on one side of the simulation box, while allowing the lateral sides to expand freely under ambient pressure (-101.325kPa) and room temperature (300 K). An impotant difference between the thermodynamic framework applied here, and the computational methods (molecular dynamics and molecular mechanics) commonly used in modeling the behavior of polymer materials is that our procedure does not suffer from non-physical boundary conditions such as fixed Poisson's ratio (molecular mechanics) or unrealistically high strain rates (molecular dynamics).Furthermore, we employ a realistic, multi-potential multi-term forcefield for modelling the complex interactions between all united atoms.In order to quantify the behavior of the microstructure, we focus on the deformation of individual atoms as well as on the deformation of entire clusters of interacting atoms.In order to quantify the deformation of clusters of atoms as a whole, we use a modified version of the procedure proposed by Falk and Langer [1] . The origial procedure is commonly used in order to quantify, by means of the residual from the least-squares fit of the deformation gradient tensor, the non-affine deformation of a group of particles with respect to a center particle.In this paper, we are primarily interested in the collective deformation of clusters of interacting atoms.Therefore, we modified the initial method such that it yields a best-fit deformation gradient tensor for a cluster of interacting atoms as a whole.To be more specific, we quantify the deformation of all atoms in a cluster with respect to the center-of-mass of that cluster.This modification results in a best-fit deformation gradient tensor that is less sensitive to the deformation of individual atoms than the initial procedure, and captures the behavior of the entire cluster.By decomposing the obtained deformation gradient tensor into a rotational part and an elongation part, we are able to determine the angle of rotation and principal stretch ratio of individual clusters of atoms. This paper is organized s follows. In Sction 2, we briefly outline the most important aspects of the computational model and the thermodynamic framework under which the simulations are performed.In Section 3, we discuss the results of the deformation simulations and the results of the microstructural analysis.We conclude this paper with Section 4, in which we will summarize the most important conclusions and give an outlook on future directions.

Computational aspects 2.1 Computational model
For the deformation simulations, we employ an orthorombic simulation box with periodic boundary conditions on all sides of the simulation box.The simulation box contains a total of four isotatic polystyrene (iPS) chains in the glassy state.Although iPS has the ability to (partly) crystallize, the cooling rate of 0.1 K ns −1 used to quench from the melt is sufficient to suppress any form of crystallization.Each individual polymer chain has a molecular weight of 8.3 kg mol −1 and is composed of 80 repeat units per chain (see Figure 1).The procedures for obtaining the equilibrated melts and subsequent equilibration of the solid-state model are described in detail in Ref. [34,35].A combination of Monte Carlo methods ensures proper equilibration of the configurations at all length scales in the melt state; torsion-angle distributions and chain dimensions are in good agreement with 2-dimensional NMR measurements and neutron scattering experiments, respectively.The quenching process progressively locks the melt configurations without altering their characteristics at long length-scales [35] .The final configurations are represented by a united-atom model (where hydrogen atoms are not explicitly modeled but mapped onto their corresponding carbon atoms).The interactions between the atoms are described by interaction potentials based on the work of Lyulin and Michels [32] , The individual terms in Equation ( 1) represent the following contributions to the total potential energy, respectively [34] : (I) bond-stretching potential for covalent bonds; (II) bending potential for all bond angles; (III) torsion potentials for all rotatable backbone bonds, torsions of phenyl rings around their stems and about all bonds connecting aromatic carbons in the phenyl ring to preserve the planarity of the ring; (IV) non-bonded interaction potential for united atoms which are more than three or more bonds apart or belong to different images of the same parent chain; (V) potential to preserve planarity of the phenyl ring and the stem (out of plane bending); (VI) improper potential to preserve chirality.All Lennard-Jones potentials are cut-off at a distance of 2σ LJ beyond which a quintic spline, that brings the interaction to zero in a continuous way, is applied up to a distance of 2.5σ LJ , where σ LJ denotes the distance parameter where the pairwise Lennard-Jones potentials vanish between atoms.Mixed potentials for the interaction between dissimilar atoms are defined according to the Lorentz-Berthelot rules [36] .We emphasize here that, given the importance of the energetic interactions between all atoms for the local plastic rearrangements, a united atom model is used which has shown to be representative for the behavior of polystyrene [34] .Besides not modeling the hydrogen atoms explicitly (i.e., using united atoms), no other form of coarse graining is employed.

Deformation framework
In this section, we present the main ingredients used in the thermodynamic framework for the deformation simulations.For a complete and detailed discussion the reader is referred to Ref. [34].The deformation procedure is as follows (see also the flowchart in Figure 2).Once we have obtained a solid state specimen in the reference state R, with volume V R , the system is deformed, in incremental steps, up to large-strains (ε > 0.3) by a uniaxial compressive deformation (no shear components).This deformation is imposed by scaling the edge length L of the simulation box in the direction of deformation, to a new length l , by a factor λ = l /L , such that the prescribed strain increment ∆ε = ln(l /L ) = −0.001.During the deformation of the box, the lateral box dimensions are fixed, that is, ε ⊥ = 0 (in both lateral directions).Scaling the box dimensions in the direction of deformation results in an affine deformation of the atoms inside the simulation box.Since the configuration after the affine deformation does not necessarily constitute a potential energy minimum, the potential energy is minimized after each incremental strain step.This potential energy minimization step (with the lateral dimensions of the box still fixed), results in a non-affine deformation of the atoms.After the potential energy is minimized, the following objective function (Legendre transform of the Helmholtz energy) is evaluated to see whether it has reached a minimum [35] : If the objective function has not reached a minimum, the lateral dimensions of the box are adjusted, resulting in a change in ε ⊥ , under atmospheric pressure, σ ⊥ = −101.325kPa, while keeping the longitudinal dimensions fixed.Equation (2) applies to the deformation protocol outlined above: Imposing a longitudinal compressive deformation (without shear components), ε , while imposing a constant (and equal) stress, σ ⊥ , on the cross-section of the simulation box.For a thorough discussion on different deformation protocols and the corresponding objective functions, the reader is referred to Ref. [34] .This procedure, i.e. optimizing the lateral dimensions of the box, under an imposed constant longitudinal strain step, is repeated until the objective function is minimal, after which a new strain increment is imposed.All these small incremental strain steps lead eventually to a large-strain deformation of the specimen.The objective function presented in Equation ( 2) is the Legendre transform of the Helmholtz energy (see Ref. [35] ).The first term in Equation (2) represents the potential energy of the energy-minimized configuration after an applied strain step.We denote this potential energy as the potential energy at the inherent structure, adopting the terminology of Stillinger and Weber [37][38][39] .The second term represents the vibrational Helmholtz energy due to thermal vibrations in the system.The vibrational contribution is calculated by applying the quasi-harmonic approximation (QHA) to the free energy, which has been successfully used in the past for several polymer glasses [34,40,41] .As far as polytyrene is concerned, the QHA is valid even at elevated temperatures (the interested reader is referred to the thorough discussion in our previous works [34,41] ).In principle, the first two terms yield the required Helmholtz energy: It denotes the Helmholtz energy of a specimen trapped in a basin on the potential energy landscape [41] .The final two terms in Equation ( 2) appear only because of applying a Legendre transform in order to obtain a Flowchart of all steps during the deformation simulations.Calculation of A vib requires the normal modes of the system, which are obtained from the eigenvalues of the Hessian matrix containing the second-order derivatives of the potential energy with respect to the atomic coordinates [34] .
Helmholtz energy expression with the applied strain, ε , and atmospheric pressure, σ ⊥ , as the appropriate control variables (see Ref. [34] ).
3 Results and analysis

Compression simulations
Figure 3 shows the results from the compression simulations.The evolution of the Legendre transform of the Helmholtz energy as function of the imposed strain (see Equation ( 2), hereafter simply referred to as free energy) is shown in Figure 3(a).Figure 3(b) shows the relation between stress and strain, which is derived from the free energy (explanation in the text below).The components of the free energy, shown in Figure 3(c)-(f) represent the following contributions, respectively: (c) Vibrational contribution to the Helmholtz energy (second term in Equation ( 2)), (d) inherent potential energy contribution (first term in Equation ( 2)), (e) internal energy, and (f) entropic contribution.
The entire free energy curve is characterized by smooth, continuous parts, where the free energy increases with increasing deformation, interrupted by sudden drops, where the magnitude of the drops in free energy varies significantly despite the fact that the imposed strain increments (∆ε = -0.001)are constant throughout the entire deformation.Similar observations regarding this saw-tooth behavior were made by other researchers [10,30,33,34,42,43] .These drops in free energy are the result of mechanical instabilities, arising due to vanishing energy barriers as a result of the external deformation.It is known that energy barriers, in the form of first-order saddle-points, separate inherent structures and decrease linearly with stress [6,8,9,44,45] .When the imposed stress is sufficiently high, the energy barrier becomes so small that any further distortion of the system results in a drift of the system into a neighboring inherent structure, which has a lower energy minimum than the initial inherent structure.This phenomenon is known as a fold-catastrophe [6,7,9] .One can derive a stress-strain relation from the free energy as function of strain by considering the fact that the stress is given as the first-order derivative of the free energy with respect to strain [34] .This implies that a discontinuity in the free energy should correspond to a discontinuity in the stress-strain curve.Figure 3(b) shows the stress as function of strain up to ε ≈ −0.016.The stresses are derived from the free energy by taking the difference quotient of the free energy versus strain, between two consecutive strain increments.The inset shows the corresponding part of the free energy as function of strain.The transition between different inherent structures is irreversible, as shown in Figure 4.For a few (randomly selected) energy drops we have reversed the deformation (from compression to tension, onset of reversed deformation shown by the red stars; the dots mark the onset of reversed deformation for reversible parts of the curve.)just after the point where the energy dropped.Upon reversing the deformation, the system follows a different energy path than it did upon compression (blue curve), indicating that the system has deformed irreversibly (i.e. has deformed plastically).The strains at which these plastic deformations occur is far below the strains at which plasticity is observed on a macroscopic scale, which for polystyrene lies somewhere around ε ≈ 0.05 [46] as indicated by the shaded area in Figure 4.A very similar observation can be made when considering rate-dependent plasticity in the form of stress-activated hopping and the Eyring viscosity: the plastic strain rate (creep rate) is non-zero even at stresses significantly below the nominal yield stress.Furthermore, we point out that the microscopic specimen represents a small, local region, of a macroscopic specimen with a relative small number of atoms (N = 2564), and that even the smallest local rearrangements are therefore immediately visible in the form of free energy drops.On the other hand, in a macroscopic specimen, which contains an enormous number of atoms (N 2564), local rearrangements are not immediately visible in the free energy, or stress, and the observed behavior is a much smoother curve of the free energy, or stress, as function of strain.However, that does not mean that local plasticity is not present in macroscopic systems at small strains, but the discontinuities are not reflected in the free energy, or stress, since these local rearrangements generally involve a relatively small number of atoms compared to the bulk.In order to show that not all incremental deformation steps lead to transitions between different inher-     4: Free energy (from hereon the free energy is shown as a dimensionless intensive quantity, normalized by the number of atoms in the box, N , the Boltzmann constant, k B , and the temperature, T ) as function of strain during compression (blue) and the reversed trajectories during tensile deformation (red).When the deformation is reversed after a free energy drop (onset of reverse deformation marked by red star), the reversed trajectory follows a different path than it did upon compression.The last two parts show reversible behavior (onset of reversed deformation marked by red dot).That is, the reversed trajectory follows the exact same energy path as in compression.
ent structures, we reverse the deformation along two continuous parts of the free energy curve (onset of reversed deformation shown by red dots), starting at ε ≈ -0.26 and ε ≈ -0.3.Upon reversing the deformation along the continuous parts, the system follows the exact same energy path as it did during compression, meaning the deformation is reversible (i.e.elastic).During these reversible deformations, the shape of the energy landscape is also distorted, but in contrast to the discontinuous parts, the energy barrier stays intact and the system resides in its current inherent state.

Atomic displacements
Our aim is to investigate the microstructural behavior during the irreversible deformations as identified in the previous section.In the remainder of this paper, we focus on the behavior during the irreversible parts in the free energy curve.As shown in the previous section, the irreversible parts are characterized by a decrease in free energy.Therefore, we define an irreversible event as two consecutive strain increments, where the free energy decreases.We start with decomposing the total displacement of each atom into two components: An affine part (the part that follows the imposed deformation) and a non-affine part (the part that deviates from the imposed deformation, recall that the non-affine deformation occurs as a result from the minimization of the potential energy, as discussed in Section 2) as follows: where ∆ r i is the total displacement of an individual atom, measured between two consecutive strain increments, r i is the position vector of each individual atom, F denotes the deformation tensor that results in the affine deformation of the microscopic specimen (not to be confused with F k , the deformation gradient tensor used for the deformation of clusters, see further below), and is of the form F = λ e i e i + λ ⊥ e j e j + λ ⊥ e k e k , where λ denotes the imposed deformation in the loading direction, and λ ⊥ denotes the change in the lateral box dimensions during the minimization of the potential energy, thereby accounting for the lateral boundary condition.
The system considered here is described by a multi-potential forcefield containing complex interactions such as bond-stretching and bond-angle bending interactions, imposing constraints on the system in terms of deformation.These complex interactions make it unlikely that complete affine deformation will ever occur.However, the non-affine deformation during the reversible parts is negligible (individual atomic displacements upon deformation being on the order of 0.01σ LJ , where σLJ denotes the average distance parameter where the pairwise Lennard-Jones potentials vanish between atom types; cross interaction terms are obtained according to the Lorentz-Berthelot mixing rules [36] ) compared to the irreversible parts.
Figure 5(a) shows the magnitude of the maximum non-affine displacement found at each irreversible event, as function of the magnitude of the irreversible event (that is, the magnitude of the drop in free energy between two consecutive strain increments).It is immediately clear that there is a (non-linear) correlation between the maximum non-affine deformation and the magnitude of the irreversible events.Since the data shows a clear monotonic increase, we use Spearman's correlation coefficient [47] , ρ, to quantify the correlation.The Spearman coefficient can take on any value between -1 (negative correlation) and 1 (positive correlation), with 0 indicating no correlation.Applying Spearman's relation yields ρ = 0.67.This indicates that the maximum value of the non-affine displacement is indeed correlated with the magnitude of the free energy drops.Instead of only measuring the contribution of the maximum non-affine displacement, which measures the influence of a single atom, we are also interested in the number of participating atoms during an irreversible event.In order to quantify the number of participating atoms during an irreversible event, we use the participation factor [10] P non−affine = 1 N where ∆ r non−affine denotes the non-affine displacement of atom i (at a specific strain increment) and N is the number of atoms in the simulation box.The participation factor is an indicative measure for the amount of atoms that show a non-affine displacement component greater than the bulk (majority of the atoms).The participation factor approaches zero if a relatively small number of atoms is participating in the non-affine deformation while the number of atoms in the bulk (N ) is large, and approaches unity when the distribution of non-affine displacements is uniform (all atoms participate equally).Intermediate values (between zero and unity) provide an indication of the fraction of atoms that participates in the non-affine displacements.The participation factor versus the magnitude of the free energy drops is plotted in Figure 5(b).Applying Spearman's coefficient to the participation factor yields a correlation coefficient of 0.02, indicating practically no correlation at all.From the preceding results it is clear that the magnitude of the free energy drops (irreversible events) are strongly correlated with the maximum non-affine displacement.In order to get a more detailed picture of which parts of the microstructure are involved in the non-affine deformation, we store (at each strain increment) the atoms with a non-affine displacement component within a pre-defined range, and group them based on their chemical composition.The results are shown in Table 1.In this table, we show two situations that differ in their normalization, as explained in the following.In the first case, denoted as "Normalized over all", we count all atoms (denoted as N tot ) that have a non-affine displacement larger than specified by the "non-affine displacement" threshold.From the total number of atoms (N tot ), we sort all the atoms based on the atom type, which we denote by N type .The fractions listed as "Normalized over all" are defined as N type /N tot .These fractions give information about how many of the non-affinely deforming atoms are of a certain atom type.The second case, denoted as "Normalized by type", differs from the first case in that we use a different normalization.Rather than normalizing by the total number of non-affinely deforming atoms N tot , we normalize N type by the total number of atoms of that particular type present in the system, denoted as N total−type .The fractions listed as "Normalized by type" are found as N type /N total−type .This second normalization provides information about what fraction of a certain atom type participates in the non-affine deformation.As an example: In the non-affine displacement range from 1.4σ LJ -2.1σ LJ , the "Normalized over all" situation gives a fraction of 1 for the amount of participating CH(phenyl ring), which indicates that all atoms with a non-affine deformation component within the specified range are phenyl-ring atoms.The "Normalized by type" situation gives a fraction of 0.0033, indicating that only a very small fraction of the total number of CH(phenyl ring) is participating.We use these two types of normalization to account for the fact that not every atom type occurs equally often.The data in Table 1 shows that almost 70% of the atoms with a non-affine displacement component between 0.2σ LJ -2.1σ LJ are located in the phenyl-rings.However, due to the large number of phenyl-rings present in the system (compared to other atom types), only a small fraction of the total number of phenyl-rings participates in the non-affine displacements.Nonetheless, from these results it follows that the bulky phenyl-rings suffer the most in terms of non-affine deformation, during the irreversible events.Normalized over all 0 0 0 1 0 Normalized by type 0 0 0 0.0033 0

Rotation of phenyl-rings
A more detailed analysis of the deformation of the phenyl-rings is performed by measuring the angle of rotation of each individual ring during deformation.In order to determine the angle of rotation, an orthonormal basis (three unit vectors perpendicular to each other) is constructed and attached to each ring, as shown in Figure 6.The angle of rotation of a basis (and therefore the phenyl-rings) between two consecutive strain increments can be extracted by constructing an orthogonal (rotation) matrix R (with R −1 = R T ) that maps the basis at the current strain increment onto the basis at the previous in- Vector c is oriented perpendicular to b and lies in the same plane as the phenyl-ring.Vector a is determined from the cross product between b and c and yields therefore a vector perpendicular to both a and b, and therefore perpendicular to the phenyl-ring.
Similar as for the non-affine displacements, one can define a participation factor for the angle of rotations.The participation factor now reads: where N rings stands for the number of phenyl-rings.The maximum angle of rotation as well as the average angle of rotation, measured over all phenyl-rings in the system, versus the magnitude of the discontinuities are shown in Figure 7(a).The participation factor for the rotation versus the magnitude of the discontinuities are shown in Figure 7(b).A similar trend is observed as for the non-affine deformations (Figure 5(a) and (b)).Applying Spearman's correlation to the results shown in Figure 7 yields a coefficient of ρ = 0.64 for the maximum rotation versus the magnitude of the discontinuities (Figure 7(a)), and a coefficient of ρ = 0.67 for the average angle of rotation of the phenyl-rings.Spearman's correlation coefficient for the participation factor versus the magnitude of the free energy discontinuities (Figure 7(b)) is ρ = 0.02.These results indicate that the phenyl-ring with the maximum rotation has a significant effect on the magnitude of the discontinuties, whereas the amount of phenyl-rings participating is almost completely uncorrelated with the magnitude of the discontinuties.
Applying the deformation on one side of the box induces anisotropy.One would expect to see this anisotropy in the orientation of the microstructure during deformation.In order to check if there is any orientation present, we extract the axis of rotation of each phenyl-ring (axis of rotation is the eigenvector belonging to the unit eigenvalue of the rotation tensor we used to rotate the bases onto each other), and use Herman's orientation factor to determine the amount of orientation between the axis of rotation and the vector along the compression direction.Herman's orientation factor is given as [49] f where φ denotes the angle between two axes; the angular brackets denote the average over all rotation axes (of the phenyl-rings).An orientation value f H = 0 means that the two axes are orientated randomly with respect to each other, an orientation value of f H = 1 indicates that they are aligned and a value of f H = −0.5 indicates that they are perpendicular.We measure the orientation of the rotation axis (of the phenyl-rings) with respect to the loading direction.The results for all three loading directions are shown in Figure 8.The results show that there is no clear orientation of the phenyl-rings during deformation.Although the orientation is not completely random (which would correspond to f H = 0), the average orientation and the average standard deviation remain constant during deformation.If there would be any preferred orientation of the rings during deformation, one would expect that more and more rings would align in the preferred direction, and the mean value of the Herman's factor would change during deformation.Since that is not the case, we conclude that the anisotropy induced during deformation is not reflected in the orientation of the phenyl-rings.

Deformation of clusters
Up till now, we have mainly focused on the motion of individual atoms and phenyl-rings.However, one of the key ingredients in our model, compared to other works, is the complex multi-term interaction potential that describes the interactions between all the individual atoms.The clusters are defined using two criteria.The first criterion is a non-affine displacement threshold that selects only those atoms with a non-affine displacement component greater or equal than the pre-defined cut-off criterion (all atoms with a non-affine displacement component smaller than the threshold value are still present in the system but are not considered in the clustering analysis).The second one is a distance-based threshold that determines which atoms form a particular cluster.The actual procedure is as follows.At a particular strain increment we define the non-affine displacement threshold and consider only those atoms with a non-affine displacement component larger or equal than the threshold value.We pick a random atom from the ones that meet the threshold value, and assign that atom to cluster 1.All atoms within a distance smaller or equal than the distance-based threshold to the atom in cluster 1, are also assigned to cluster 1.This process is repeated until, at a certain moment, no more atoms are within the pre-defined distance of any of the atoms in cluster 1.From the remaining (not clustered) atoms, we again pick a random atom and assign that to cluster 2. We repeat this process until all atoms that meet the non-affine displacement threshold are assigned to a specific cluster.When all atoms are assigned to a certain cluster, we quantify the deformation of each individual cluster (with respect to the previous strain increment) by a best-fit deformation gradient tensor that maps the atoms of a certain cluster in the reference state, onto the positions of the atoms in that cluster in the deformed state.Consider the deformation of an arbitrary cluster, k, composed of N k atoms, in an incre-  Herman's orientation factor applied to the angle between the rotation axis of the phenyl-ring (red) and the loading direction (black).Different plots show different loading directions.The blue line (left axis) shows the mean Herman's orientation factor, averaged over all phenyl-rings at every strain increment.The orange line (right axis) shows the standard deviation.Despite the fact that unidirectional deformation is imposed, the axis of rotation does not show any particular alignment during deformation.mental strain step ∆ε = −0.001from ε to ε + ∆ε .Departing from the procedure developed by Falk and Langer [1] , we use for the deformation gradient tensor of cluster k the expression It is noted that the expressions in the parentheses are (summations over) dyadic products of vectors, and therefore F k is indeed a tensor.With Equation ( 8), we try to map the positions of the atoms in cluster k at strain level ε , onto the positions of the atoms in cluster k one strain increment later, ε + ∆ε .Comparing this expression to the expression given in Ref. [1], the following modification is made: instead of using an atom inside the cluster as a reference, we use the center-of-mass (COM) of that cluster as a reference point, and map all relative position vectors δ r ε k,j (we define the relative position vector as δ r k,j ≡ r k,j − r k,COM ) onto the relative position vectors one strain increment later, δ r . The reason for introducing the center-of-mass as the reference point is that, when using the center-ofmass as reference, the result of the fit is not as sensitive to individual outliers (individual atoms with a strong heterogeneous deformation compared to the rest of the atoms in the cluster) as would be the case when using an individual atom as a reference point.Note that, opposed to the non-bonded Lennard-Jones glasses often encountered in literature, we are considering a complex system, composed of bonded and non-bonded contributions, where strong heterogeneous motion of (individual) atoms is likely to occur.Note that Equation ( 8) is derived from a least-squares fitting procedure [1] .Therefore, the atoms can, in general, not be mapped to their exact positions.In order to quantify the quality of the fit, we use the normalized residual Since the best-fit deformation tensor is a homogeneous deformation tensor, the normalized residual serves as a measure for how homogeneous the deformation of a cluster actually is.A residual of D2 k = 0 means the deformation gradient tensor describes the exact deformation of a cluster between two strain increments (i.e. the deformation is completely homogeneous).In practice it turns out that the maximum value for the normalized residual is D2 k = 1, which corresponds to a worst-case fit (i.e. the deformation is heterogeneous).The residuals for each cluster at each individual strain increment are shown in Figure 9.We consider two types of deformation: Short-range deformation and long-range deformation, which are defined by means of the distance-based threshold for clustering.In case of the short-range deformation we specify the distance-based threshold value to be 0.2σ LJ , which is slightly larger than the covalent bond length, so we only select atoms that are a bond length away from each other.For the long-range deformation the distance-based threshold is 1.7σ LJ , which is equal to the non-bonded interaction distance.The difference between the two types of deformation is that in case of the long-range deformation, cluster might (but not necessarily) consist of atoms from multiple chains (recall that the specimen consists of four different chains), whereas in the case of short-range deformation the cluster primarily contains atoms from a single chain.The residuals shown in Figure 9 can take values between zero (homogeneous, perfect fit) and unity (heterogeneous, worst case).We define the deformation to be homogeneous for D2 ≤ 0.1 (see red line in Figure 9).That is, an order of magnitude smaller than the worst case.
From the results we see that the majority of the clusters show a strong heterogeneous deformation ( D2 > 0.1), for both the long-range and short-range deformation.However, in the case of short-range defor-  mation, the clusters with a residual D2 ≤ 0.1 contain more atoms (up to 20) than in the case of longrange deformation (up to 10).The more homogeneous motion for the short-range deformation can be attributed to the fact that the motion in a single chain is more concerted because of the covalent bonds.From these results we see that the deformation gradient fitting procedure works reasonable for clusters that contain up to 10 atoms, when one does not discriminate between short-range and long-range deformation.The heterogeneous motion of the microstructure explains also why we determine the deformation inside a cluster with respect to the center-of-mass, instead of a single atom.In the original procedure, proposed by Falk and Langer [1] , an individual atom is used as reference.However, the fit of the deformation gradient tensor would be strongly affected by the motion of a single, heterogeneously moving atom.This influence is strongly reduced by using the center-of-mass as a reference, although the heterogeneous motion still affects the fit to a significant extend.For the clusters that deform sufficiently homogeneously (i.e.clusters with a D2 ≤ 0.1) the fit of the deformation gradient tensor is considered to be accurate enough to provide meaningful results in terms of deformation.Using the polar decomposition on the obtained deformation gradient tensors that are sufficiently accurate ( D2 ≤ 0.1) allows to decompose the total deformation into its components related to rotation and stretch.The results of this analysis are presented in Figure 10 and Figure 11, which can be found the Appendix.From this analysis it follows that there is no clear correlation between a certain type of deformation (rotation or stretch) and the magnitude of the discontinuities in the free energy.What the results in the Appendix do show is that some of the clusters experience strains that exceed unity.In view of the rather small value of the imposed strain increment, ∆ε = −0.001,this shows how large some of the local rearrangements are.

Conclusions
In this study we have shown, by means of deformation simulations on a glassy polystyrene system under physically relevant boundary conditions and described by a realistic multi-potential forcefield, that discontinuities in the free energy of a glassy polymer are representative of irreversible plastic events which are strongly correlated with non-affine deforming atoms located primarily in the bulky phenyl-rings.The phenyl-rings show a strong rotational motion during plastic deformation.Although we impose a unidirectional deformation, the phenyl-rings show no particular orientation during deformation.
We have proposed a modification of the procedure of Falk and Langer [1] in order to assess the deformation of clusters of atoms, which has shown that the majority of clusters shows strong heterogeneous motion during the irreversible events.Our results also show that when applying this procedure to bonded polymer systems, in contrast to the relatively simple Lennard-Jones glasses studied in literature, the deformation cannot be described accurately.For the few clusters that did show a reasonable homogeneous deformation, it was shown that the deformation (in terms of strains and rotational angles) is orders of magnitude larger than what is imposed on the system during deformation, leading to significant local rearrangements and in turn to plastic deformation.The relevance of this study lies primarily in the fact that we have shown which parts of the microstructure are contributing most to the plastic deformation, namely the phenyl-rings.This kind of information (i.e. which parts of the microstructure are responsible for plastic deformation) is very valuable in order to understand the molecular mechanisms of plastic deformation and in designing materials that can better withstand plastic deformation.Secondly, we have shown that the modification of the procedure by Falk and Langer [1] (presented in the previous section) is not able to quantitatively describe the deformation of clusters if the deformation is strongly heterogeneous.Finally, we have illustrated the potential of a rigorous thermodynamic framework which allows us to mimic deformation experiments under physically relevant boundary conditions, and which enables the investigation of the microstructural behavior of glassy polymers in great detail.Furthermore, the ability to perform numerical simulations on realistic polymer systems under physically relevant boundary conditions, as presented in this paper, provides great opportunities in bridging the gap from the microscopic level to the macroscopic level.

Appendix
In this Appendix we show the results of the polar decomposition, applied to the clusters that deform sufficiently homogeneously, i.e., clusters that show a residual D2 ≤ 0.1 as described in Section 3.4, Figure 9.Using the polar decomposition on the obtained deformation gradient tensors that are sufficiently accurate allows to decompose the total deformation into its components related to rotation, R k (with R −1 k = R T k ), and stretch, U k, (with U T k = U k ), from which the total deformation tensor can be re-constructed as From the rotation tensor we can determine the angle of rotation with Equation ( 5) by replacing R with R k .From the eigenvalues of U k one can extract the principal stretch ratio's λ k , where i = 1, 2, 3, which can then be used to determine the logarithmic strain in all three principal directions for each individual cluster.We emphasize again that the deformation for each individual cluster, at a certain strain level ε , is determined with respect to the previous strain increment ε + ∆ε .The results for the long-range clusters are shown in Figure 10, where Figures 10(b) and 10(c) only show the results for the clusters in Figure 10(a) that have a residual D2 ≤ 0.1.Figure 10(d) shows the maximum strain versus the rotation angle of each individual cluster.Note that some of the clusters (data points) show large rotations and strains.The reason for this is that clusters can consist of non-bonded atoms, making the deformation of such clusters less concerted, and a large displacement of a single atom can therefore already significantly influence the deformation of that cluster as a whole.The different colors represent different threshold values for the non-affine displacement of the atoms.From these results it follows that there is no relation between any type of deformation (stretch or rotation) and the magnitude of the discontinuities, as indicated by the Spearman coefficients, ρ (although the correlation coefficients are mathematically correct, one has to interpret the correlation coefficient carefully since some of the datasets contain relatively few datapoints).That holds also for the relation between a specific type of   deformation and non-affine displacement threshold.It is noted that one can observe large rotational angles and large strains for some of the clusters (as compared to the strain increment ε = −0.001); a significant amount of the clusters in the case of long-range deformation yields strain values of ε (i) k ≈ 1, and in some cases even strains that exceed unity.This shows how large the local rearrangements are, even when very small strains are applied.Similar conclusions hold for the short-range deformation, as shown in Figure 11.Although the deformation is somewhat more homogeneous (more clusters with a residual D2 ≤ 0.1), there is no clear correlation between any type of deformation and the magnitude of the discontinuities, or between deformation and non-affine displacement threshold.That is, in both the longrange and short-range deformations, the local deformation is strongly heterogeneous and extremely magnified compared to the applied deformation.

Figure 1 :
Figure 1: Model used for the deformation simulations.Periodic boundary conditions are employed on all sides of the box.Different colors indicate different chains.The close-up shows the united-atom representation, where the hydrogen atoms are mapped onto their corresponding carbon atoms.

Figure 2 :
Figure2: Flowchart of all steps during the deformation simulations.Calculation of A vib requires the normal modes of the system, which are obtained from the eigenvalues of the Hessian matrix containing the second-order derivatives of the potential energy with respect to the atomic coordinates[34] .

Figure 3 :
Figure 3: Free energy and its contributions, and the stress, all as function of the imposed strain.Note that the horizontal axis shows, from left to right, increasing compression, i.e., increasingly negative strain.(a) Free energy.(b) Stress as a function of strain, for small values of strain.The inset shows the free energy as function of the applied strain and shows the two discontinuities that correspond to the (vertical) dotted lines in the stress-strain figure.A more detailed explanation of this figure is given in the main text.(c) Vibrational contribution to the free energy.(d) Inherent contribution to the free energy.(e) Internal energetic contribution to the free energy.(f) Entropic contribution to the free energy.

Figure
Figure4: Free energy (from hereon the free energy is shown as a dimensionless intensive quantity, normalized by the number of atoms in the box, N , the Boltzmann constant, k B , and the temperature, T ) as function of strain during compression (blue) and the reversed trajectories during tensile deformation (red).When the deformation is reversed after a free energy drop (onset of reverse deformation marked by red star), the reversed trajectory follows a different path than it did upon compression.The last two parts show reversible behavior (onset of reversed deformation marked by red dot).That is, the reversed trajectory follows the exact same energy path as in compression.

Figure 5 :
Figure 5: (a) Maximum non-affine displacement, per strain increment, plotted versus magnitude of the free energy drop, ρ = 0.67.(b) Participation factor plotted versus magnitude of the free energy drop , ρ = 0.02.Figures contain results for three different deformation experiments: compression in x-,y-and z-direction.

Figure 6 :
Figure 6: Illustration of how the orthonormal basis is constructed onto the phenyl-rings in order to determine the angle of rotation.Vector b is oriented along the stem of the phenyl-ring and therefore lies in the same plane as the phenyl-ring.Vector c is oriented perpendicular to b and lies in the same plane as the phenyl-ring.Vector a is determined from the cross product between b and c and yields therefore a vector perpendicular to both a and b, and therefore perpendicular to the phenyl-ring.

Figure 7 :
Figure 7: (a) The filled circles represent the maximum angle of rotation of a phenyl-ring (left axis), per strain increment, as function of the magnitude of the discontinuities, with ρ = 0.64.The open triangles represent the average rotation of all phenyl-rings (right axis), as function of the magnitude of the discontinuities, with ρ = 0.67.(b) Participation factor applied to the rotation of the phenyl-rings versus the magnitude of the discontinuities, with ρ = 0.02.Figures contain results for three different deformation experiments: compression in x-,y-and z-direction.

Figure 8 :
Figure8: Herman's orientation factor applied to the angle between the rotation axis of the phenyl-ring (red) and the loading direction (black).Different plots show different loading directions.The blue line (left axis) shows the mean Herman's orientation factor, averaged over all phenyl-rings at every strain increment.The orange line (right axis) shows the standard deviation.Despite the fact that unidirectional deformation is imposed, the axis of rotation does not show any particular alignment during deformation.

Figure 9 :
Figure 9: Normalized residual for the fitted deformation of clusters.Each data point corresponds to an individual cluster.Each fit is determined at a certain strain increment (horizontal axis) with respect to the previous strain increment, as discussed in the text.The colorbar indicates the number of atoms inside each cluster.(a): short-range deformation.(b): long-range deformation.Note that the number of clusters found in the case of short-range deformation (a) is significantly higher than in the case of long-range deformation (b).This is due to the lower value for the distance-based threshold used in clustering the atoms.

Figure 10 :
Figure 10: Long-range deformation.Figure legends show the corresponding Spearman correlation coefficient, where (-) indicates that there are not enough datapoints for a correlation coefficient.Colors indicate threshold values for non-affine displacement: 1.4σ LJ ( ); 1.2σ LJ ( ); 1σ LJ ( ); 0.7σ LJ ( ); 0.5σ LJ ( ); 0.2σ LJ ( ).(a) Normalized residual as function of the magnitude of the free energy drop.The red line indicates the value below which the residual is considered to be small enough in order for the fit to be acceptable.(b) Angle of rotation for those clusters that have a normalized residual below the red line in Figure (a).(c) Maximum principal logarithmic strain for those clusters that have a normalized residual below the red line in (a).(d) Angle of rotation versus logarithmic strain.

Figure 11 :
Figure 11: Short-range deformation.Figure legends show the corresponding Spearman correlation coefficient.Colors indicate threshold values for non-affine displacement: 1σ LJ ( ); 0.7σ LJ ( ); 0.5σ LJ ( ); 0.2σ LJ ( ).(a) Normalized residual as function of the magnitude of the free energy drop.The red line indicates the value below which the residual is considered to be small enough in order for the fit to be acceptable.(b) Angle of rotation for those clusters that have a normalized residual below the red line in (a).(c) Maximum principal logarithmic strain for those clusters that have a normalized residual below the red line in (a).(d) Angle of rotation versus logarithmic strain.

Table 1 :
Fractions of atoms participating in a certain range of non-affine displacements, specified by atom type.Two different normalization procedures are used, which are explained in the main text.