Boson peak in disordered materials under shear deformation

During a shear process the vibrational mode structure of a non‐crystalline model material will change under load. Thus, we expect an effect on the characteristic boson peak, which correlates with numerous features of disordered materials. In this paper, we perform shear deformation on two‐dimensional random network materials and investigate the distribution of their vibrational density of states (VDOS). Furthermore, the spectra of eigenvalues are studied in detail using similar approaches to investigate the eigenvectors and specifically their change due to load and plastic rearrangements.


INTRODUCTION
The boson peak is ubiquitous in disordered materials [1].It is expected to be deeply related to the glassiness of materials [2] and appears as an excess of modes in the low-energy regime relative to the Debye theory [3].Although extensive research has been done during the last two centuries regarding its origin and composition [4][5][6][7][8][9], the reason and the nature of these modes is still a topic of discussion.Current theories include descriptions using 1D stringlets [10,11], higher order excitations [4,12], broadening of the van Hove singularity [9] and fluctuations in the local elasticity [13].In this paper, we want to contribute to the comprehension of this complex problem by performing athermal quasistatic shear (AQS) simulations on vitreous, silica-like model materials in order to understand the evolution of this phenomenon during mechanical load.Additionally, we will briefly investigate the influence of disorder on the modal structure of vitreous materials.

TWO-DIMENSIONAL DISORDERED NETWORK MATERIAL
The chosen model material is a representation of two-dimensional vitreous silica.It was manufactured and investigated using transmission electron microscopy (TEM) by Huang et al. [14].In the crystalline state it looks similar to graphene.However, the binary structure and isolating properties are major differences that define the modal structure and influence the fracture behavior.The interaction between all particles is evaluated using a Yukawa-type potential [15]: F I G U R E 1 Bond switch process in a crystalline environment.
In this equation, (  ) denotes the pair potential,   is the distance between atom  and atom , and   ,   and  are fitting parameters taken from Roy et al. [16].

Generation
The generation of vitreous silica is based on the assumption, that the origin of disorder is represented by random bond switches in the structure.Following this idea leads to the Monte Carlo random bond switching algorithm proposed by Morley [17].Starting from an initially crystalline structure, we iteratively induce random bond switches to the dual lattice.Two connected silicon atoms with their neighboring oxygen atoms are rotated.The outer four bonds are detached and reconnected to the next unbound silicon atom, as shown in Figure 1.This procedure takes four neighboring rings and changes their ring sizes, such that the two rings placed in the row of the two connected central silicon atoms become larger by one, whereas the size of the two rings at the side of the two silicon atoms is reduced by one.Observation of the configurations generated by totally random bond switches revealed that the configuration does not represent a physically meaningful structure [14,18].Further constraints are needed to prevent clustering of rings of similar sizes, leading to the introduction of the Aboav-Weaire law [19][20][21].Their law introduces a condition for the ring size distribution under which bond switches are accepted or rejected: In Equation ( 2),   denotes the mean ring size around all rings of size ,   is the mean of the ring statistics and  2  is the variance.The parameters   are taken from ref. [22].The difference   between target ring statistics,  () 4 = 0.31 and  () 5 = 0.36, and current ring statistics,  4 and  5 , is evaluated as: This equation is used as the objective function for the Monte-Carlo bond switch algorithm.The coefficients  1 are  2 weighting factors between 0.1 and 10.0 [23] and are set to 2.0 in this work.The target values for the parameters  4 and  5 are chosen to approximate the left hand side of the ring statistics (ring sizes smaller than six) of a physically meaningful configuration.To prevent getting stuck in a local minimum while minimizing   , upwards steps from switch  to  + 1 Clip of true shear deformation protocol for a configuration with heterogeneity of 1.00.are accepted with a probability P using the Metropolis condition, written as: In this equation, the energy equivalent T is chosen to be 10 −4 to achieve satisfying results.

Deformation
We use athermal quasistatic simulation protocol as deformation method to apply true shear to the samples.This shear protocol is fully deterministic and allows to investigate isolated transitions between mechanically stable states [24].First, affine deformation is applied to each atom.This results in changes of the potential energy landscape and the state of the system.In most cases this will result in a non-equilibrium state for non-crystalline materials.From this state, the system is transferred to the closest local minimum using the conjugate gradient method.The trajectory to the minimum after the affine deformation shows the nonaffine displacement field.In order to keep the shear stress state consistent, the sample is relaxed every five timesteps.During the relaxation, the stress is minimized along the  11 ,  22 , and  12 direction.True shear is implemented by elongating the box in one direction and shortening the box in the other direction, such that its area remains constant, as shown in Figure 2. The step size is set to 5 Å and the maximum number of steps is 500.The nonaffine displacement contains information about the material failure [25][26][27] and is derived from the potential energy landscape using the approach from Maloney et al. [24,28], written as: In this equation, •   ∕ denotes the displacement field, and is the inverse of the Hessian, as defined in Equation ( 6) (the Goldstone modes are excluded to ensure invertibility) and   represents the affine deformation.The quantity  is the shear angle.The Hessian is evaluated as: In this equation,   denotes the force acting on atom .The stress-strain curve is determined to calculate the elastic and shear modulus and to find the position for the first stress drop, as shown in Figure 3.

F I G U R E 3
Stress strain curve of a sample with heterogeneity 0.60.The stress is given for the horizontal direction, the strain is given in terms of the shear angle .

Vibrational Density of States
The vibrational density of states (VDOS) is a characteristic of a vibratory system defined via Equation (9).It can be used to visualize the distribution of modes along the frequency axis.Beginning from the potential of an arbitrary model material, we calculate the VDOS by determining the eigenvalues of the corresponding Hessian of the system.While an analytical approach over the integration of the state function is not possible for disordered materials, the concept of AQS allows us to directly solve the eigenproblem of the dynamical matrix   , which can be derived from the Hessian by: The eigenproblem with the eigenvalues   , eigenfrequencies   and eigenvectors   is defined as: The modes linked to each eigenvalue are investigated in detail for crystalline materials.However, in disordered materials the excess of modes over the Debye approximation in the low energy regime is a most intriguing observation that has been the focus of extensive research during the last years.The investigation of the distribution of modes along the energy axis, is evaluated as: Here, () is the VDOS,  is the number of degrees of freedom,  is the frequency and   is the −th eigenfrequency of the system.In our case, we have relative units and therefore calibrate the -axis, since we use a dimensionless potential.We chose a calibrations referring to the van Hove singularity of the system.Additionally, every vibratory system can be associated with a Debye frequency of states, which was evaluated as: In Equation (10),  is the number of atoms.The shear modulus  was determined using the AQS protocol and a linear fit as an approximation to the first 20 deformation steps in the elastic regime of the stress strain curve, exemplary shown in Figure 3 for one sample with a heterogeneity of 0.60.Further,  is the total mass of the system.The Debye frequency is:

RESULTS
Performing the AQS simulations described in Section 2.2 gives an exemplary stress-strain curve, as shown in Figure 3.
The VDOS described in Section 2.3 is, for comparison, divided by , yielding to the reduced VDOS.The reduced VDOS is depicted for the crystalline configuration and two disordered configurations in Figure 4.
Looking at this figure reveals that the crystalline configuration has higher and sharper peaks in its reduced VDOS than the configurations with disorder.Furthermore, it seems that the amount of disorder in the system broadens the spectrum and flattens out peaks, meaning that eigenfrequencies, which are clustered around a certain frequency, spread apart as a response to the increase randomness in the force coupling constants.
Having discussed the characteristics of the reduced VDOS for various heterogeneities in an initial state, we proceed to the investigation of the reduced VDOS during the shear process until the break.These results are shown in Figure 5.This figure shows from top to bottom the reduced VDOS of samples with increasing heterogeneity.The left side shows the evolution of the reduced VDOS over the whole shear process until the first rearrangement.The red plot on the left side denotes the reduced VDOS immediately before the first rearrangement.The right side shows the last seven spectra right before the break.The -axis shows the reduces VDOS at the same scale for all plots and the -axis is rescaled to the van Hove frequency of the crystalline configuration which is determines as the first low-frequency peak in the crystalline spectrum shown in Figure 4.The inlays show the time-step of each spectrum.The first observation is similar to Figure 4, the overall reduced VDOS is decreasing and becoming blurred with increasing heterogeneity (from top to bottom).The distinct peaks from the first row with heterogeneity 0.60 are decreasing.The peak at lower frequencies eventually diminishes to a level below most of the higher frequencies.The higher frequencies seem to transition into lower frequencies and the maximum frequency decreases.An interesting property of the reduced VDOS is that the less disordered sample seems to have a continuous increase of vibrational modes in certain regions over the whole shear process, whereas the increase in the samples with heterogeneities 1.00 and 1.40 appear more sudden immediately before the first rearrangement.
The transition of the modes to lower frequencies indicates an increase of noninteracting modes in the configuration [29], which correlates with the applied strain in the -direction.For more ordered configurations, this transition proceeds slowly during shear.However, for disordered materials, the mode transitions (in terms of their frequency) happen on a shorter timescale.The reason for this phenomenon can not be described by the spectrum only, but needs further investigation on the picture of the modes themselves, which will be part of a future work.

CONCLUSION AND OUTLOOK
In this paper, we performed AQS simulations on configurations of different heterogeneities.The corresponding Hessian matrices were determines in order to calculate the (reduced) vibrational densities of state.The VDOSs were compared regarding their qualitative structure.The characteristics were an overall decrease in energy for the modes with increasing shear deformation.Furthermore, the modes of more disordered configurations did not show sharp peaks, but flattened and broadened extrema.The stress strain curve for configurations has already been addressed in previous works [30].In this work, we used the stress strain curve to determine elastic constants and to find the first stress drop of the protocol.Several topics open up for further investigation.The most obvious step might be the investigation of the modes themselves.The distinction between different classes of localization and diffusion give rise to look deeper into the mode sequence.Further research may establish a mode transition matrix in order to track the evolution of the modes along the frequency (energy) scale during shear and reveal details on the mode sequence and energy change.Another point may be the investigation of local Hessian matrices and their corresponding local VDOSs, originating from an approach similar to the local yield stress method proposed by Barbot et al. [31].

F
I G U R E 4Reduced VDOS with the -axis scaled to the van Hove singularity (first peak in the crystalline spectrum).The graphs show the reduced VDOS for heterogeneities 0.00 (crystalline), 0.60 and 1.40.

F I G U R E 5
Various VDOS for heterogeneities 0.60, 1.00, and 1.40 during shear.The axes are scaled equally in all plots.

A
C K N O W L E D G M E N T S This research was funded by the Federal Ministry of Education and Research (BMBF) and the Ministry of Culture and Science of the German State of North Rhine-Westphalia (MKW) under the Excellence Strategy of the Federal Government and the Länder.Open access funding enabled and organized by Projekt DEAL.O R C I D Tobias Focks https://orcid.org/0000-0001-6582-3356Franz Bamer https://orcid.org/0000-0002-8587-6591R E F E R E N C E S