Scaling up the lattice dynamics of amorphous materials by orders of magnitude

We generalise the non-affine theory of viscoelasticity for use with large, well-sampled systems of arbitrary chemical complexity. Having in mind predictions of mechanical and vibrational properties of amorphous systems with atomistic resolution, we propose an extension of the Kernel Polynomial Method (KPM) for the computation of the vibrational density of states (VDOS) and the eigenmodes, including the $\Gamma$-correlator of the affine force-field, which is a key ingredient of lattice-dynamic calculations of viscoelasticity. We show that the results converge well to the solution obtained by direct diagonalization (DD) of the Hessian (dynamical) matrix. As is well known, the DD approach has prohibitively high computational requirements for systems with $N=10^4$ atoms or larger. Instead, the KPM approach developed here allows one to scale up lattice dynamic calculations of real materials up to $10^6$ atoms, with a hugely more favorable (linear) scaling of computation time and memory consumption with $N$.


I. INTRODUCTION
For the case of elasticity of centrosymmetric crystalline solids, Born and Huang developed a theory which can straightforwardly predict and compute the elastic moduli from the atomistic structure 1 . Unfortunately, the task becomes considerably more complex in the case of amorphous materials which lack atomic-scale centrosymmetry. Only recently it was shown that so-called non-affine corrections to the original Born and Huang approach offer a pathway for the prediction of glass viscoelasticity 2-4 . These corrections account for additional relaxations of atomic positions in non-centrosymmetric cases and result in an overall softening of a material. In our previous work, we examined the non-affine lattice dynamics theory (NALD) against the results produced by bead-spring MD simulations for the case of polymer glasses and found excellent agreement between the two 4,5 .
However, the application of lattice dynamics calculations into the context of materials science has proven more difficult. Lattice dynamical calculations have been demonstrated as a promising path forward to relate the chemical composition of amorphous materials, as defined through atomistic models, with the full range of frequency-dependent viscoelasticity 6,7 . Similar success has been achieved for the lattice dynamics of simpler systems, such as monoatomic glasses 8,9 , granular and jammed systems [10][11][12][13][14] and topological materials 15 . Yet, two key bottlenecks remain to be solved before a broad class of "real" material compositions can be examined. First, realistic atomistic simulations often require relatively large systems on the order of 10 5 atoms or more 16,17 . A key component of lattice dynamics computations is the analysis of the spectral density of dy-namical matrices and their eigenvector characteristics or eigenmode spectrum 4,9,18 . For this relatively large number of atoms, direct diagonalization (DD) of the Hessian matrices ceases to be a viable method to obtain the eigenfrequencies and eigenmodes, since typically DD becomes prohibitive for N ≥ 10 4 due mainly to memory requirements. Thus, a method to treat larger systems must be established in order to solve the multi-scale problem in computational materials science. Fortunately, within the framework of the NALD approach, it suffices to get the distributions of eigenvalues of the Hessian matrix, which are directly linked to the vibrational density of states, rather than the exact discrete set of values. This allows the use of approximate approaches for direct computation of the vibrational densities of states (VDOS) as well as a quantity computed from the eigenvectors, the correlator of the affine force Γ. Second, the original theory 2,3 was developed for single-mass material models. The chemistry of most solids requires a multi-mass representation, and particularly so for models with atomistic detail. Hence, this paper generalizes the NALD framework for the case of multi-component solids. It will then be shown that a Kernel Polynomial Method, based on Chebyshev approximants of the eigenvectors of the Hessian matrix, can yield not only the VDOS, which was shown in previous work, but also the eigenvector-based quantities that are critical in evaluating the viscoelastic moduli of the material within the NALD approach.
All in all, the framework allows lattice dynamic calculations of viscoelastic moduli of amorphous solids to be performed, for the first time, on systems larger than N = 10 5 . From the analysis of the computational performance, it is clear that calculations for N = 10 6 are now possible. The proposed framework thus provides a working solution to the problem of bridging time and length-scales in the molecular simulation of materials mechanics. The values of storage modulus G (ω) and loss modulus G (ω) respectively, from direct diagonalisation (blue), KPM (red) and MD simulations (black squares) for the polymer system. The well equilibrated glassy system is probed at the temperature T = 0.1 T g ≈ 0.4. The total number of monomers in the system is equal to 5 × 10 3 for DD and MD, and 1 × 10 5 for KPM. More details about the simulations can be found in 20 including Refs. 5,21-24

II. VISCOELASTIC RESPONSE FROM NONAFFINE LATTICE DYNAMICS
The nonaffine lattice dynamics (NALD) approach 2 assumes that the deformation can be represented as a sum of two contributions: (1) the affine deformation typical for centrosymmetric materials, and (2) a nonaffine relaxation of atomic positions towards new equilibrium positions. The latter part produces a negative correction to the elastic free energy. Here we specialize to shear deformations, keeping in mind that the equations can be easily rewritten for other types of deformation. The expression for the free energy of deformation reads 4 where F A is the affine contribution and −F NA is the nonaffine contribution. The latter can be expressed through derivatives of net force due to affine deformation, the derivatives of the radius-vector and the shear strain amplitude (angle) γ. Under constraint of mechanical equilibrium and small deformations (γ → 0) Equation 1 can be written as 25 : where we introduce the variable Ξ i for the affine force field, defined through the force acting on atom i, f i = Ξ i γ, while H ij is the Hessian of the system (a 3N × 3N matrix), and summation over repeated indices is implied. The assumption of small deformations means that the system resides in the vicinity of a local energy minimum (effects of anharmonicity and temperature are taken into account via tension terms and negative eigenvalues of the Hessian 4,7 ). With dissipation at the molecular level, the equation of motion for a particle i of mass m can be written in damped harmonic oscillator form 2 , with inertial, dissipative and harmonic force terms on the left hand side and the affine-force field on the right side.
This key equation in the non-affine formalism has to be modified for the multi-component case. In the multicomponent case the particles or atoms can have different masses. Hence, the equation (2) can be rewritten in the following form 26 , where M is a mass matrix (N × N block matrix, where each 3 × 3 block assigns the mass m i to the particle with label i), and r(t) represents the full configuration of the system, i.e. it is a 3N -element vector. In order to solve Eq. (3), let us consider an auxiliary generalized eigenvalue problem, with φ p and ω 2 p being eigenvectors/eigenvalues correspondingly. Transforming this equation by multiplying with M −1/2 from the left gives After inserting the unit matrix I = M −1/2 M 1/2 on the right hand side where we have defined H = M −1/2 HM −1/2 andφ p = M 1/2 φ p . Now we reformulate the generalized eigenvalue problem given in Eq. (3) in a matrix form by introducing an auxiliary matrix Φ with columns being made of the eigenvectors of Eq. (4). This matrix is related to M and H as where we have defined another auxiliary diagonal matrix Ω of the eigenfrequencies, Ω = diag(ω 1 , . . . , ω 3N ). In order to solve the Eq. (3) we replace r = Φq and multiply the equation from the left with Φ T obtaining 27 with g = Φ T f being the transformed driving force. The first and the third term now can be substituted from Eq. (7) and simplified as, The second term contains the matrix product Φ T CΦ, which makes the general analytical solution of Eq. (9) impossible. This can be overcome by assuming that the damping is not correlated across different eigenmodes, i.e. Φ T CΦ is a diagonal matrix. The frictional drag force is proportional to the mass C ∝ M, which decouples the equations (9) of motion 27 . If the friction matrix C has non-zero off-diagonal elements one could approximate it with a diagonal matrix and check under which conditions the off-diagonal elements are small enough. This allows us to use index-independent notationν for the friction since (Φ T CΦ) kk = (Φ Tν MΦ) kk =ν for any k. Hence we obtain a system of decoupled equations, Applying a Fourier transform maps the equation to the frequency-space where q k , g k are the corresponding Fourier transforms of q k and g k . In Ref. 2 the general relation between the stress response ∆t η (ω) of the system to a strain η and the displacement fields r is where the summation extends over all particles. The vectors of 3N -dimensional affine force Ξ and the Fourier transform of displacement field r are functions of the driving frequency ω. In our case the second term in Eq. (12) can be transformed by using the definition r = Φq and Eq. (11), In frequency space the generalized force vector can be written as g p = j Φ T pj f j . For small deformations one can assume that the contributions of components of driving force with different frequencies are independent. Hence it is conventional to consider the case of the driving force defined as f (t) = Ξ η sin ωt 2 . The previous expression can be modified further, The matrix product i Φ T pi Ξ i = Ξ p and its transposed counterpart represent the basis transformation of the affine force field into the generalized eigenbasis. Thus, Since ∆ t η (ω) = G * (ω) η(ω) in the linear regime we get the final expression for the complex viscoelastic shear modulus of a multi-component disordered system In the thermodynamic limit, it can be rewritten as 4 , (17) where C is an integration contour which includes negative eigenvalues (imaginary frequencies), widely known as instantaneous normal modes (INMs) 4,28-30 , ω denotes the eigenfrequency as continuous variable, ρ(ω ) is the vibrational density of states (VDOS), and the correlator Γ(ω ) is defined in the following (note that in the last expression the mass dependence enters through Γ(ω )). Importantly, one can see that, in the multi-component case, the expression is very similar to the single component one 2,4 , however, the Γ(ω ) in Eq. (17) has a dimensionality difference of M −1 with respect to Γ(ω ) from Refs. 2,4 .
The results for the comparison between the VDOS computed using the Direct Diagonalization (DD) of the Hessian method and the KPM method (that will be introduced below) for the multi-mass KG polymer can be found in the Appendix F.

III. SIMULATION DETAILS
We have used the Kremer-Grest model 19 of a coarsegrained polymer system consisting of linear chains of 50 monomers which were equilibrated using LAMMPS 22 . The polymer chain under consideration consisted of two different types of masses, where the two masses were chosen as m 1 = 1 and m 2 = 3. The geometry of the chain is such that the masses are placed in alternating fashion, as illustrated in Fig.1a.
The polymer chains are embedded in a threedimensional box subject to periodic boundary conditions. In the Kremer-Grest model each constituent monomer is allowed to interact via a Lennard-Jones potential where the parameters are chosen as = 1, σ = 1. The cutoff radius of the potential is set to r c = 2.5. In addition, in the Kremer-Grest model, the covalent along-thechain bonds are represented by a finite extensible nonlinear elastic (FENE) potential given by 19 The interaction parameters of the FENE interaction are K = 30 and r 0 = 1.5. A Langevin thermostat was used for the molecular dynamics simulations where particles experience a viscous damping force proportional to the velocity. The corresponding damping constant ξ, which is related to the damping term appearing in the lattice dynamical equation of motion by ξ = mν. Using dimensionless LJ units in terms of the mass M , length d and energy , we set σ = 1 and = 1, which results in a fundamental unit of time given by τ = mσ 2 / . To apply the theory described above, one must first obtain a low-energy configuration of the solid. All of the quantities can then be extracted from this snapshot of the system and the interaction potentials. We will use the same simulation procedure as in Ref. 4 . In brief, the snapshots of the system are obtained using the LAMMPS simulation package 22 . After a sufficient number of equilibration steps in a melted state at T * = 1 the system is slowly quenched, maintaining zero external pressure using a Nose-Hoover barostat, below the glass transition temperature (T = 0.1 T g ≈ 0.4). The timescale of the cooling is τ c 10 5 τ . Ten replica configurations were constructed, and all results are averaged over these ten structures. Each glassy configuration is used as an input for the calculation of the Hessian. The latter is then diagonalized directly for comparison with the distributions of the VDOS and Γ correlator obtained by KPM. The viscoelastic moduli are also extracted from direct mechanical spectroscopy simulations 4,18 and compared then with the theoretical predictions.
For the eigenanalysis via direct diagonalisation a system consisting of N = 5000 particles. The analysis using KPM was done on a system with the same parameters but of significantly larger size of N = 1 × 10 5 particles which would be challenging for direct diagonalisation. We have checked that the results for the viscoelastic response of the system obtained with the KPM from the small and large system give similar results. Clearly, there will be slight variations in the VDOS and Γ(ω ) due to the fact that different snapshots (samples) of the glass are considered.

IV. KERNEL POLYNOMIAL METHOD TO COMPUTE THE EIGENMODES
The relevant quantity appearing in the expression for complex viscoelastic shear modulus is the product of the Γ-correlator and the VDOS ρ(ω ). Γ(ω ) is defined as a squared norm of the projection of the eigenvectors of the Hessian matrix onto the affine force field of the disordered particle system, i.e. | Ξ|p | 2 , where |p represents an eigenvector of the Hessian matrix 2 .
We will now discuss how this quantity can be evaluated by means of the KPM methodology developed here, which allows one to scale the calculation up to much larger systems than previously possible with the direct diagonalization method. We should emphasize here that the KPM does not need the assumptions of small γ and mechanical equilibrium, used in the theory. It does not care about the physical origin of the Hessian H or of the affine force field Ξ.
Since the ansatz of KPM starts with the decomposition of the VDOS expressed as a sum of δ-functions 31 , we will directly compute the product of | Ξ|p | 2 and ρ(ω ) with KPM. The Γ correlator can then be computed by dividing this quantity by the VDOS. The exact expression for the product reads 2 The KPM approximation, as shown in Appendix A, gives where the Chebyshev expansion coefficients µ k are calculated as follows.
Using the notation |u k = U k ( H)|u 0 , where U k are Chebyshev polynomials of the second kind and H is a rescaled Hessian, as shown in detail in Appendix A, we obtain that is the correct approximate Chebyshev moment which stochastically converges to µ k , i.e. m k → µ k . Here, u 0 are random vectors and the average is taken over a certain number of realizations of random vectors. These expressions are valid for a one-component system, but they have been extended here to the multi-component case. Full details of the derivation can be found in the Appendices B, C, and D.

V. KPM PERFORMANCE AND COMPARISON WITH DIRECT DIAGONALIZATION OF THE HESSIAN
The two key parameters which control the convergence and the accuracy of the KPM method are: K, which is the degree of the Chebychev polynomial at which the sum over k for J (ω ) is truncated, and R, which is the number of random vectors u 0 realizations over which the average for m k is taken. The full details about the analysis of convergence of the KPM procedure can be found in the Appendix E.
To be more precise, our realisation of KPM method has exactly R × K dot products of the sparse matrix and vector, and R×K dot products of two vectors. The speed of computation of these dot products depends on the details of the linear algebra libraries used and the sparsity of the Hessian matrix. However, each iteration in R-cycle of KPM is independent, hence KPM can be parallelized over R cycles. We apply this parallelization in the following calculations of the viscoelastic shear moduli G (ω) and G (ω) with KPM.
In Fig. 1(b)-(c), the results from MD simulations of oscillatory deformation for the polymer glass in Fig. 1(a) performed in LAMMPS are compared with the theoretical calculations using NALD with (i) ρ(ω )Γ(ω ) evaluated with direct diagonalization (DD) of static MD snapshots and (ii) evaluated with KPM. With DD we use N = 5000 while with KPM we use N = 1 × 10 5 . An excellent agreement is observed across the entire frequency range. Note that the NALD equations are fully predictive with no adjustable parameters. Even the friction parameter is taken to be identical with the friction value set in the Langevin thermostat of the MD simulations (see Section III). The successful comparison validates the derivation of multi-mass NALD above.
In the perspective of using NALD for atomistic calculations, it is important to evaluate how the computation time and the memory usage scale with the number of particles, N . Panel (a) of Figure 2 shows the dependence of the computation time for DD and KPM methods performed on multi-mass systems with different N , and KPM parameters K = 50 and R = 500. We can see that the DD method scales almost as N 3 , in contrast to KPM which exhibits linear dependence on N , with C being the number of cores, d the density of the Hessian matrix, and one should also take into account the N -dependence of d (see below). Panel (b) of Figure 2 shows the comparison of the memory usage of DD and KPM for the same systems. As expected, the memory requirement for DD is proportional to N 2 , whereas the KPM memory usage is proportional to N , or to be more exact, dN 2 ∼ N . This is mostly due to the fact that KPM uses sparse matrices and in our case the density of the matrix d is proportional to the inverse of the system size. This is the consequence of the cutoff introduced in our potential. This cutoff limits the number of interactions each atom can have to a certain number N inter (≈ 68.5 in our system). Thus, the total number of non zero Hessian elements is 9N N inter and the density d = 9N Ninter N . In our case, we have a simple system with no angular (bond-bending) or dihedral potentials ( Fig. 1(a)). In general, the N inter depends not only on the cutoff but on the complexity of the potentials. Other interactions mean just a slightly different form of Hessian, but do not change the idea or algorithm of KPM. Actually, introduction of simple angular and dihedral potentials does not even change the Hessian density (since with our cutoff these particles already interact via LJ potential, and hence the Hessian elements are already non-zero), thus the performance of the KPM remains unaltered.

VI. CONCLUSIONS
In conclusion, we have developed a new KPM-based multi-mass lattice dynamics method for computing the mechanics of real materials, which scales linearly in time with N , as opposed to the standard direct diagonalization (DD), which scales as N 2.7 . The new method is also much more efficient in terms of memory storage, with a memory consumption that scales as ∼ N as opposed to ∼ N 2 found for DD. This methodology may prove key to solve the longstanding time-scale bridging problem of atomistic simulations, which can access only the extreme high-rate (∼ 10 10 Hz) response of materials, due to the shortness of time-step. With the new method proposed here it will be possible to compute the viscoelastic response of large systems (N = 10 6 ) at atomistic or coarse-grain resolution down to deformation rates that are experimentally accessible. Furthermore, the method is general and can be applied to any solid, including perfect crystals 32 and real crystals 33 . ACKNOWLEDGMENTS A.Z. and I.K. gratefully acknowledge financial support from US Army Research Office through contract nr. W911NF-19-2-0055. Dr. Johannes Krausser is gratefully acknowledged for discussions and input during the early phase of this work.

Appendix A: Derivation of the Kernel Polynomial Method
We start by shortly describing the basics of the KPM algorithm starting from Ref. 31 (which summarizes the method as it was originally developed in the context of Fermionic particles). We consider a real-valued function f (x) on the interval [a, b]. The key idea behind the kernel polynomial approximation lies in the expansion of the function f (x) into a series of Chebyshev polynomials of the second kind U k (x) 31 , i.e.
The polynomials of the second kind U k (x) are used, because they show better convergence properties than the polynomials of the first kind 31 .
Introducing the weighted scalar product on the interval [−1, 1], we have The Chebyshev polynomials of the second kind are orthogonal with respect to the weight w(x) = π √ 1 − x 2 , i.e.
where δ k,l represents the Kronecker delta. Hence, the expansion coefficients α k appearing in Eq. (A1) are given by The Chebyshev polynomials U k (x) can also be computed using the recurrence relations or, equivalently, can be defined through their trigonometric representation Appendix B: KPM for computation of the vibrational density of states One of the first applications of the KPM algorithm in physics was the computation of the eigenfrequency spectrum of a generic Hessian matrix 9,26,31 . The vibrational density of states (VDOS) can be defined as Here and in the following we use ω to denote the eigenfrequency as a continuous variable, and ω p to denote the eigenfrequency as a discrete variable. For the KPM we have to express it as a series of Chebyshev polynomials. The function ρ(ω ) is the distribution of eigenfrequencies which result from the generic eigenvalue problem Hx p = λ p x p . Usually the matrix H represents the Hessian matrix of an interacting particle system, where the eigenvalues represent the vibrational eigenfrequencies, i.e. λ p = ω 2 p . Since the set of eigenvalues {λ p } p∈1,...,3N of the underlying 3N × 3N Hessian matrix are just the squared eigenfrequencies, we can use the variable transformation λ = ω 2 and write the DOS as In order to be able to apply the KPM algorithm, the support of the function ρ(ω ) has to be mapped onto the interval [−1, 1], to allow the expansion in terms of Chebyshev polynomials. We thus need to express the VDOS in terms of a rescaled variable λ, such that the original support of eigenvalues [λ min , λ max ] is mapped onto [−1, 1] λ. This can be achieved by a linear transformation of the eigenvalue problem, which is given by 31 where ε is a small parameter which has the function of stabilising the convergence of the kernel polynomial method against unwanted fluctuations at the edges of the support of the eigenvalue spectrum, known as Gibbs oscillations 31 . The extremal eigenvalues λ min and λ max can easily be found by standard Lanczos or Arnoldi algorithms. Using the above transformation we can express the VDOS as We now just have to expand the δ-function appearing in Eq. (B7) in terms of the Chebyshev polynomials U k ( λ).
Making use of the relation f (y)δ(x − y)dy = f (x) we can express the δ-function as 9 Using the trigonometric definitions of the Chebyshev polynomials one can write the series expansion 9 , where we have introduced the Chebyshev moments defined by The approximation then essentially consists of truncating the infinite series at a finite order, At this point the damping factor γ k has to be introduced to counteract thed Gibbs oscillations. The damping induced by γ k effectively truncates the series expansion gradually to avoid the oscillatory fluctuations which would appear if the sum were truncated abruptly 9,31 . By substituting this truncated series into the expression for the VDOS in Eq. (B7) we obtain the approximate VDOS as where K denotes the degree of the approximation which basically sets the resolution of the algorithm for approximating the δ-peaks which constitute the VDOS.
The moments µ k can be found from a modification of Eq. (B10), where |p represent normalised eigenvectors of the rescaled Hessian matrix H. The central point of the KPM is that the above trace can be approximated stochastically very accurately if the matrix H becomes very large 31 . Thus, instead of evaluating the trace over the full set of all eigenvectors, we initialise a number of normalised Gaussian random vectors |u 0 , which we want to use for the evaluation of the above trace. As an example, let us first expand one realisation of the Gaussian random vector in terms of the eigenvectors of the matrix H, i.e.
Hence, using this expansion we obtain the matrix elements which holds due to the orthonormality of the eigenvectors. The components of the random vector |u 0 in an arbitrary basis, i.e. both the components u 0,i and α p , are independently and identically distributed. They have zero expectation value and unit variance, i.e. α p = 0 and α * p α q = δ p,q /(3N ), where . . . denotes the expectation value with respect to the Gaussian probability distribution. Therefore, taking the expectation value of Eq. (B14), we obtain since the random vectors u 0 are normalised to one, i.e. |α p | 2 = 1/(3N ) 9 . Hence we can stochastically approximate the Chebyshev moments as µ k ≈ u 0 |U k ( H)|u 0 . Upon setting |u k = U k ( H)|u 0 and m k = u 0 |u k we see that, after averaging over many realisations of the random vector |u 0 , m k will converge to µ k , i.e. m k → µ k 9 . The relative error of the stochastic approximation of the trace is of the order O( √ RN ) 31 , where R is the number of random vectors drawn from the Gaussian ensemble. Therefore, starting from |u 0 , we can subsequently compute the Chebyshev moments µ k by applying the recurrence relation defining the Chebyshev polynomials. In the first iteration |u 1 is obtained by using Eq. (A5) and by applied the procedure recurrently Appendix C: KPM algorithm for the nonaffine correlator Γ(ω ) The relevant quantity appearing in the expression for complex viscoelastic shear modulus Eq. (17) is the product of the Γ-correlator and the VDOS ρ(ω ). The function Γ(ω ) is defined as a squared norm of the projection of the eigenvectors of the Hessian matrix onto the affine force field of the disordered particle system, i.e. | Ξ|p | 2 , where again |p represents an eigenvector of the Hessian matrix.
Since the ansatz of KPM starts with the decomposition of the δ-function, we will directly compute the product of | Ξ|p | 2 and ρ(ω ) with KPM. The Γ(ω ) correlator can then be computed by dividing this quantity by the VDOS. The exact expression for the product reads (C1) The KPM approximation naturally looks similar to (B9), where the expansion coefficients µ k take the form One can pull the Chebyshev polynomial U k ( λ p ) into the scalar product above in order to make use of the relation U k ( λ p )|p = U k ( H)|p (cf. Eq. (B13)): In this form we obtain a trace and can deal with it by using the stochastic approximation introduced earlier. We expand a random vector |u 0 with respect to eigenvectors of the transformed Hessian H and write the statistical average of the trace as (C5) where the asterisk denotes the complex conjugation. Since the components of the random vector fulfill α * p α q = δ p,q /(3N ), the above equation is reduced to (C6) which takes us back to Eq. (C3) and is therefore the desired result. Using the same notation as for the KPM discussion of the density of states, i.e. |u k = U k ( H)|u 0 , we conclude that the expression is the correct approximate Chebyshev moment which stochastically converges to µ k , i.e. m k → µ k .

Appendix D: Kernel Polynomial Method for multi-atom systems
Various types of particles contribute differently towards the viscoelastic properties of a solid. Hence it is useful to identify the contributions of different mass species. In this subsection we define eigenvector weight functions of different mass types and partial densities of states (pDOS).
The eigenvalue distribution can be written as a sum of delta-functions, ρ(λ) = 1/(3N ) p δ(λ − λ p ). Assuming that |p is a complete set of orthonormal eigenvectors, i.e. p i |p j = δ(p i − p j ), we can express the eigenvalue distribution as (D1) where, in the last equality, we have projected the eigenvectors of the particle basis using the projection operator I = |i i|. Thus, p|i i|p = |p i 2 |, the projection of the eigenvector onto the particle coordinate i provides the proportionality factor of the contribution of the vibrational motion of the i th degree of freedom to the full vibrational density of states. To define the generalized eigenvector weight function correctly, let us first introduce an index set M n for each mass type. The set M n denotes the set of labels of the particles with mass type n. Expressing the norm of a generalized eigenvector φ p in terms of the eigenvectorsφ p , we have where φ p |i =φ p,i represents the i th component of the eigenvectorφ p . Hence, in order to define the correct weight function in terms of the generalized eigenvectors, we need to normalise the contributions φ p |i i|φ p with φ p |φ p . If one uses projections on the vector components which belong to the different mass species given by the index set M n , n P n = I, the Eq. (D2) can be written as Thus, for instance, the correct weight function of mass species 1 reads The generalized frequency-dependent weight functions sum up to unity, i.e., n i∈Mn This gives a method for splitting the full density of states (VDOS) into different mass contributions or partial density of states (pDOS) In order to compute χ n (ω ) with KPM, we have to modify the scheme used in the single-mass density of states.
In our context, the starting point for the KPM always involves summations over δ-peaks. In the case of the quantities χ n (ω ) and Γ(ω ), the sum contains an additional weighting factor depending on eigenvectors to expand and approximate the desired function. The weight has the numerator φ p |P n |φ p and the denominator φ p |φ p .
To implement this idea, we start with the numerator in Eq. (D4) and define an auxiliary function χ n (ω ) 26 which is amenable to the expansion in terms of Chebyshev polynomials. Performing the analogous steps as for the single-mass VDOS above, we now map the support of the eigenvalue spectrum of the generalized eigenvalue problem onto the interval [−1+ε/2, 1−ε/2] and continue with where the corresponding Chebyshev moments are now To pull the polynomials U k ( λ p ) inside the scalar product we now have to be careful since we are dealing with generalized eigenvectors. Only the eigenvector components of the mass species n are projected out. Hence, it is possible to write the term φ p |i i|φ p appearing in the above summation as 1 mn φ p |i i|φ p , because the vectors i are eigenvectors of the mass matrix M with eigenvalue m n . This allows us to use the relation U k ( λ p )|φ p = U k ( H)|φ p and obtain by reabsorbing the factor m n into the generalized eigenvector, where H is a rescaled Hessian H as given in Eq.
B3. Subsequently, making use of the stochastic evaluation of this trace with the Gaussian random vectors u 0 , the Chebyshev moments in Eq. (D9) are approximated by averaging the quantity where P n represents the projector of the particle species n. The approximate Chebyshev moments converge to the actual Chebyshev moments µ k , i.e. m k → µ k . Due to the fact that random vectors u k are supposed to represent the generalized eigenvector of the multi-component system, it would be incorrect to use normalised Gaussian random vectors as before.
To achieve the correct stochastic approximation of the Chebyshev moments µ k , we multiply a normalised random vector ξ 0 by the inverse square root of the mass matrix M. As a result, the initial random seed of the KPM algorithm in this case is the random vector u 0 = M −1/2 ξ 0 . The same reasoning is applicable to the denominator of Eq. (D4). The first step is an auxiliary function given by where we use the subscript to signal that this is the KPM approximation function for the normalisation factor of the weight function χ n (ω ). Going through the same steps as for χ n (ω ), the final result in terms of the associated approximate Chebyshev moment is which converges to the true Chebyshev moments appearing in the expansion of χ norm (ω )  Figure 3: Convergence of the KPM with the number of random vector, the difference with the reference 5k system is plottted. K = 1000 for both DOS and Γ(ω )ρ(ω ).
in the statistical average as m k → µ k . Having set up the KPM approximation for these two components we subsequently obtain the weight functions χ n (ω ) as the ratio of the two converged auxiliary functions, i.e.
which is defined on the support of the eigenvalue spectrum of the Hessian matrix with the condition that χ norm (ω ) = 0. denotes the reference values of the exact solution. Using the 5k system as the reference KPM, Figure 3 shows the convergence of the KPM with the number of random vectors. We note that J (ω ) has much slower convergence than ρ(ω ). Hence, a significantly larger number R of sample random vectors has to be drawn in order to achieve a good approximation. The cause for this difficulty stems from the fact that the random vector used in the approximation of the Chebyshev moments µ k is projected on the affine force field vector Ξ, which itself is an inherently random quantity due to the structural disorder of the polymer configuration.
As a consequence, larger fluctuations occur which need more iterations to be smoothed out. To achieve a good approximation using the KPM for the VDOS, usually between 10-100 averaging iterations are required. In the case of the non-affine correlator Γ(ω ) estimation, however, between 10 3 to 10 4 iterations are needed to converge the algorithm to a reasonable degree, depending also on the desired resolution.  polymer system consisting of N = 1 × 10 5 particles. In addition to the large Lennard-Jones peak at lowfrequencies, we notice that in comparison to the shape of the VDOS of the single-mass polymer system 23 , the highfrequency FENE peak has split into two smaller peaks. In Fig. 4d we show how the two partial densities of states sum up to the full VDOS. It is interesting to observe that this double peak is comprised almost exclusively of modes from the lightest masses m 1 = 1, likely representing fast oscillations of the m 1 with respect to m 2 , which is a factor of three heavier (see also 23 for the discussion of the physical origin of the VDOS peaks). Figure 6 shows the product Γ(ω )ρ(ω ), which is the direct output of the KPM algorithm, obtained for the multi-component non-affine correlator. Again the agreement between the direct diagonalization and KPM results is excellent.
Finally, we use Γ(ω )ρ(ω ) to calculate the components of the complex viscoelastic shear modulus obtained by different methods in Figs. 1b-c of the main text. Following 7 we introduced here a low-frequency cutoff ω cut = 1, i.e. the frequencies with absolute value ω ≤ ω cut were excluded from the integration in Eq. (17). This cutoff eliminates spurious contributions from poorly sampled low-frequency regions, which otherwise lead to large uncertainties in the low-frequency modulus. The exact value can be estimated as ω min = cs L0 , where c s is the shear wave sound speed and L 0 is the box length. It also means that since KPM give possibility to approach larger system it also able to sample lower ω . The agreement between KPM and direct diagonalisation is excellent. The MD results also show excellent agreement with the DD/KPM results for G and G .