Mutual Information in Molecular and Macromolecular Systems

The relaxation properties of viscous liquids close to their glass transition (GT) have been widely characterised by the statistical tool of time correlation functions. However, the strong influence of ubiquitous non-linearities calls for new, alternative tools of analysis. In this respect, information theory-based observables and, more specifically, mutual information (MI) are gaining increasing interest. Here, we report on novel, deeper insight provided by MI-based analysis of molecular dynamics simulations of molecular and macromolecular glass-formers on two distinct aspects of transport and relaxation close to GT, namely dynamical heterogeneity (DH) and secondary Johari–Goldstein (JG) relaxation processes. In a model molecular liquid with significant DH, MI reveals two populations of particles organised in clusters having either filamentous or compact globular structures that exhibit different mobility and relaxation properties. In a model polymer melt, MI provides clearer evidence of JG secondary relaxation and sharper insight into its DH. It is found that both DH and MI between the orientation and the displacement of the bonds reach (local) maxima at the time scales of the primary and JG secondary relaxation. This suggests that, in (macro)molecular systems, the mechanistic explanation of both DH and relaxation must involve rotation/translation coupling.


Introduction
The nature of the solidification process observed at the glass transition (GT) temperature T g by cooling supercooled viscous liquids is a topic of intense research. However, a clear characterisation of all the phenomena related to the glass transition as well as a deep microscopic understanding of the glassy state remains an open challenge [1][2][3].
Both frequency response and time relaxation of viscous liquids close to GT exhibit clear nonlinear features, thus motivating intense experimental [4][5][6], simulation [7,8], and theoretical [3] research. Within the framework of the familiar linear response, the fluctuationdissipation theorem evidences that two-time correlation functions incorporate all of the information to evaluate the susceptibility of a system in thermal equilibrium [9]. In principle, as outlined by Kubo in his seminal paper on adiabatic linear response theory, a formal treatment of the adiabatic nonlinear response in terms of a series expansion involving two-, three-, four-, and higher-order time correlation functions may be developed [10]. However, these expansions are exceedingly difficult to translate into a useful, experimentally verifiable form and are nowadays believed to be not a viable approach for most transport processes and irreversible relaxation processes [11]. This poses the question of alternative approaches to revealing nonlinear aspects of relaxation and response.
Starting with Galton's first expositions in 1888 on how to accurately characterise the dependence (correlation) between two or more random variables, a plethora of measures have been developed [12]. Perhaps, the most popular is the Pearson correlation coefficient: where σ X is the standard deviation of X. A disappointing issue of the Pearson coefficient, however, is that it measures only linear dependence.
To overcome this limitation, mutual information (MI) attracted quite a large interest. MI is defined as follows [13]: where p(x, y) is the joint probability distribution of the random variables X and Y with distributions p(x) and p(y), respectively. Two random variables with no MI are independent. MI can be thought of as the reduction in uncertainty (ignorance [9]) about one random variable given knowledge of another [12,14]. Interestingly, the data processing inequality states that manipulation of the data does not increase MI [15]. Even if developed in information theory and applied to maximise the amount of information transferred in communication, the notion of MI extended influence to many other fields [15]. A basic motivation is the property that MI assumes no sort of underlying distribution or dependence between random variables [16]. In particular, it is known that the nonlinear dependence leads to rather subtle effects [12]. In this respect, Figure 1 offers further illustration.
The sensitivity of MI to nonlinear dependence is certainly appealing in studies on GT [17][18][19] as well as in other fields of physics, including topological transition in the XY model [20], phase transition in a 2D disordered Ising model [21] and evaluation of the configurational entropy of liquid metals [22]. On the other hand, unlike the Pearson correlation coefficient, the MI evaluation requires knowledge of the distributions of random variables. If these are not known, their estimation is a delicate (and sometimes uncertain) procedure leading to higher computational cost [13,15]. Building on previous studies carried out by Molecular Dynamics (MD) simulations, the present paper reports novel insight provided by MI-based approaches on two distinct aspects of transport and relaxation close to GT, namely dynamical heterogeneity [23,24] and its influence on secondary relaxations [25,26]. Details about the model systems and numerical methods are found elsewhere [23][24][25][26].

Dynamical Heterogeneity and Mutual Information between Particle Displacements
The transition from a liquid to a glass is accompanied by the growth of transient domains that exhibit different mobility. This phenomenon is usually dubbed "dynamical heterogeneity" (DH) and has extensively been studied [3,[27][28][29]. The size of the DH domains is relatively small, involving about 10 molecule diameters [28], corresponding to a few nanometers [30] and strictly related to the possible presence of characteristic length scales in glass-forming systems [31]. Even if growing static length scales have been reported by experiments [32] and simulations [33], there is still debate about whether they control the glass transition [34]. On the other hand, it is not clear to what extent dynamic correlations are a consequence or the primary origin of the slow dynamics occurring close to the GT [35].
In a pioneering work in 2015, Dunleavy et al. [18] proposed a MI approach to investigate DH in a polydisperse mixture of hard spheres. They considered the MI between pairs of particles displacements, obtained a dynamic length scale growing on approaching GT, and revealed two distinct fractions of particles, so-called early and late fractions. To test its robustness, we extended the MI approach of Reference [18] to the context of molecular liquids and investigated by MD simulations [23,24]. For readers' convenience, these previous results are briefly outlined in Section 2.1. Then, Section 2.2 presents novel results proving that the so-called early and late populations are arranged in localised clusters.

Mutual Information Reveals DH in a Molecular Liquid
MD simulations of a model molecular liquid made of fully flexible trimers have been carried out with technical details given elsewhere [23,24]. Relaxation and heterogeneous transport are characterised by the intermediate scattering function (ISF) and the non-Gaussian parameter (NGP) [36], respectively. ISF is defined as where r j (t) is the position of the jth particle at time t. The brackets denote suitable ensemble averages, and N is the total number of particles. In an isotropic liquid, ISF depends only on the modulus of the wavevector q = ||q|| and provides a convenient relaxation function to study the rearrangements of the spatial structure of the fluid over the length scale ∼ 2π/q. We define the structural relaxation time τ α by the relation F s (q max , τ α ) = e −1 where q max is the maximum of the static structure factor. In the present system, q max 2π/σ, with σ being roughly the particle diameter. NGP is defined as follows: α 2 (t) vanishes if the particle displacement δ r(t) is Gaussian. NGP is a well-known DH metrics [1,3,29]. By suitably adjusting both the density ρ and the temperature T of the liquid, six states were prepared to be grouped in three pairs, labelled A, B, and C, with coinciding ISFs (and then τ α ) and NGP; see Figure 2 and References [23,24] for details. Figure 2(left) shows that ISF has a characteristic plateau region at intermediate times.
Increasing the relaxation time, the plateau widens, signalling that particles are trapped for longer times by the surrounding particles (cage effect). The subsequent escape from the cage leads to the structural relaxation and the ISF drop at t τ α . The escape from the cage is strongly non-Gaussian, i.e., heterogeneous, as evidenced by the NGP growth around τ α , revealing DH, i.e., the presence of particles with different mobility. The data were taken from Reference [23].
We now show how MI can be employed to characterise DH. Additional details are given in References [23,24]. We consider the MI between the displacements of two generic particles i and j I ij (t) = I(δ r i (t), δ r j (t)), where δ r m (t) is the displacement of the mth particle within a time interval t, and MI from Equation (2) is estimated by the Kraskov-Stögbauer-Grassberger (KSG) estimator [37]. The particles are said to be correlated at time t if I ij (t) > I 0 with I 0 = 0.2 [18]. We define n i (t) as the number of particles that are MI-correlated to the generic ith at time t. The MI-correlated particles exhibit a tendency to group around the generic one; see Figure 3(left). We denote with p(n, t) the distribution of particles being MI-correlated to other n particles at time t. The distribution is evaluated over four iso-configurational ensembles (ICEs) of the liquid. Each ICE is represented by the particle displacements generated starting from the same initial overall particle configuration, assigning about 10 3 different initial velocities to each particle according to the pertinent Maxwell distribution and evaluating the subsequent 10 3 distinct time evolution of the liquid. It has been shown elsewhere that the average over the four ICEs are, within the errors, equivalent to the customary ensemble average, i.e., it adequately samples the phase space of the system [23].
Representative plots of p(n, t) at different times are given in Figure 3(right). At vibrational times, t ∼ 1, the distribution peaks around the average number of surrounding particles hit by the tagged particle during the oscillatory motion within the cage where it is trapped. At longer times, still shorter than the structural relaxation time, t τ α , the tagged particle establishes further correlations through collision and p(n, t) spreads at higher n values. At t τ α , when (on average) the central particle escapes from the cage, p(n, t) narrows. Later, at t τ α , the distribution widens again. Finally, at t τ α in the diffusive regime, most correlation is lost and p(n, t) approaches a time-independent narrow shape, peaking at the number of particles with permanent MI-correlation with the tagged particle, i.e., all of the particles but one belonging to the same molecule, n = 2.
Insight is provided by n(t) and σ(t), the average and the standard deviation, respectively, of the number of MI-correlated particles to a tagged one at time t. They are evaluated by p(n, t). The results are shown in Figure 4. Figure 4(left) shows that the average number of MI-correlated particles n(t) reaches a peak at t τ α and drops to 2 at long times. Notably, contrary to what happens in atomic liquids [18], the height of the peak does not increase with τ α and the shapes of both n(t) and σ(t) of states with equal τ α nearly coincide [23].  particles at time t, p(n, t), for a state of the C set (same shape for both states). Moving from short to long times, i.e., from the top to the bottom of the stacked plots, the second and the fourth plot correspond to t = τ early and t = τ late , respectively, when the standard deviation σ(t) of p(n, t) is at a local maximum (see Figure 4(right)). The fractions of particles with higher correlation at these two times are highlighted. They are referred to as early and late fractions [18]. Their exact definition is given in main text. The data were taken from Reference [23]. The two peaks of the standard deviation, also reported in atomic liquids [18], identify two time scales that are referred to as τ early (orange dot) and τ late (violet dot). The data were taken from Reference [24].
Intuitively, the standard deviation σ(t) offers an alternative DH metrics. Figure 4(right) shows that it exhibits a peculiar bimodal structure similar to the one observed for the polydisperse mixture of hard spheres [18]. Following Reference [18], we define two characteristic times, τ early and τ late (τ early < τ late ), corresponding to the position of the two maxima of σ(t). In Figure 3(right), the second and the fourth plots from the top are p(n, τ early ) and p(n, τ late ), respectively. We single out the fractions of particles that are more correlated at these two times; see the highlighted regions in Figure 3(right). These fractions, henceforth referred to as early or late, respectively [18], correspond to the condition n ≥ n(t) + 2 σ(t), t = τ early , τ late .
Notably, the early and late fractions of particles show no intersection among them [18,24]. The analysis of the mobility reveals (i) that the early fraction has higher mobility and faster relaxation than the late fraction [24] and (ii) that the time positions of the two peaks of the standard deviation are simply the relaxation time of the two fractions [24]. In addition to this, it has been verified that the noted scaling between fast vibrational dynamic and relaxation observed in bulk systems (see for example Reference [38,39]) also holds for the early and late fractions [24]. Early and late fractions show a correlation with local structure [18,24]. In particular, the early population has a tendency to be located in regions of the sample with lower local density and to be arranged in unstable topological structures. On the other hand, the late one is correlated with the denser region of the sample and stable topological structures.

Clustering of Early and Late Fractions
We now focus on the spatial distribution of the particles of the early and late populations and inspect the possible tendency to group in clusters.
The cluster analysis comprises three steps: (i) first, given a single microscopic configuration of a given state of the liquid, the particles that belong to the early or late fractions are identified according to the ICE procedure outlined in Section 2.1; (ii) then, within each fraction, clusters of particles are searched according to their positions in the given configuration; and (iii) finally, suitable averages over four distinct initial configurations are performed to provide statistical significance.
To detail the second part of the cluster analysis, we start by defining a cluster as a set of M particles of the same population (M ≥ 2) with the distance from at least one other member of the cluster being less than . We choose as the first minimum of the radial distribution function, which, depending on the sample and the state, lies in the range 1.44 ≤ ≤ 1.46. Within this range, the exact value of does not affect the results. The clusters are identified by a bottom-up algorithm performing the following steps [40]:

1.
A particle of a given population is chosen and included as first member of a possible cluster; 2.
Particles of the same population within distance from the first member are searched and, in the positive case, added to the cluster; 3.
Step 2 is repeated for all new members of the cluster; 4.
If no new members are found and M ≥ 2, the cluster is completed and its size M is defined; 5.
All of the particles involved in the already identified clusters are removed, and the procedure is restarted from step 1 by considering one left particle. Figure 5 provides an example of the outcome of the first and second steps of the cluster analysis for a single initial configuration.
After completion of the cluster analysis by performing the third step, i.e., the average over the four ICEs, we start to characterise the clusters of particles belonging to the early and late populations. To this aim, we first consider the average cluster size and the corresponding gyration radius R g . The latter is defined for a cluster with size M as where r cm is the cluster center of mass. The results are plotted in Figure 6. They show a moderate increase of both the average size and radius with increasing relaxation time, which is fully consistent with the increasing DH as signalled by NGP; see Figure 2(right).
Notably, the clusters of states with equal relaxation times have equal average size and gyration radius within the errors.  Average gyration radius where c is a suitable constant. Figure 7 shows the good correlation between the size and the gyration radius of all of the clusters in all states considered in this work. A best-fit of the data by using Equation (7) yields d = 1.87(2) and d = 2.45(3) for the early and late clusters, respectively. The results suggest that the clusters of the early fractions are more filamentous, reminding us of the low-dimensional structures with correlated motion observed in glassforming liquids [41][42][43]. Conversely, late particles appear to be arranged in compact, sphere-like clusters with d = 2.45 (3). These conclusions are strengthened by limiting the fit procedure to the largest clusters with M ≥ 3 and M ≥ 5, for which One finds d = 1.96(2) and d = 2.00(1) for the early fraction, respectively, but for which, one finds d = 2.63(2) and d = 2.98(2) for the late fraction, respectively.

Mutual Information and Johari-Goldstein β-Relaxation in a Model Polymer Melt
When approaching GT, molecular rearrangements occur via both primary modes, referred to as structural or α relaxation, and via the faster secondary (β) processes, as evidenced by the mechanical, electrical, and thermal properties of the materials [44][45][46]. Although it has been the topic of a large number of phenomenological and theoretical studies as well as of experiments and simulations [46][47][48][49][50][51][52], there is still no definitive microscopic description available for the β relaxation. There is a special class of secondary relaxations involving the translation or reorientation of the molecular unit as a whole, as opposed to the trivial ones that are usually related to intramolecular degrees of freedom, such as the motion of pendant groups in polymers. This special class of β processes, called Johari-Goldstein (JG), to honor the researchers that first noticed it [47], is universal in glass formers [46,48,49] and argued to be the precursor of the structural relaxation [46,[48][49][50]. Recently, the investigation of JG relaxation using MD simulations has been facilitated by the development of coarse-grained models of linear polymers overcoming the need for complex chain architectures [53][54][55]. It is also worth mentioning the studies of the JG relaxation in asymmetric diatomic molecules [55,56].
We performed extensive MD simulations of a model linear polymer with nearly fixed bond length l 0 and bond angles constrained to 2π/3; see Figure 8(left). We considered N c = 512 linear chains made of M = 25 monomers each, resulting in a total number of monomers N = 12, 800. Further details are given elsewhere [25,26]. An interesting aspect of the model is that, even if the force field does not include a torsional interaction, an effective torsional barrier may hinder the dihedral rotations depending on the bond length l 0 due to the LJ repulsion between two non-adjacent monomers. We studied the two cases l 0 = 0.48 σ and l 0 = 0.55 σ, which present considerable or missing torsional barrier, respectively; see Figure 8(right). Henceforth, all lengths are expressed in units of σ.

Mutual Information and Bond Reorientation
Since in linear chains the JG process involves local motion of the polymer chain backbone [46,49], we focus on the most elementary relaxation process, that is, the reorientation of the bond linking two adjacent monomers of the same chain. To this aim, we consider the unit vector along the mth bond of the nth chain at time t b m,n (t) = 1 l 0 (r m,n (t) − r m+1,n (t)), (8) where r m,n (t) denotes the position of the mth monomer in the nth chain at time t. Hence, we define the bond correlation function (BCF) C(t) [57] The brackets · · · denote a suitable ensemble average. It is worth noting that BCF is a kind of Pearson correlation and then senses linear dependencies only. Figure 9 plots BCF of the two polymer melts with chains having different bond length. In agreement with previous studies [53,54], we found that the chains having shorter bond length exhibits a characteristic two-step decay after a first fast decay at t ∼ 0.1, which is a signature of the presence of two distinct relaxation processes. We ascribed the faster one to the JG β relaxation and the slower one to the structural α relaxation [25,26].

MI Correlation in Time Domain
We now show that replacing BCF from Equation (9) with the corresponding MI quantity provides sharper information on the bond reorientation and, in turn, on JG relaxation. This suggests that the JG process has significant nonlinear contributions, as anticipated in a dense glassformers [58,59]. The MI correlation of the bond reorientation is expressed by the quantity I(b(t + t 0 ), b(t 0 )), i.e., Equation (2) with X = b(t), the bond orientation at time t, and with Y = b(0), the bond orientation at the reference time t 0 . For simplicity we set t 0 = 0.  Figure 9. Temperature dependence of the BCF of the chains with different bond lengths. If l 0 = 0.48, a two-step decay-evidencing two distinct relaxations-is observed. We ascribed the faster one to the JG β relaxation and the slower one to the structural α relaxation [25,26]. The data were taken from Reference [25].  (t), b(0)) for the polymer melts with (left panel) and without (right panel) JG relaxation. In both cases, it is seen that about 50-60% of MI is lost quite fast within t ∼ 0.1, corresponding to a few picoseconds. This correlation loss is significantly higher than in the case of BCF; see Figure 9. Nonetheless, MI gives a hint of a clearer multiple decay pattern at intermediate and long times in the presence of JG relaxation. The increased resolution of MI to JG relaxation is substantiated for the state with slowest relaxation by Figure 11 comparing BCF and the normalised MI correlation function: (t), b(0))/I(b(0), b(0)) (10)

MI Correlation in the Frequency Domain
Having shown the increased resolution of MI to JG relaxation with respect to BCF in the time domain, it is legitimate to wonder how this is reflected in the frequency domain. Here, the quantity of interest is the spectral response function, or susceptibility, to be defined as Our plan is to compare χ BCF (ω) and χ MI (ω), i.e., Equation (11) with Φ BCF (t) = C(t) (Equation (9)) or with Φ MI (t) =Î(b(t), b(0)) (Equation (10)), respectively. Figure 12 compares the imaginary part of the susceptibility, χ BCF (ω) and χ MI (ω) for the systems with and without JG relaxation. It is seen that the presence of JG relaxation is weakly detectable in χ BCF (ω). Unfortunately, given the investigated time window, due to the slow decay rate of BCF, it is not possible to properly compute χ BCF (ω) down to the lowest studied temperature. Differently, a well-defined double-peaked pattern of χ MI (ω), ascribed to splitting of the two peaks associated to the α and β processes, develops when decreasing the temperature. The inset shows the Arrhenius plot of the two time scales, τ i = 1/ω i,max , i = α, β, with ω α,max and ω β,max being the local maxima of χ MI (ω), (ω α,max < ω β,max ). The maxima have been searched by best-fitting χ MI (ω) with a weighed sum of two Gaussians. On cooling, the increasing splitting of the time scales τ α and τ β of the primary and the secondary JG relaxation is apparent. χ MI (ω) (bottom) for the systems with (l 0 = 0.48) and without (l 0 = 0.55) JG relaxation. Inset: Arrhenius plot of the quantity τ i = 1/ω i,max with ω i,max , i = α, β, and the frequency of the local maxima of χ MI (ω) (ω α,max < ω β,max ). The color codes are as in Figure 9.

Dynamic Heterogeneity and Mutual Information between Displacement and Rotation of the Bond
In the polymer model under study, Figure 8, a shorter bond length yields an effective torsional barrier. We argued that the role of the barrier is twofold [26]: (i) it is responsible for a partial averaging of the DH at the JG relaxation time scale, giving rise to a peculiar double-peaked time evolution of the NGP with two local maxima located at the JG and primary relaxation time scales, respectively, and (ii) it slows down the bond reorientation and, jointly, the monomer translation, thus suggesting a roto-translation coupling, which, indeed, has been recently reported in glassforming systems with strong JG relaxation [60].
Here, we expose the roto-translation coupling by quantifying the correlation between the bond reorientation and the monomer translation. We find that the correlation is higher at the time scales of the primary and JG secondary relaxation when DH is also at a local maximum.
To provide evidence, we consider the bond reorientation in a time t, θ(t), and the corresponding squared modulus, δr 2 b (t)), of the bond displacement, δr b (t), i.e., the displacement of the center of mass of the two monomers linked by the bond; see Figure 13.
We evaluate the equal-time Pearson correlation C(cos θ(t), δr 2 b (t)) and the analogous MI correlation according to Equations (1) and (2) with X = cos θ(t) and Y = δr 2 b (t). Figure 13. Schematic representation of the angle covered by a given bond in a time interval t, θ(t), and δr b (t), the displacement of the bond position, i.e., the center of mass of the two monomers linked by the bond. Figure 14 plots the above quantities. The correlation grows with time up to the JG time scale and then decays at long times in the presence of structural relaxation. In the presence of JG relaxation, two major facts are apparent. First, MI reveals a bimodal structure peaking at the time scales of JG and α processes. The structure is quite attenuated by Pearson correlation. Second, MI correlations reach a local maximum when the NGP reaches a maximum too. The finding points to the conclusion that the primary and the secondary relaxations have in common not only the DH presence, as already noted elsewhere [26], but also the significant correlation between the bond reorientation and displacement. This suggests that, in molecular liquids, the mechanistic explanation of both DH and relaxation involves the rotation/translation coupling.  [26]. Color codes as in Figure 9.

Conclusions
We reported on extensive MD simulations concerning two coarse-grained models of dense fluids, a molecular liquid and a polymer melt. The focus was on the assessment of the additional insight offered by MI with respect to the usual Pearson correlation when dealing with two major features of transport and relaxation close to GT, namely DH and JG secondary relaxation.
In the molecular liquid, significant DH was evidenced. MI detected two distinct fractions of particles, so-called early and late fractions, which were spatially grouped in clusters with filamentous and compact structures, respectively, and exhibited different mobility and relaxation properties.
In the polymer melt, the JG relaxation was tuned by the bond length and was better revealed by MI analysis. We found that both DH and MI between the reorientation and the displacement of the bonds reach their (local) maxima at the time scales of the primary and JG secondary relaxations. This suggests that, in (macro)molecular systems, the mechanistic explanation of both DH and relaxation must involve the rotation/translation coupling.