Network topology of the states probed by a glassy polymer during physical aging

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers. Link to publication


DOI: 10.1002/mats.201900036
Following Goldstein's idea, Stillinger and Weber [7][8][9][10] introduced the key feature of configurational mapping whereby arbitrary sets of molecular positions are referred, essentially uniquely, to local minima on the PEL for a many-body system. These authors introduced the term "quenching" for describing the steepest-descent path leading from an instantaneous configuration of the system to the corresponding local minimum, the so-called inherent structure (IS). In this way, the configuration space of the system can be split into basins, where a basin is defined such that a local minimization of the potential energy maps any point in the basin to the IS contained within it; the resulting tiling into nonoverlapping basins of attraction simplifies the thermodynamic description of the system. [11,12] In the PEL framework, one can focus on the ISs, at the bottom of the basins, wherein the system is trapped for most of the observation time. Owing to the separation of timescales, one is able to precisely describe many features of glass formers by only considering the properties of the ISs. [13] However, transitions are equally important, since relaxation (approach to equilibrium, i.e., "physical aging") occurs only via occasional jumps to neighboring basins. [6,[14][15][16][17] In dynamics, configurational changes can thus be separated into fundamental structural changes (transitions from an IS to the next IS [7,8,[18][19][20][21][22] ) and vibrational motions within the potential well (basin) of each IS. In the "landscape" picture of slow dynamics in supercooled liquids, [4,6,23] large-scale structural relaxation can be attributed to the complex pathways that connect basins.
Wales and co-workers, [23][24][25] have performed extensive studies on the morphology of the PEL of model glass-forming liquids. In particular they have observed that the landscape of Lennard-Jones liquids has a multifunnel structure. There exist many states, residing inside a funnel, which are separated by small barriers, while barriers separating funnels may by large. Heuer and co-workers, [12,26,27] have shown that the system may find itself trapped in a "metabasin" (cluster of connected basins) for long periods of time, making frequent small hops within the metabasin, and infrequent transitions from one metabasin to another. Thus, a set of basins communicating through fast transitions, relative to the residence time in the set, is grouped into the same metabasin. [5] The developed picture resonates well with the multifunnel observation. De Souza and Wales, [28] have correlated transitions between metabasins in PEL with specific changes in the configuration space, more precisely with the extent of cage-breaking (i.e., the molecular mechanism where the first neighbors of individual atoms are changing).

Introduction
One of the long-standing goals of theoretical and experimental research on (polymer) glasses has been to understand the relation between the time-dependent structure of the glass and the observed dynamics. Dynamical behavior of glassy systems can be considered in terms of the transient localization in basins of the potential energy, and transitions among them on longer time-scales. [1][2][3][4][5] Goldstein was the first to make the connection between the potential energy landscape (PEL) framework and the glass transition in a heuristic way, reasoning that the potential barriers separating basins dominate the dynamical properties at low temperatures (so-called β-processes). [6] We envision physical aging as the redistribution of the glassy specimen on its potential energy landscape. The ability of the sample to explore its PEL depends on two factors. The first being the morphology of the PEL, that is, the existence of valleys (local minima) and passages (saddle points). The second one dictates the time-dependent exploration process on the PEL, that is, the connections between minima and the way the glassy sample evolves as a function of aging time. In order to address the challenge of revealing pathways (giving rise to structural changes) connecting basins in the PEL and defining rigorously metabasins, we switch our attention to a rather remote field, that of graph theory. [29] The ISs visited by a glassy polymer specimen during physical aging are envisioned as the states (vertices) of a continuously expanding network (graph). Transitions from one basin to another naturally correspond to the directed edges of the graph. This level of abstraction allows us to focus on the topological characteristics of the exploration of the PEL, instead of its topography (which is fully defined by the chemistry and interactions of the model). Based on the concepts and methods developed in graph theory, a plethora of interesting features of that network can be revealed, e.g. the distribution of (incoming and outgoing) connections per state, percentage of bidirectional transitions and existence of specific patterns. [29] During the last decades, there has been a growing interest in modeling complex systems as networks (graphs), [30] deeply inspired by the discovery of Watts and Strogatz that many real-life networks share common collective characteristics and behave as "small worlds". [31] The designation "smallworld" reflects the observation that neighbors of a specific vertex are very likely to be neighbors to themselves (c.f. "six degrees of separation" of social networks). The feature refers to the connectivity of the network, irrespectively of the existence of any kind of physical distance between the vertices, and sets these networks aside from randomly-connected networks, where no emerging patterns can be identified. Another interesting characteristic, shared by a diverse range of networks, for example, the World Wide Web, [32] is the "scale-free" character of their topology. The distribution of the number of connections to each node, the node's degree, follows a power law. This topology results from the dynamics of the network growth in these systems. [33] In our study, the growth process is the process of physical aging, striving to bring the system to (thermodynamic) equilibrium. Formal definitions and criteria for quantifying the "small-world" and "scale-free" characteristics of the networks will be provided in the following section. A network may display neither of the two, only one or even both characteristics.
The statistical properties of networks have been the topic of considerable attention in the physics literature in recent years. Scala et al. [34] presented evidence that the (finite) conformation space of a 2D lattice homopolymer chain is a small-world network. In their approach, the vertices of the network had one-toone correspondence to all allowed conformations of the chain, and a link between two vertices indicated the possibility of switching from one conformation to another by a single Monte Carlo (MC) move of the chain. These authors found that the geometric properties of the network were similar to those of small-world networks, while locally the network appeared to have low dimensionality. However, the connectivity distribution was found to be Gaussian. [35] Doye analyzed the topology of the network formed by the minima and transition states of a 14-atom Lennard-Jones cluster. [36] The author found that this network exhibited both a small-world and scale-free character. The topology of the network, as in the case of Scala et al., [34] was static, that of the underlying PEL. Moreover it was observed that low-energy minima having large basins of attraction act as the highly connected hubs in the network. Massen and Doye [37] provided a detailed characterization of the inherent structure networks for a series of small Lennard-Jones clusters. Those networks showed a mixture of local order (as probed by clustering coefficients) and small average separation between nodes, that is, behavior characteristic of both small-world and scale-free networks. [38][39][40][41] Later, these authors examined how the potential energy landscape of a 13-atom cluster evolves with the range of the potential, thus tracking the transformation of the energy landscape from one with just a single minimum to a complex landscape with many minima and a scale-free pattern of connections. [42] They found that during this growth process, new edges in the network of connected minima preferentially attach to more highly connected minima, thus leading to the scalefree character. [43] However, their topology could not reflect the rules governing the dynamics of network growth, because they were static spatial networks.
Morgan et al. [44] created a database of minima and transition states for small clusters bound by the Morse potential. These authors treated the networks as being unweighted and undirected (thus allowing the possibility of bidirectional transitions between minima). They found that small-world character was present in those networks but they were not scale-free. They then compared the Morse networks with molecular and atomic glass-forming systems (ortho-terphenyl [45] and binary Lennard-Jones [28,46] ). These glassy systems did not show small-world properties, suggesting that such behavior is linked to specific aspects of the landscapes of the Morse clusters.
In this paper, we report the results of molecular dynamics (MD) simulations which provide information about the network of states visited by a sufficiently large polymer model system during its time evolution (i.e., physical aging). Regular quenches of the system are performed during the MD run, in order to map the MD trajectory onto a set of inherent structures. The interbasin dynamics can then be represented as a walk on a network whose nodes correspond to inherent structures and where edges link those minima which are directly connected by feasible transitions. The topology of the network emerges from the dynamics of the network growth. Graphtheoretical concepts are employed in order to characterize the topology of the network formed by the ISs visited by the system during physical aging. Using this approach can lead to substantial progress when combined with the tools of statistical mechanics and molecular simulations as it has already been demonstrated for other glass-forming systems, like the binary Lennard-Jones mixture. [28,46] Despite the fact that the PEL network is static, that is, it is simply determined by the number of atoms and the form of interatomic interactions, its time-dependent exploration starting from an initially localized distribution in phase space is a dynamical process, ultimately connected to the physical aging of (polymer) glasses. We do not aim for an exhaustive exploration of the energy landscape of our model system; rather, the network of visited states is augmented by the ISs discovered as time elapses. The sampling procedure focuses on a single route through the high-dimensional landscape, thus generating a limited subset of connections between minima. The graph is a representative sample of the space explored through aging. The system is not ergodic, thus we have to start from different initial configurations and take the arithmetic (not equilibrium-Boltzmann) average of the observables, as Theodorou and Suter proposed. [13] As far as physical aging and glassy relaxation are concerned, there are two features which are considered important: apparent infinite dimensionality and low connectivity of the network of accessible states. [15,23,[47][48][49][50][51][52] Small-world networks combine both of them and emerge as perfect candidates for the description of the network of states visited by a glassy specimen in the course of physical aging. On the contrary, it has been widely hypothesized that, in order to reproduce stretched-exponential decays (e.g., Kohlrausch-Williams-Watts [53,54] ), the space of accessible conformations must have a sparse random graph structure. [48,[55][56][57][58][59] In this work, we set out to study the characteristics of the network that consists of the ISs visited by a glassy atactic PS specimen in the course of aging. A better understanding of the exploration of the PEL can lead to useful physical conclusions, providing insight into complex mechanisms such as cage breaking and string motion in atomic systems. [28,[60][61][62] 2. Theory

Probing the Inherent Structures along the Trajectory
An atomistically detailed configuration of the system in the glassy state is obtained by quenching from the melt using MD. The well equilibrated melt configuration, of an arbitrary system can be obtained by any combination of MC or MD methods at elevated (above T g ) temperatures. Then, a MD run under constant pressure and progressively lower temperature is employed in order to drive the system into the glassy state, where it gets trapped within a specific basin of its PEL. The so prepared system is then used as an initial configuration for long MD trajectories. This procedure is completely general since classical MD simulations of all kinds of materials (e.g., polymers, liquid crystals, inorganic glasses) are nowadays widely available.
After achieving steady-state behavior (since the system is glassy there is no "true" equilibrium) at a specific temperature, we produce a discrete time series of configurations, R R t ( ), sampled along the trajectory of a standard MD simulation. Each of the configurations, R t ( ), is then mapped to its corresponding inherent structure, R t ( ) inh , by locally minimizing the potential energy in configuration space. In their classic work, [7] Stillinger and Weber used a steepest-descent method to find potential minima along a trajectory and demonstrated how to separate changes of the inherent structures from the vibrational motion. These authors used the term "quenching" for that procedure and should not be confused with the one used by Kirkpatrick et al. [63] in annealing processes of amorphous materials to their global minimum. The time-sequential quenching-procedure furnishes a tool to extract fundamental structure changes in liquid and disordered solid dynamics, by removing its vibrational components. It is meaningful only when the character of instantaneous structures is reflected in the corresponding inherent structures, which holds for polymers. [64] We perform the local minimizations by using a variablestep steepest-descent method and the quasi-Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, [65][66][67][68][69][70][71] and its limited-memory variant (L-BFGS). [72] The quasi-Newton methods are at least an order of magnitude faster in terms of CPU wallclock time, than the steepest-descent. However, the quenching trajectories leading from the instantaneous configuration to the IS produced by the quasi-Newton minimizers are not strictly the same as steepest-descent paths, because the Newton-like step is not parallel to the gradient vector. The steepest-descent path is a well-defined property of the potential energy surface, independent of coordinate system, as shown by Banarjee and Adams. [73] Sun, Ruedenberg, and Atchity [74,75] developed a steepest-descent method that employs second derivatives and tested it on some model surfaces. These authors pointed out that quasi-Newton paths to minima can deviate significantly from steepest-descent paths and can even converge to the wrong minimum on the 2D Muller-Brown surface. [76] Later, Asenjo et al. [77] have shown that deviations of the ISs probed by quasi-Newton algorithms from basins of attraction defined by steepest-descent pathways, can be controlled to some extent by adjustment of the maximum step size during minimization. In all cases, we employ a very conservative step of |ΔR j | = 10 −2 Å for the change of any degree of freedom, R j . Imposing a maximum step size to the quasi-Newton minimizers, can mitigate the discrepancy between the produced minimum-energy paths and the steepest-descent paths, at the cost of slightly higher computational effort. An appropriate value for the maximum step size can be chosen based on the intrinsic length scales of the system: for atomic systems, a good choice is about 1% of a carbon-carbon bond. In our calculations, quasi-Newton and steepest-descent minimizations converged to the same minima for a subset of configurations randomly sampled throughout all trajectories and conditions. The combination of the short step and the steepness of the boundaries of the PEL basins render both methods equivalent in our study. Since we are interested in probing the minimum-to-minimum connectivity of the expanding network of visited states, the exact location and characteristics of the minima are not of crucial importance. At every minimum, the Hessian (matrix of the second derivatives of the potential energy with respect to the atomic Cartesian coordinates) is calculated and diagonalized in order to verify that the stationary point is actually a minimum and not a transition state. Indeed, all quenching steps have ended up in local minima of the PEL, thus allowing us to define the network of the minima.
Since it is very time-consuming to determine many inherent structures of large systems, the time interval used in sequential quenching along a trajectory must be carefully chosen. If it is too small, many sequential configurations are quenched to the same inherent structure before jumping to the next inherent structure, and so it is computationally inefficient. On the contrary, if it is too large, many inherent structures can be missed. An optimal time interval depends on the temperature and size of the system. For a low-temperature trajectory of a polystyrene specimen we found it to be on the order of a few picoseconds. Thus, we employ the shortest interval of 1 ps in our calculations.
After having produced the quenched configurations R t i ( ) inh at equidistant time-points t i , we take a step further in order to probe the elementary hopping events. If , we cannot know whether , since many other minima could have been visited by the system between t i and t i+1 . In this case, we employ the interval bisection method of Doliwa and Heuer, [27,78] which involves further minimizations within the interval [t i , t i+1 ], till the resolution of a single MD step (10 −15 s). By storing all positions and velocities at time t i , we reconstruct the trajectory by integrating the equations of motion in the interval [t i , t i+1 ] during post-processing. We can then apply the interval bisection method for this short trajectory and probe all relevant IS hopping events that have occurred.

Graph-Theoretic Representation of the Network of Explored States
The time evolution of a glassy system can be considered as a sequence of infrequent events, those being hops from one basin of its PEL to another. Thus, the process of exploring the PEL can be cast in the form of a graph (network), where a set of vertices (nodes), v i ∈ V, are joined together in pairs by edges (links), e ij ∈ E (with i, j denoting specific vertices, v i and v j , respectively). The nodes represent the ISs of the basins visited by the specimen (as a function of time) and the edges the transitions occurring among them. We use lower-case letters with indices for enumerating vertices and edges (v and e, respectively) while using upper case characters for sets (V and E) and global properties related to the whole graph. The counts of vertices and edges are denoted as |V| and |E|, respectively. Research on reallife networks (ranging form social networks [31] to power grids [35] and the Worldwide Web [32] ) has focused on a number of distinctive statistical properties that most networks seem to share. [79] Networks can be classified based on their expansion with respect to time as "static", "quasi-static," or "growing." [80] In a static network both the number of nodes as well as the number of connections remain fixed and do not change with time. In the case of a quasi-static network, a fixed number of nodes are present at the initial stage. The connections are then introduced one after another between pairs of nodes using some specific probability distribution. Therefore, in a quasi-static network the number of nodes in the network is fixed, but the number of edges grows in time. Finally in a growing network the nodes are introduced in the network one after another. After introducing the node, a link is introduced connecting this node to one of the already introduced nodes. Therefore, in a growing network both the number of nodes as well as connections grow with time. In the case of the PEL exploration, the network of states can be classified as growing during physical aging.
In order to analyze the network of IS of the glassy polymer, a few concepts from graph theory are introduced. The degree of a vertex, k i , is the number of other vertices to which it is connected. In our case, each vertex of the PEL graph that evolves in the course of time, is associated with an incoming, k i in and outgoing, k i out , degree, referring to the number of edges (transitions) entering into and departing from the node, respectively. A particularly important feature of a graph is the shortest path between two vertices, d, defined as the smallest number of edges that must be followed for the system to navigate from one vertex to another. Based on the shortest-path construction, we can introduce the characteristic path length of the network, 〈d〉, defined as the number of edges in the shortest path between two vertices, averaged over all pairs of vertices (which can be connected by a path). Moreover, we can introduce the diameter of the graph, [81] d max , as the length of the longest shortest-path between any two graph vertices, that is, the length of the longest chain of transitions the system is forced to use to get from one vertex to another, excluding paths which backtrack, detour, or loop.

Identifying Small-World and Scale-Free Networks
In terms of their connectivity, graphs can be classified in four groups as shown in  of them are first neighbors. In this case, the diameter of the graph, that is, the greatest distance between any pair of vertices is unity, d max = 1. Reducing the number of edges by specifying only a limited count of neighbors per vertex, we end up with a regular graph, as in Figure 1b. Each vertex is connected to its k = 4 first neighbors and the corresponding diameter of the graph, d max , increases. The diameter of a k-regular graph is given by the greatest integer less than or equal to (3|V| − k − 3)/(k + 1) − 2. [82,83] In the example of Figure 1b, the diameter of the graph is equal to 3. By specifying the same number of edges and vertices, as those of the regular graph, we can create a random graph, whose edges (i, j) are randomly sampled from the population of all possible combinations of i and j. This gives the graph in Figure 1d, where there is a distribution of neighbors, k, per vertex, and the diameter of the graph increases even more. In ordinary applications, the connection topology of a graph is assumed to be either completely regular or completely random.
In between regular and random graphs, there exists a special group of graphs, the so-called "small-world" (SW) networks. A SW network is a type of graph in which most vertices are not neighbors of one another (like in regular graphs), but the neighbors of any given vertex are likely to be neighbors of each other ( Figure 1c). Thus, it is possible to connect any two vertices in the network through just a few edges. The local connectivity suggests the network to be of finite dimensionality. The average distance between vertices in a network is short, [84] usually scaling logarithmically with the total number of vertices. [31,85] SW networks can be highly clustered (like regular lattices), yet having small average path lengths, like random graphs. The small-world network of Figure 1c is created by applying the random rewiring algorithm of Watts and Strogatz [31] to the random network of Figure 1d. [86] SW networks exhibit clustering, or network transitivity, which refers to the increased propensity of pairs of vertices to be acquainted with one another if they have another acquaintance in common. The transitivity of a graph is closely related to the clustering coefficient of a graph, as both measure the relative frequency of triangles. Watts and Strogatz [31] defined a clustering coefficient that measures the degree of clustering of a graph. For our analysis we will refer to it as global clustering coefficient and define it by: [87,88] where "triangles" refers to triples of vertices each of which is connected to both of the others, and "connected triples" are triples in which at least one is connected to both the others. The factor of 3 in the numerator accounts for the fact that each triangle contributes to three connected triples of vertices, one for each of its three vertices. By introducing the factor of 3, the global clustering coefficient is normalized; its value, C, lies strictly in the range from zero to one. [88] It becomes unity on a fully connected graph (every vertex is connected to every other vertex) and has typical values ranging from 0.1 to 0.5 in many small-world real-life networks. [89] Moreover, the local clustering coefficient, c i , of a vertex v i can also be defined. [90] It is given by the proportion of edges, e jk , between the neighbors of v i , relative to the total number of possible edges between them: [31,88] (2) where k i out is the outgoing degree of vertex i and e jk is an edge formed by a pair of vertices v j and v k which are both neighbors of vertex i, that is, v j , v k ∈ N i with N i being the set of neighbors of i, N i = {v j : e ij ∈ E}. For a directed graph, e jk is distinct from e kj , and the degree of outgoing connections is used in the denominator.
The clustering coefficients introduced above are employed to tell apart the several classes of networks. For a random network, [31] it is expected that C ≃ 〈c i 〉 ≃ k/n where k is the average vertex degree of the network and n the number of vertices (c.f. Figure 1). In contrast, small-world networks have values of C, 〈c i 〉 on the same order of magnitude as those of regular lattices, because only a small percentage of connections are different from those in the lattice (c.f. Figure 1). From a quantitative point of view, a network exhibits "small-world" characteristics if the clustering coefficients, defined by equations (1) and (2), are significantly larger than zero, for example, on the order of 0.1. The actual values of the clustering coefficients can be interpreted in the sense that for example, C = 0.1, one out of every ten possible connections is present between the neighbors of a vertex.
Finally, irrespectively of being SW or not, many networks exhibit scale-free hierarchy, that is, the presence of hubs, or a few nodes that are highly connected to other nodes in the network. One finds that there are typically many vertices in the network with a low degree, and few of them with a high degree. The precise distribution follows a right-skewed power-law or exponential form, indicating that the network exhibits scalefree characteristics at several levels of organization. [30,33,35] For an undirected network, we can write the degree distribution as The power law P deg (k) remains unchanged (other than a multiplicative factor) when rescaling the independent variable k, as it satisfies P deg (ak) = a −γ P deg (k). Scale-free networks may emerge in the context of a growing network, in which new vertices connect preferentially to the more highly connected vertices in the network. If time is used to describe the evolution and expansion of the network, the fractal nature of these networks is sustained over all timescales. Empirical data obtained for several social, biological, and computational networks have confirmed the existence of such probability distributions decaying as power laws. [91][92][93] The exponent, γ, for these networks, varies between 2 and 3.
Clauset et al. [94,95] have suggested a statistical test based on the method of maximum likelihood for determining the best fit power-law to distributions on empirical data. The maximumlikelihood estimator for the exponent γ of equation (3) is [96] where n kmin is the number of nodes with degree at least k min and the sum runs over only those nodes. This formula gives an estimate of the best value γ, but gives no indication of the quality of the fit. Typically, a power law is only obeyed in the tail of the distribution, so a cutoff, k min , above which the power law applies needs to be determined. Clauset et al. also have a method for determining the best value of k min . [95] A fit is calculated for all possible values of k min and the value which gives the smallest Kolmogorov-Smirnov statistic is chosen. [97] Thus, based on the procedure of Clauset et al. we can have a reliable estimate of the parameters γ and k min , that is, the degree after which the distribution has the best power-law (or any other) fit.

Systems Studied
In this work, monodisperse glasses of atactic PS chains with 50% meso diads obeying Bernoullian statistics and a chain length of 300 (30 kg mol −1 ) repeat units (coarse-grained sites) (diads) were generated and simulated. It is a well-established fact that the glass-transition temperature, T g , of linear polymers increases with molecular weight, M w , and saturates at high M w . Hintermeyer et al., [98] studied the dielectric response of a series of polystyrenes covering a wide range of molecular weights including the monomeric limit. They found that polystyrene chains of M w > 30 kg mol −1 exhibit asymptotically low relaxation dynamics and above that threshold, the dynamics is not affected by the length of the chains. Bearing in mind the limitations of the atomistic MD simulations and for simplifying the presentation of the results, the main part of this work is based on systems composed of a single chain with a molecular weight of 30 kg mol −1 , although systems with five and ten chains have also been simulated. In general, larger systems are preferred, however, according to the current literature, smaller systems should be employed for the analysis of dynamics in configuration space, because otherwise many interesting effects may be averaged out. [26,99,100] Periodic boundary conditions are employed in all three directions of the simulation cell.

Polymer Melt Equilibration
For long-chain polymer materials, the computational prediction of structural, dynamical, and mechanical properties is very challenging because of the broad spectra of length and timescales governing structure and molecular motion. As in our previous works, [101,102] we have addressed the challenge of producing faithful PS melts by employing two interconnected levels of representation. First, a coarse-grained representation is employed, wherein each polystyrene repeat unit is mapped onto a single bead. [103] The smoother effective potential energy of the coarse-grained representation permits its equilibration at length scales larger than the repeat unit by using powerful connectivity-altering Monte Carlo algorithms. [104] Initial configurations were generated by a quasi-Metropolis gradual insertion procedure [101,105] at a temperature of 500 K, at which the coarsegrained forcefield has been developed. [103] Coarse-graining and reverse-mapping between the two levels of representation is accomplished in a manner that preserves tacticity and respects the detailed conformational distribution of chains. [101] Reverse mapping from a well-equilibrated coarse-grained configuration to a detailed atomistic configuration is accomplished in three stages. A detailed description of this procedure is presented in Ref. [101]. Generally, the first stage involves an iterative quasi-Metropolis introduction of the atomistic sites of the polymer, obeying the atomistic potential to be described below. In the second stage the generated configuration is optimized using local MC moves. In the final stage, the energy of the atomistically detailed configuration is minimized before initiating MD integration. Throughout the reconstruction procedure, CH 2 united atoms containing the achiral carbons of the chains are kept fixed at the positions of the superatoms of the coarse-grained configurations. Minimization is performed in rounds, by gradually increasing the atomic radii. [105] One starts with atoms of reduced size (sigma equal to one tenth of its actual value), adjusting the size in stages so that the atoms reach their full size in the end. The reverse mapping scheme just described was designed in order to prevent locking of the local configuration in torsional states which are inconsistent with the unperturbed conformational statistics adopted by PS in the melt, without departing at all from the well-equilibrated configurations provided by the coarse-grained simulations.

Atomistic MD Simulations
The detailed configurations obtained by reverse mapping can be described by a united-atom model without partial charges, based on the work of Lyulin and Michels. [106] This united-atom (UA) model will be referred to as "atomistic" in the following. Despite the fact that hydrogen atoms are not included in the description, it can fully capture the structure, dynamics, and mechanical properties of PS. [101,104,107] It takes into account the following contributions to the system potential energy: i) Lennard-Jones nonbonded interaction potential between all united atoms that are three or more bonds apart or belong to different images of the parent chain; ii) bond stretching potential for every covalent bond; iii) bending potential for all bond angles, including those in the phenyl rings; iv) torsional potential for all rotatable backbone bonds; v) torsional potential for the torsions of phenyl rings around their stems; vi) out-of-plane bending potential to preserve the coplanarity of the phenyl ring and the phenyl stem; vii) torsional potential about all bonds connecting aromatic carbons in the phenyl ring to preserve the planarity of the ring, and viii) improper torsional potential to preserve the chirality of all carbons bearing a phenyl substituent. [108] All Lennard-Jones potentials are cut at an inner cutoff distance of 2σ beyond which force smoothing to zero is applied using a quintic spline up to a distance of 2.5σ. [105,109] Ensuring continuous derivatives of the potential is crucial for successful convergence of the energy minimizations. We pay special care to ensure that discontinuous functions in any potential term are avoided in order to eliminate numerical instabilities during energy minimization. No tail corrections are used for the nonbonded interaction potential. Our experience has been that this united-atom model does a reasonable job in predicting structure, volumetric properties, elastic constants, and stress-s train behavior in the melt and glassy states.
All MD simulations have been conducted using the large-scale atomic/molecular massively parallel simulator (LAMMPS), [110] extended with the non-standard terms of the united-atom force field of Lyulin and Michels. [101,104] In all cases, a timestep of 1 fs was used. Initially, the reverse-mapped configurations were subjected to 0.5 µs of isothermal-isobaric (NpT ) MD under T = 500 K and p = 101.325 kPa, using the intrinsic barostat of LAMMPS which implements the equations of motion proposed by Shinoda et al. [111] The barostat couples all three diagonal components of the stress tensor when the pressure is computed, and dilates/contracts the dimensions together, maintaining the parallelpiped symmetry. The final configuration from the melt at 500 K was subjected to further NpT simulations with the set temperature T lowered by an effective cooling rate of 0.1 K ns −1 down to a final temperature of 200 K, through the glass transition, which was found to be 373 K for the system under consideration. [101] During the cooling procedure, configurations of the system were recorded at 300 and 200 K. These configurations were used for NpT equilibration under atomspheric conditions at both temperatures and subsequent NVT simulations.
The simulation results will be reported for constant volume (NVT ensemble) simulations at T = 200 K and T = 300 K. Conceptually, the analysis at constant volume is easier to perform than at constant pressure, since the appropriate free energy (Helmholtz energy) to be minimized does not contain any terms involving stresses applied to the system. Speedy [112] has shown that the number of minima on the potential energy landscape of a soft-sphere system is independent of the volume, under constant number of particles, and the shape of each basin varies with volume in a predictable way, assuming that the basins are harmonic. Tarjus et al. [113] examined the respective role of density and temperature in the viscous slowing down of glass-forming liquids and polymers. They concluded that density does not affect the fragility of glass-forming systems. Moreover, if the volume is kept constant in the quenching procedure, the structures are kept homogeneous, thus, the inherent structures must be a good representation of the corresponding instantaneous structures and there will be no discrepancy between the findings under constant volume and those under constant pressure.

Trajectory of the Inherent Structures
During the NVT runs, we monitor the time series of ISs, characterized by their respective potential energies, t V ( ) inh . A representative trajectory of the t V ( ) inh is presented in Figure 2. The trajectory of the ISs is characterized by a specific spectrum of potential energy values, opposed to the instantaneous potential energy of the MD trajectory, t V ( ) pot , which is presented in the inset to Figure 2. The latter quantity exhibits reasonable energy fluctuations, around 0.2k B T per interacting united atom (k B T ≃ 0.6 kcal mol −1 ), around its average value. During a thermostatted Molecular Dynamics run, the system can in principle overcome high energy barriers. However, accomplishing this requires time, proportional to the height of the barrier. Thus, within the finite simulation time, primarily low energy barriers connecting a basin to its neighboring ones are explored. The longer the simulation time, the higher is the chance of surpassing increasingly high energy barriers as well. Interestingly, there exist several time intervals during which the system is jumping back and forth between a limited subset of the available inherent structures, for example, in the interval [0.1, 0.2]µs. This is marked by the structured form made up of quantized bands of energy values in Figure 2.
At this temperature (T = 300 K), the details of the transitions between successive ISs are examined in more detail. We identify such transitions by quenching the MD configurations, as discussed in the previous section, and looking for structural signatures of a transition from one inherent structure to another. We are considering three observables in order to differentiate between inherent structures. First, we monitor the IS energy, V inh as a function of time, as shown in Figure 2. Second, we estimate the distance between two successive quenched configurations, where ∆ ≡ . [19] Based on the convergence limits and the numerical accuracy of our minimization method, we consider two ISs distinct, if the configurational distance, ∆R t N ( )/ inh is greater than or equal to 10 −7 nm. In the rare event where a transition occurs between two inherent structures with the same energy, configurational distance will still exhibit a sign of the transition, and for this reason we mainly use ∆R t N ( )/ inh to identify transitions. Third, we also use the Hessian to discern different minima. The Hessian matrix, that is, the matrix of the second derivatives of the potential energy with respect to the Cartesian atomic coordinates, is very sensitive to changes in shape of the basin.
Macromol. Theory Simul. 2019, 28,1900036  Comparison of the lowest eigenvalue suffices to judge whether two minima are the same, in case the criteria introduced above lead to contradictory results. When evidence of a transition was found within a time interval, this time interval was further divided by following the careful bisection method of Doliwa and Heuer, [27,78] as described in the Methodology section.
Having the trajectory of inherent structures at hand, we can then assign unique indices to every one of them. The corresponding results, as a function of time, for one trajectory are presented in Figure 3. It is evident that the number of different ISs grows almost linearly with elapsed physical time. However, the system still visits already explored ISs (once an IS is visited for the first time it is classified as "explored"), as this can be seen from the horizontal plateaus of the tags of the ISs (c.f. detailed inset to the Figure 3). These states seem to form distinct clusters, or in the graph-theoretic language "communities" of highly connected states. This level of organization can be connected to "meta-basins," which can be thought of as valleys on the PEL. [26] Clusters appear at multiple time scales as the system can escape from one of them, spend some time in a neighboring one and then return to the original, as happens between 0.1 and 0.2 µs. In that particular time framework the system keeps visiting several times two different subsets of states, with indices around 30 000. This traversal of the system gives rise to "islands" of states in Figure 3, one of them is magnified in the inset to the figure. Over longer timescales, the system finally escapes from a cluster of basins, revealing the irreversible nature of the dynamics in the glassy state.
Taking a step further, we keep track of how many times a specific IS is visited by the system. The occupation probabilities of ISs are presented in Figure 4, for three different time horizons. The trajectory of the IS is the same as that of Figure 2. The system evolves in time by starting from an initially narrow distribution of states, around the basin of the IS reached after quenching from the melt state. Even at t = 0.01 µs, the system has already discovered a few thousand states. However, these are not evenly populated. The probability distribution contains spikes at specific ISs, whose occupation probability exceeds by an order of magnitude the occupation probability of the rest. At longer timescales, t = 0.5 µs, the system finds its way to lowerenergy ISs, being the first signature of physical aging. The specimen does not visit all minima on the energy landscape; it just samples them along a single-dimensional stride through the high-dimensional space. The probabilities presented in Figure 2 correspond to the evolution of a single (notably long evolving) glassy specimen.

Graph Representation of the ISs Visited
Our analysis so far has been focused on characterizing the statistics (count and occupancy probabilities) of the ISs visited by the system during its time evolution. This has led us to the interesting observations that the system spends most of the time by occupying a limited subset of ISs, while it probes new states at a rate proportional to the elapsed physical time. However, in order to gain better understanding, it is important to characterize the topology, that is, the connectivity between ISs at different scales, of the underlying network of ISs. To this end, we envision the visited ISs as a directed graph, including the ISs as its vertices and the transitions between them as its edges. Our edges are directed, since the transition e ij  e i→j is inherently different from the reverse transition, e ji (different free energy barrier and therefore different rate constant).
A planar view of the corresponding graph is presented in Figure 5. The planar layout of the network enables visual inspection of the connectivity and provides a fast insight into the topology of the network of states. The usefulness of the representation is dependent on the aesthetics of the drawing. While there are no strict criteria for aesthetics, it is generally agreed that a helpful planar view should minimize edge crossing, and should evenly distribute the vertices on the plane. This problem has been studied extensively in the literature [114] and many approaches have been proposed. In this work, we choose to produce a layout of the graph by using the force-directed algorithm Macromol. Theory Simul. 2019, 28,1900036   of Hu. [115] A force-directed algorithm models the graph-drawing problem through a physical system of bodies with forces (harmonic springs) acting between them on the plane. The algorithm finds a good placement of the bodies in two dimensions by minimizing the energy of the system.
There exist some profoundly connected states, serving as hubs, and some states in funnel-like structures which bear only one incoming and one outgoing connection. The system is trapped in some stable configurations (those with high occupation probability of Figure 4), a subspace of its configuration space for long times, during which a small number of minima are visited over and over again. Obviously, these minima form superstuctures, which following Stillinger, [4] we denote as metabasins (MBs). Based on the graph-theoretic representation of the network, community structure arises naturally, as agglomerates of ISs, and well developed community-finding algorithms, like the one of Newman and Girvan, [116] can be employed in order to determine MBs as community structures in the graph. It is anticipated that the links between different "neighborhoods" of states will be more important for the relaxation of the system, than transitions within them. Moreover, the underlying configurational changes may be characterized by different levels of cooperativity. [117] The evolution (expansion) of the network of states is also important. Based on Figure 3, one could anticipate that new states are discovered at a constant rate, with their count increasing almost linearly in time. Based on Figure 6, the set of unique minima discovered at T = 300 K grows roughly as t 4/5 , the exact exponent of the best fit to the data-points in Figure 6 being 0.854. At short time scales, there are some fluctuations, mostly due to populating the same minima but eventually the exploration reaches the 4/5 exponent. The lower-than-linear dependence can hardly be discerned in Figure 3 due to the linear abscissa axis; the difference becomes more pronounced in the logarithmic axes of Figure 6. On the contrary, the count of different transitions (different in the sense that the system may discover several pathways connecting a specific pair of states that has already been explored) probed by the system is proportional to the elapsed time (the exact exponent of the best fit to the data-points in Figure 6 being slightly less than 1, i.e., 0.949). At lower temperatures the expansion of the network is slower, with the dominant scaling-law being t 3/5 at 200 K.
In the inset to Figure 6, the percentage of bidirectional edges is presented as a function of time. Initially, as more and more connections are discovered and the network is populated, the percentage increases up to 20% which appears to be an asymptotic value. We can draw the conclusion that most transitions can occur only in a specific direction for the time horizon examined in this study, the reverse transition being difficult (i.e., unlikely), implying the heterogeneous nature of the barriers transversed in the PEL.

Small-World Characteristics
In order to probe whether the network of ISs has "smallworld" characteristics, we calculate the global, C, and average local, 〈c i 〉, clustering coefficients according to Equations (1) and (2), respectively. Both coefficients are calculated at regular time intervals during the expansion of the network. They  both converge to the same asymptotic value (around 0.17, see Figure 7); the average local clustering coefficient exhibits more fluctuations, due to its local character. The current accepted definition of a small-world network is that it has clustering similar to a regular lattice, that is, higher than an equivalent random network, c.f. discussion in the Methodology section. In order to prove that the network of ISs has "small-world" characteristics, we compare it with a random network of the same size (number of vertices) and connectivity (average number of edges per vertex). To this end, for every point in the physical time, we create an ensemble of a few random networks with the same size and connectivity as the IS network at that point in time, following the procedure presented in the Methodology section, and calculate the average global and local clustering coefficients, which are shown in the inset to Figure 7. The measured clustering coefficients of the PEL network (main part of Figure 7) are at least two orders of magnitude higher than the ones obtained for equivalent random networks (inset to the Figure 7). Thus, we infer the existence of a small-world structure (i.e., neighborhoods of states and only few connections between different neighborhoods) and rule out a purely random structure for the network of explored states. [118] To verify whether or not the networks indeed exhibit smallworld behavior, the average shortest-path length is also considered. The average shortest-path lengths for two trajectories at T = 200 and T = 300 k are presented in Figure 8. Their dependence on the number of states discovered is not monotonous. Initially, the system discovers new states, out from the quenched state in the very beginning, increasing linearly the average shortest-path length of the network. While the system remains trapped within a specific metabasin, connections are established between already existing states; thus the average shortest-path length decreases, for example, between 75 and 100 states for 300 K or between 20 and 75 for 200 K. Once the system finds its way out of a metabasin, the average shortest path increases again. The overall behavior is on average sublinear, but faster than logarithmic. The study of the average shortest-path length is not conclusive, in the sense that we cannot fully match the behavior of the exploration network to either random (linear), or small-world (logarithmic) networks (behavior). Based on Figure 8, the average shortest-path length seems to scale with the square root of the number of states of the network.

Scale-Free Characteristics
Within a network, some nodes may be more highly connected than others are. To quantify this effect, we determine the probabilities P in (k) and P out (k) that an IS has k incoming and outgoing links (via transitions), respectively (k being the degree of the node as defined above). The results presented in Figure 9 show that both P in (k) and P out (k) follow a power law over several orders of magnitude, remarkably different from the Poisson distribution predicted by the classical theory of random graphs, [119] and the bounded distribution found in models of random networks. [31] The points in Figure 9 tend to fall well along the power-law line. However, there exist fluctuations, especially the horizontal points for large degrees, k > 10 2 (k is the abscissa of the Figure 9), as there are few points (ISs with so many transitions) to average out the noise. Although most nodes have a very small degree, there are a few nodes with a degree above 500. The presence of hubs that are orders of magnitude larger in degree than most nodes is a characteristic of scale-free networks. Thus, the time-dependent PEL traversalconnectivity is dominated by highly-connected ISs. The same holds if incoming connections are considered instead of the Macromol. Theory Simul. 2019, 28,1900036   outgoing ones. The overall system obeys scaling laws characteristic for highly interactive self-organized systems. [120] We now examine the exponent of the power-law tail of the distribution of the number of connections (transitions) for each node. The exponent of the power law, γ ≃ 2, is similar to other scale-free networks. [91] By comparing the main part of Figure 9 to its inset, we can readily observe that the power-law exponent of the distribution is independent of the physical time elapsed (i.e., aging time). This decay has also been observed for the Internet, [121] metabolic reaction networks, [122] the telephone call graph, [123] and the World Wide Web. [124] Remarkably, the exponent, γ, lies in the range 2.1−2.4 for all of these cases. Barabàsi, Albert, and Jeong, [33,125] call these networks "selfsimilar", by analogy with fractals, phase transitions, and other phenomena where power laws arise naturally and there exist self-similarity over multiple scales (or no single representative characteristic scale can be defined).
Following Clauset et al., [94] an additional goodness-of-fit test should be employed to check whether the proposed power-law is a plausible fit for the measured data (c.f. Methodology section). Datasets are generated from the power law, fitted with the same procedure and their fit is compared to the fit of the measured data. The p comp value is the fraction of generated datasets that have a worse fit than the measured data. Clauset et al. [94] suggest that the fit should be rejected if the p comp value is less than 0.05. Higher values, that is, p comp ≥ 0.05, do not necessarily indicate that the power law is a good fit, merely that there is not sufficient evidence to reject it. By following the procedure of Clauset et al., we find that the exponent of the degree distribution is γ = 2.37 ± 0.08 and the power-law distribution is a good fit for k > k min , with k min = 52 (Equation (4)). The probability, p comp , is 0.76 indicating that there is no reason to reject the power-law fit of the degree distributions.

Conclusions
In this work, we have developed a graph-theoretic analysis of the network of inherent structures explored by a glassy polymer specimen during physical aging. Summarizing, our results reveal the following points: i) the network of explored states exhibits clear community structure with neighborhoods of closely connected states which can be identified as metabasins, ii) there exist pathways involving strings of states leading from one neighborhood to another, iii) the distribution of connections per state of the network has a scale-free character, and iv) the network possesses small-world network characteristics. In other words, the network of explored states is extremely heterogeneous with a few well-connected states, acting as hubs, but with the majority of nodes only connected to a relatively small number of ISs. Models of network growth [125] and studies on the time evolution of real networks [126] strongly suggest that the heterogeneity at the heart of the scale-free topology develops as a result of new nodes preferentially linking to those nodes which have many connections. The intuition is that the existence of short paths, due to highly connected states, contributing to the clustering of the network, is responsible for high-speed communication-channels between ISs of the PEL, thereby facilitating any dynamical process that requires global coordination. Highly connected states bear a lot of connections to neighboring states within the same cluster; however, only one of them links the hub of a neighborhood to the hub of another neighborhood. This single transition may be qualitatively different from intra-neighborhood transitions, for example, relying on collective and concerted rearrangements of the particles instead of uncorrelated single-particle displacements. This may help elucidating the missing links between chemical structure, cooperativity, and physical aging.
The observations drawn within this work can pave the way for a methodology able to deal with the, otherwise inaccessible, timescales of the glassy state. The representation of physical aging as a hopping process between discrete ISs on the PEL can yield computationally tractable simulations of the extremely large relaxation times of (organic and inorganic) glasses based on atomic-level structure and interactions. Instead of letting a MD simulation spend most of the time being trapped within a limited volume of the phase space, an exploration process can be developed which will create a network of connected states by overcoming local barriers. [64] The first step toward this direction was accomplished by the authors, [127] and our present work provides design rules on how networks of states for glassy materials can be prepared and analyzed. The characteristics of the network of states studied in this work should be reproduced by the sought procedure in order to faithfully reproduce the true dynamics. Figure 9. Distribution of the number of states having k in incoming and k out outgoing connections of the minima discovered by a glassy PS specimen aging at T = 300 K, during a time horizon of 0.5 µs. In the inset to the figure, the same distribution is presented for a time horizon of 10 ns.