Characterizations of Network Structures Using Eigenmode Analysis

We introduced an analysis to identify structural characterization of two-dimensional regular and amorphous networks. The analysis was shown to be reliable to determine the global network rigidity and can also identify local floppy regions in the mixture of rigid and floppy regions. The eigenmode analysis explores the structural properties of various networks determined by eigenvalue spectra. It is useful to determine the general structural stability of networks that the traditional Maxwell counting scheme based on the statistics of nodes (degrees of freedom) and bonds (constraints) does not provide. A visual characterization scheme was introduced to examine the local structure characterization of the networks. The eigenmode analysis is under development for various practical applications on more general network structures characterized by coordination numbers and nodal connectivity such as graphenes and proteins.


Introduction
In nature, all materials can be considered a collection of atoms in specific patterns, and their properties are determined by the interactions between the atoms.The structures of materials, even in continuous limits, can be approximated by discrete nodes and bonds [1,2].A simple model to describe these properties in the set of many atoms and their interactions would be a network model.Research on the network structures has a long history mostly back to the Lagrange era for mechanical properties [3].The analogy of network structures to real materials in small scale has been widely employed to address OPEN ACCESS material characterizations both in analytic and computational approaches [3][4][5][6].Mathematical formalisms in many different forms have been developed to predict statistics and kinematics of the networks by solving linear sets of equations of motion [7].One of the fundamental mathematical analyses has been developed in the form of spectral graph theory to explore the topological characteristics of matrices constructed by nodal connectivity [8,9].In the network structures, topological connectivity between each structural unit can construct a Laplacian matrix, which includes intrinsic structural information.The spectral analysis on the networks can provide structural stability by determining vibrational modes by matrix diagonalization.Afterwards, many other methods of the eigenmode analysis have been developed originating from the spectral graph theory.However, network analyses in distinct details and application fields shared similarity of their mathematical frameworks, which is originated from the spectral graph theory.
The network models have been successfully employed in various characterizations of materials in nanoscale, and have been investigated due to interest in their applications in mechanical, electrical, and biological devices.It has been shown that the network models in specific atomic arrangements either in certain lattice or in amorphous can reflect most of the nanomaterials in infinitesimal variation limits.The material properties such as structural stability and energetics can be dependent on the specific network structures [10].For instance, covalent bonds in atomic structures can be strong enough to maintain topology in low temperature regime.Graphenes or proteins can be typical network objects on which the spectral graph theory can be applied to address structural characteristics.One of the widely-known network models would be Gaussian network model (GNM), which decomposes normal modes by the diagonalization of Laplacian (Kirchhoff) matrix.Residual fluctuations and correlations of biomolecules such as B-factor (or Debye-Waller factor) of proteins can be determined by GNM and shown consistent with the X-ray crystallography [11].It has been generalized to be an isotropic network model (ANM) to determine three-dimensional characteristics of directional preferences [12].A computationally efficient model called flexibility and rigidity index (FRI) has been recently proposed to determine localized properties of the networks without performing full eigenmode decomposition [13].The FRI method utilizes only local geometries to determine localized properties such as B-factors in proteins and chemical reactivity in atomic molecules.By utilizing geometry of networks and specific environment, the FRI method is considered highly suitable and efficient for the prediction of localized biological/chemical functions and stabilities.
In this study, we employed a full eigenvalue decomposition method to address the structural characterization of the elastic networks of nanomaterials, such as graphenes not limited to specific biomolecules like proteins, and compared with previous rigidity percolation studies [5,6].From the numerical analysis, amplitude variations of nodes for eigenvalues (frequencies) were obtained for various lattices and amorphous networks with defects.Due to its intrinsic dynamic nature, the eigenmode analysis provides structural stabilities of the networks in dynamic environment.It was applied to determine global rigidity characterization of the networks, and also to identify local structural characteristics such as floppy and rigid regions.Pristine lattices and inhomogeneous networks with void, topological defects, and bond removals were numerically analyzed to obtain variations of the eigenmode spectra.Using the eigenvalue spectra, structural stability of the networks was visually characterized.Similarly, vibration energy of specific local regions was successfully identified to predict actively responsive to external dynamic inputs.The eigenmode analysis provided insight on geometric effects on the vibrational density of states and structural stability in global scale and local scale as well.

Rigidity Analysis
One of the genius characterizations of network stability would be Maxwell counting to determine rigidity of networks [4].It is based on a simple constraint counting on nodes (degrees of freedom) and bonds (constraints).Extra zero frequencies other than translational and rotational motions of whole body can allow virtual motion without energy variation, which corresponds to floppy modes [5].Rigidity percolation study is a fascinating subject providing physical understanding on continuous random network models for glasses by relatively simple geometric models [2,5].As numerical capability to solve large network models is enhanced, various numeric algorithms replace the analytic approach to investigate generalized network models (with defects, amorphous, and complex connectivity).One of the recently developed, simple and fast combinatorial algorithms would be a pebble algorithm originated from ancient children's play [14].It was applied to large central-force and bond-bending networks to determine rigidity percolation characteristics and to identify locally floppy, rigid, over-constraint regions.It was shown to agree well with neutron scattering measurement of the density of states [6].The structural stability sometimes had been investigated in dynamical point of view by considering realistic interactions between the nodes [15].It can provide more versatile characterization tools, e.g., assessment of network rigidity and mechanical strength which the Maxwell's constraint counting cannot provide.For rigorous study, until recently, it still requires high computational cost because enormous large set of linear equations should be resolved.
One characteristic of the rigidity of networks would be mechanical strength to resist under certain deformation.A network structure made of nodes and bonds can be characterized by certain geometric parameters, e.g., coordination number and nodal connectivity.Rigidity percolation is a phenomenon that cohesive rigidity occurs for a structure beyond a certain threshold of these parameters.A general network structure can have floppy, rigid, and over-constraint regions [16].Nodes in the floppy region can freely move without any energy cost in infinitesimal displacement.If the infinitesimal displacement of the nodes changes the energy, it is called rigid.In an over-constraint region, there are more constraints to be marginally rigid.
The network shown in Figure 1 is globally rigid by the Maxwell counting rule with more constraints than the degrees of freedom, but it is noted that the network can have a locally floppy region.Two nodes in the red dotted box are only in the floppy region and the rest nodes in the rest part of the network are completely rigid.It is shown that one can remove a bond without destroying the rigidity of the whole network.For example, the removal of a diagonal bond connecting A and B does not change the structural property of the network.Thus, this bond is an excessive bond in the over-constraint region.It should be noted that the Maxwell counting rule is simply to determine the rigidity, but it does not provide any detail information of locally floppy, rigid, and over-constraint.
Suppose that one may add a bond between C and D, and an over-constraint bond at the center is removed.We can have a completely rigid structure without any locally floppy region.This new network structure has the same number of bonds as the original network.There are four zero-frequency modes for the structure deformation without any energy cost.Two transitional modes and one rotational mode are all trivial modes in two-dimensional networks with open boundary conditions.The one nontrivial zero-frequency mode is determined that for which the floppy region moves horizontally with infinitesimal amount.Here we regard the zero frequency mode is the case where the longitudinal length is not changed for all bonds in the first order approximation.

Methodology: Eigenmode Analysis
We employed an eigenmode analysis to investigate structural properties on various types of networks [7,17].A set of linear Lagrange equations induced for various network structures was constructed and implemented for computational analysis of their eigenmode spectra.The networks consist of nodes with mass and bonds by Hook's law subject to infinitesimal deformation limits [15].It was applied to analyze the structural properties and energetics of networks.

Equation of Motion
The eigenmode analysis has been introduced for the study of normal mode properties of some two-dimensional network structures.For this, we have derived a set of linear equations for small oscillation of d × n coordinates where n is the number of nodes and d is the spatial dimension.In two dimensions (d = 2), we establish a 2n × 2n matrix with the coefficients of these linear equations and, hence, the normal mode vibrations were investigated in two-dimensional networks.Force equation for i-th node exerted by j-th node for the updated positions ( ⃗  ,  ⃗  ) of two nodes from the initial positions ( ⃗ ,0 ,  ⃗ ,0 ) is given by where   is the spring constant of the bond between i-th and j-th node.For small displacements of  ⃗  =  ⃗  −  ⃗ ,0 and some derivations under the first order approximation, this becomes where  ⃗⃗ ,0 =  ⃗ ,0 −  ⃗ ,0 and  ⃗  =  ⃗  −  ⃗  .All nodes in a normal mode are subject to vibrate with the same resonant frequency ω from the equilibrium positions as the following (only x-component is under investigation and y-component is in the same form) Then, the equation of motion of the x-component of the i-th node (   ) surrounded by the j-th nodes is represented by Here,   is the mass of the i-th node and all   are considered equal as M. For the spring constant,   =   in general.The matrix form of the linear equations is represented by where Here, Myy and Myx can be determined in a similar manner with Mxx and Mxy.The sum of ∑ Element  〈〉 j for i-th node is taken for all connected node (j).In the matrices η  is 1 when two nodes of i and j are connected and 0 otherwise.These matrices are sparse because only the nearest neighbors are connected.
Solutions of the above Equation ( 5) can describe entire normal mode motions of the network system.Eigenmodes with zero eigenvalue, equivalently zero frequency, are those which require no energy for the vibrational motions.Under open boundary conditions, there are three trivial modes, two translational modes and one rotational mode.Under a periodic boundary condition, the rotational mode does not exist.If there are zero energy modes other than these trivial modes, they are floppy modes of the network structure and show that the network is floppy.Here, it should be noted that if the network contains a single locally floppy region, the analysis generates zero energy mode and the network is considered floppy.Thus, it may not distinguish global rigidity and local rigidity.However, as we observed in Figure 1, it is possible for globally rigid network to possess a locally floppy region.Since the global rigidity (or structural stability) might be crucial in many real structural applications, one leaves this study for an important future work.
Figure 2a shows a triangular network with randomly removed bonds.Before one removes any bond, the perfect triangular network is completely rigid, and, hence, the number of nontrivial zero energy mode is zero.One may expect that as more bonds are removed, the number of nontrivial zero energy modes increases.If nearly all bonds are removed, that is, all nodes are nearly independent, the number of nontrivial zero energy modes reaches its maximum value and this result is shown in Figure 2b.Obviously, this agrees with the results of previous studies including the pebble game algorithm.The value of bond impurity fraction, at which the second derivative of the nonzero energy mode fraction with respect to the bond impurity fraction reaches the maximum value, is the threshold of the rigidity percolation for this structure.The value of the fraction of connected bonds is reported as 0.66 from the previous studies [6], and the result from our eigenmode analysis on relatively small sized samples is shown to be well matched with the prediction.

Nodal Amplitude
Normal mode behavior of a network structure is well represented by the eigenvalue distribution.Since the eigenvalue is a square of the normal mode frequency, the eigenvalue distribution shows some aspects such as the bandwidth, that is, the maximum frequency and the probability distribution of the frequencies, which is related to the degeneracy of the specific eigenvalue.As in our previous studies [7,17], the eigenvalue distribution, that is, the vibrational density of states can be represented by a continuous form as a sum of Gaussian functions as Equation (8).The width of the Gaussian function σ is arbitrarily chosen and each Gaussian function represents a single eigenmode and, therefore, eigenvalues with large degeneracies are represented with high peaks.
For each normal mode, the nodes in the network structure vibrate in some patterns and with certain amplitudes.In order to represent a normal mode in two-dimensional structure with n nodes, we have 2n amplitudes.Therefore, for a given eigenvalue λ and node j, there are two amplitudes in two directions, A , (λ) and A , (λ).The two amplitudes in each axis can be combined to a radial amplitude A , (λ) as follows: This radial amplitude A , (λ) is the level of the motion of j-th node at the normal mode with the eigenvalue λ.In order to represent this in a continuous form, we multiply A , (λ) by the Gaussian function g(λ, σ) as While the above represents the response of the j-th node to the normal mode with certain eigenvalue λ, one may want the response of the node to many normal modes with the eigenvalues within a certain domain λ a ~λb .We just sum the value in the Equation ( 10) over this eigenvalue domain as follows: This is a form of amplitude spectrum of the j-th node, which represents how the level of the motion of j-th node changes according to the eigenvalue in this domain of the eigenvalue.For the full amplitude spectrum, we set

Structural Analysis
The radial amplitude can be simultaneously determined for all nodes in the network structure.By calculating and showing a scalar value for all nodes, one can show certain properties for the structure.One can add up squared values of radial amplitude for an eigenvalue domain as follows: This value F  (λ 1 , λ 2 ) can be calculated for all node j and be shown together in order to represent node responses for the eigenvalue range of λ 1 ~λ2 .We will call this quantity as a squared amplitude sum.
With this squared amplitude sum, we can identify individual node or set of the nodes which are more reactive under external dynamic inputs with a certain range of eigenvalue.

Triangular Network with Bond Defects
In a network structure, the floppy region costs no energy to be deformed while others cost a certain amount of energy.If we choose the eigenvalue domain only for zero eigenvalue, we may find the floppy regions with the squared amplitude sum.Other than the trivial modes such as translational modes, only nodes in floppy regions may move for zero eigenvalues, that is, for zero energy modes.In general, there are degeneracies for zero eigenvalue and there are many zero energy modes other than the trivial modes for large networks.Therefore, the squared amplitude sum for zero eigenvalue has much larger value for nodes in floppy regions than for the rests.Figure 3 shows the triangular network with randomly removed bonds.All nodes in the network have the same mass and the spring constant is uniform.In the figure, the size of the node denotes the squared amplitude sum for zero eigenvalue.Only two nodes show large amplitude for the zero energy modes representing floppy nodes.This method can be applied not only for zero energy modes but also for the finite energy modes.By calculating this squared amplitude sum for various eigenvalue domains, one can find structural properties of a network regarding normal mode dynamics.

Randomized Hexagonal Networks
As mentioned earlier, the squared amplitude sum for higher eigenvalue domains enables us to determine stiffer regions.Nodes in these regions are more responsive for normal modes with higher eigenvalues.Here we introduce a hexagonal lattice network as an illustrative example.Two-dimensional hexagonal network is geometrically identical to the graphene in real world [18].Numerous characterizations were performed on the superior properties of the graphenes with realistic defects.One of the typical defects in graphene is the Stone-Wales defect [19,20].Figure 4a shows a structure of 800 node hexagonal network with five Stone-Wales defects.Here, we calculated the amplitude sum for the eigenvalues larger than 3.0, F  (λ > 3) and showed as the color for all nodes.Nodes with large values of F  (λ > 3) are mostly in the Stone-Wales defects and nearly all other nodes have very small values of F  (λ > 3) except for a few.In addition, these few nodes are all rather close to the defects.The reason to choose the eigenvalues above 3.0 is that adding these Stone-Wales defects to the hexagonal network creates eigenvalues above 3.0 as shown in Figure 4c. Figure 4c represents the vibrational density of states for various networks plotted with the sum of the Gaussian function g(λ, σ) with σ = 0.05 according to Equations ( 8) and ( 12) [17].By adding more Stone-Wales defects and relaxing with molecular dynamics, one can get a random network as shown in Figure 4b [1,5,17] and here we also showed the squared amplitude sum F  (λ > 3).Nodes with large values of F  (λ > 3) are distributed in an irregular pattern; however, we find these are located near the regions with more deformed local structure.Especially, nodes in the square box have highest values of F  (λ > 3) and three pentagons are neighboring here, that is, this is the mostly deformed region.The center point where these three pentagons meet has a very small value of F  (λ > 3).Normal modes with λ > 3 create radial motions of nodes in each pentagon, and the destructive interference of these three radial motions at a given frequency destroyed the amplitude of the center node.From Figure 4c, one can see that there are a lot more eigenvalues larger than 3.0 for the random network than for the network with five Stone-Wales defects.This is consistent with the fact that there are more deformed regions and more nodes are involved with the motions with the eigenvalues larger than 3.0.

Energy Analysis
While the squared amplitude sum describes the dynamic property of a network for a specified eigenvalue range, one may need a more universal quantity obtained through whole eigenvalue domain to describe overall network properties.One can try to add up the amplitude sum for all eigenvalues λ 1 to λ 2 , that is, F  (λ 1 , λ 2 ).From Equations ( 9) and ( 12) And this is a sum of squared amplitudes of two orthonormal eigenvectors and hence is just two in any case.
In the above Equation ( 15), λ is multiplied to the squared amplitudes before we sum up.We will consider this quantity as an averaged eigenvalue.For a normal mode in a system with uniform mass for all nodes,  = 1 and λ = ω 2 = .Then, the above has a form of ∑ A , 2 (λ)  and has a meaning of a sum of energy stored in the springs connected to the j-th node.Unless the network is uniform both for the coordination number and the spring constant, this quantity C  (λ 1 , λ 2 ) is nontrivial.By calculating this value for all nodes, one can show an energy contour to describe certain network properties.

Triangular Network with Bond Removals
Figure 5 shows a triangular network with randomly removed bonds and the averaged eigenvalue C  (λ 1 , λ 2 ) denoted for all nodes.This network has uniform spring constants, and the nodes have various coordination numbers varying from one to six.As one might conjecture, the averaged eigenvalue is the coordination number in this case since C  (λ 1 , λ 2 ) has a sense of sum of spring energy connected to a j-th node mention above.

Hexagonal Networks
Figure 6 shows a hexagonal network with inhomogeneous spring constants with Gaussian distribution / = 0.1.The homogeneous spring constant case is rather trivial, that is, the averaged eigenvalue C  (λ 1 , λ 2 ) would be three for all nodes, but this inhomogeneous case is different.Nodes connected with springs with higher spring constant have higher values of C  (λ 1 , λ 2 ).Regions containing these nodes  For these extreme nodes A and B, Figure 7 shows the amplitude spectra F  (λ 1 , λ 2 ) in Equation ( 12).F A (λ 1 , λ 2 ) has a large value for high eigenvalues and a small value for middle and low eigenvalues.F B (λ 1 , λ 2 ) has very small value for high eigenvalues and large value for middle and low eigenvalues.
The amplitude spectra F  (λ 1 , λ 2 ) represent physically more active frequencies and less active frequencies for a specific node.12) of the hexagonal lattice for nodes with highest and lowest values of C  (λ 1 , λ 2 ).

Conclusions
In this study, we have various network properties by using the eigenmode analysis.A full eigenvalue decomposition scheme was introduced to explore structural information of the networks in a wide range of normal modes.This method was shown to determine and visualize various characteristics of normal modes for specific frequency ranges or a full range of spectrums.By using the amplitude analysis, we were able to identify floppy regions in a network without energy cost under infinitesimal deformations.For this, we searched those regions containing nodes with large amplitudes for the zero eigenvalue.By searching the same for higher eigenvalues, more deformed regions were visually identified for the Stone-Wales defects added on the hexagonal network and the random network.In general, with this structural analysis, we can find dynamic response behavior of the networks for normal modes in any frequency domain.Along with the change of the network structure such as adding structural impurities and random networks, inhomogeneity of the spring constant is also a possible generalization for the network problem.The energy analysis was utilized for this problem and hence an energy contour map was visualized.
It is clear that our eigenmode analysis may be implemented to calculate the vibrational density of states for a given system, which is comparable to the phonon spectrum by experiments such as Raman spectroscopy and other calculation methods [21].Introducing longer-range bonds and angular force terms may be possible along with various structures with different masses to simulate real materials.Unlike the quantum mechanical calculations, our method has the capability to explore much larger sizes of the system.This is because one can solve the eigenvalue problem for tens of thousands of matrix dimensions even without using the sparse matrix technique [22].The amplitude spectrum may replace the vibrational density of states in order to investigate the locality dependence of the structural property.Hence, this enables us to predict certain results of some other methods.

Figure 1 .
Figure 1.Square lattice with diagonal bonds added with open boundary condition.Nodes in the dotted box denote a floppy region.

Figure 2 .
Figure 2. (a) Triangular network with randomly removed bonds.The dotted box denotes the unit box for periodic boundary condition.(b) Rigidity percolation behavior is shown in the plot for the fraction of nonzero energy mode vs. bond impurity fraction.

Figure 3 .
Figure 3. Triangular network with randomly removed bonds.Dotted box denotes the unit box for periodic boundary condition.Size of the node denotes amplitude for zero energy modes, that is, floppy nodes.Note two floppy nodes have longer radii than the other rigid nodes.

Figure 4 .
Figure 4. (a) A hexagonal network with five Stone-Wales defects.(b) Nearly random network of the hexagonal network with 65 Stone-Wales defects.Color of the node denotes sum of the amplitude for eigenvalue larger than 3.0.(c) Vibrational density of states of hexagonal networks: pristine graphene, five Stone-Wales defects, and random network.

Figure 5 .
Figure 5. Triangular network with randomly removed bonds.Dotted box denotes the unit box for periodic boundary condition.Color of the node denotes the averaged eigenvalue.
are stiffer and have more mechanical energy stored.Similarly, nodes connected with springs with lower spring constant have lower values of C  (λ 1 , λ 2 ).Node A and B have highest and lowest values of C  (λ 1 , λ 2 ), respectively.

Figure 6 .
Figure 6.Two-dimensional Hexagonal lattice with inhomogeneous spring constant, Δk/k = 0.1.Dotted box denotes the unit box for periodic boundary condition.Color of the bond denotes value of the spring constant and color of the node denotes the averaged eigenvalue.