Interaction of defects in disordered materials – Insights from a nonlinear decomposition of the nonaffine displacement field

In this paper, we link the unique mechanical properties to structural aspects using a splitting method in order to disentangle nonlinear interactions of Stone–Wales defects as fundamental perturbations in a crystalline environment. We investigate the spatial features of the displacement field achieved by this splitting method and couple it to macroscopic material properties. The nonlinear interaction broadens the understanding of material failure at the onset of disorder and may also play a role in the emergence of the Boson peak.


INTRODUCTION
In this work, we want to deepen the knowledge of the interaction of Stone-Wales (SW) defects.Such defects are induced by the rotation of bonds, as shown in Figure 1, and constitute one of the basic defects that lead to disorder present in non-crystalline networks [1].
According to this figure, one bond switch is realized by a rotation of a compound of seven atoms around its center.Iteratively inducing random bond switches creates disorder, which plays a major role in the breaking process of 2D materials, such as graphene, silica bilayers, or other hexagonal model materials [2][3][4][5][6].We consider the SW defects as the main source of disorder in the model material.Under this assumption, we generate certain defective hexagonal configurations that fulfill the conditions of Zachariasen's theory for network glasses [3,[7][8][9][10][11] from a crystalline initial state.We perform uniaxial tension, under athermal quasistatic conditions with the created structures.In athermal quasistatic simulation (AQS), the systems are driven slowly, such that the atomic responses are considered as instantaneous.In doing so, high frequency modes are suppressed and Brownian motion is neglected.This allows for an isolated treatment of the atomic dynamics without any disturbances generated by noise.For our simulations, we generate selected configurations with two defects, which are separated in loading direction varying their distance.All simulations were performed using periodic boundary conditions.We investigate the nonaffine displacement fields of these selected configurations.The nonaffine displacement field is known to carry crucial information about material stability, elasticity, and fracture behavior in noncentrosymmetric lattices [12][13][14].We apply a method to further split the nonaffine displacement field into linear and nonlinear components in order to detect possible origins of material stability.The residual of the nonlinear displacement field splitting helps us to understand the dynamics caused by the interaction of defects and their effect on stability of certain regions in the simulated configuration.Herein, we will focus on two types of SW defect arrangements: one configuration, where the defects are aligned with the loading direction and another configuration where the two defects are rotated relative to the loading direction.

Model potential
The defective hexagonal structures investigated in this work are model network materials generated from a lattice originating from a honeycomb structure with a total number of 4410 atoms.The interaction is defined by a Yukawa-type potential, written as: with the fitting parameters   ,   , and  taken from [15].The potential energy landscape (PEL) is then the sum of all interatomic interactions within a predefined cutoff radius  cutoff , written as In this paper, we chose a value of 10.0 Å for  cutoff .The atoms are placed in a honeycomb-shaped hexagonal lattice in a way such that the system is located at its global minimum (crystalline state).In order to create a defective structure, we use a bond switching algorithm [16][17][18].A corresponding dual lattice is defined based on the initial configuration [19].In the dual lattice, one of the bonds is reconnected and the structure is relaxed.Afterward, the dual lattice is transformed back to the physical structure.Using this procedure guarantees that after each switch, the number of connected atoms remains the same.Additionally, the structure is relaxed with respect to its potential energy.
The applied loading is a uniaxial tension protocol [13,20].Each atom in the structure is affinely set to a new position with respect to the pulling axis.In doing so, the neighborhood of every atom will change and the net force on each atom might not add up to zero, as it is the case in Bravais lattices [21,22].In a second step, the PEL of the structure is minimized.The closest local minimum is computed using the conjugate gradient method [23].

Displacement field
The displacement field of Bravais lattices under athermal quasistatic loading is correctly described by the Born approximation [24].However, this does not hold for non-bravais lattices or structures with a defective or disordered neighbourhood.Due to the nonzero forces resulting from a displacement, these forces need to be compensated by a non-affine displacement.The equation of motion was derived by Maloney and Lemaître [24,25].We evaluate the force on each atom  as: where  denotes the potential energy which depends on the atomic position   and the deformation state which we summarize by the box elongation Δ = 0.05 Å.In Equation ( 3),   denotes the position of particle .Since this is a quasistatic scheme, the force does not change during loading and remains zero: The circle • ⋅ refers to the particle's position in the reference state.The term following the equal sign in Equation ( 4) represents a virtual work   , that emerges when taking the atom out of its equilibrium position.This virtual force remains in equilibrium with the corresponding force caused by the change of position of atom  with respect to the Hessian of the configurations.When we remove the two translational modes (Goldstone modes [7]) from the Hessian, the Hessian becomes invertible and Equation (4) can be rewritten as: The displacement caused by the update of the atoms with respect to the loading protocol is defined as the affine displacement field.The trajectory of the minimization in the PEL can be identified as the nonaffine displacement field [14].The decomposition of the atomic motion into the affine and the nonaffine part is linear and written as: We further split the nonaffine displacement field for two separate defects, that is, their linear and nonlinear superposition using the proposed method from [26].Therefore, we take the displacement field representing the trajectory of the atoms during the bond switch algorithm and the subsequent minimization as a field variable: The atom positions refer to the configurations after a single bond switch, including the minimization and the configuration of the underlying honeycomb lattice.The nonlinear interaction is then the residual of the linearly superimposed nonaffine displacement fields of the configurations with fewer bond switches that we call less defective: F I G U R E 2 Relative tensile strength of selected configurations loaded under uniaxial tension in the armchair (AC) direction.
In this equation, the variable In order to underline that this quantity is not directly a result from a simulation, but in this case calculated as the nonaffine displacement field of defective structure minus the hexagonal structure, these quantities are named

Configurations
It was shown that the orientation and the distance of defects in the lattice and relative to the loading direction are crucial for the dynamics of the atoms during tensile loading [27].Consequently, we call these SW defects as type SW1, which have one of their axes of symmetry aligned with the loading direction, and type SW2, defects which are rotated by 30 • relative to the loading direction.Furthermore, we call the directions normal to the atom bonds zigzag (ZZ) direction, which is parallel to the -direction as shown in Figure 1.Furthermore, we call the direction parallel to a vertical bond armchair (AC) direction, which is the -direction as shown in Figure 1.This denotation allows for a definition of the configurations with respect to the type of their defects and loading direction.The denomination of the configurations takes into account the types of the first and the second defect as well as the loading direction.In the following section, we will find an explanation of why the relative tensile strength of some configurations with two SW defects is higher than their single defect counterpart [28], as shown in Figure 2. In doing so, we build our theory upon the nonaffine, nonlinear interaction displacement field defined by Equation ( 8).

RESULTS
Figure 2 shows the results for the relative tensile strength of four selected configurations, replotted from Focks et al. [26], corresponding to Fig. 14 strengths of the SW2-SW2-AC configuration as a function of the distance between the two defects.As mentioned, the relative tensile strength for this configuration turns out to be higher than the tensile strength of a single SW2-type defect.
For larger distances between the defects, the tensile strength seems to approach the line of the single SW2-type defect form above, indicating a decoupling phenomenon.The square marked line shows the tensile strength of the SW1-SW1-AC configuration that starts considerably lower than each of the single defect configurations.The line approaches the SW1-type tensile strength for larger distances, indicating the decrease of weakening interaction between the two defects.
In Figure 3, the result of the separation of the displacement fields for strains (shown in percent in the inlays) shortly before the break can be seen.
In Figure 3A, we present the nonaffine interaction displacement field of the SW2-SW2-AC configuration with a distance of four rings in Figure 3A, the nonaffine simulated field of the same configuration in Figure 3B, the nonaffine interaction field of the SW2-SW2-AC configurations with a distance of 10 rings between the defects in Figure 3C, and the nonaffine interaction field of the SW1-SW1-AC configuration with a distance of four rings between the defects in Figure 3D. Figure 3A shows the interaction displacement field of the SW2-SW2-AC configuration at different strains.It can be seen that the overall interaction is increasing close to the break point.The movement of the atoms is larger along the preferred axes.These movements are along the ZZ direction away from the center point between the defects.The arrows themselves point roughly toward the defects, creating a stabilizing environment for the weak spots, that is, the defects.Comparing the interaction field with the nonaffine displacement field of the simulation, one notices that atoms far away from the defects and atoms that are not located on the preferred axes undergo much larger displacements relative to the preferred axes.Furthermore, the increase in displacement before the break is smaller.However, the most intriguing phenomenon is that the arrows point in different directions in this plot, that is, away from the defects almost aligned with the -axis.Figure 3C shows the interaction displacement field for the SW2-SW2-AC configuration for a larger distance between the defects.In this case, the interaction spreads throughout the material in a quadrupolar pattern as strain increases.The displacement points in the direction of the center of each defect, similar to Figure 3A. Figure 3D shows the interaction field of the SW1-SW1-AC configuration for the same distance as Figure 3B,C.Here, the axes of symmetry of the defects are aligned with the loading direction.It is noticeable, that in the first figures, larger arrows are abundant in the y-direction of the defects.The arrows point in the direction of the defects, exerting a stabilizing effect on the weak spots.During further strain, these stabilizing branches vanish, resulting in an interaction mainly perpendicular to the loading direction.As the interaction cannot compensate for the destabilizing loading, this configuration is, generally, weaker than the SW2-SW2-AC configuration and even weaker than a single SW1-or SW2-type defect, as shown in Figure 2.

CONCLUSION
In this work, we performed a separation of displacements into affine and nonaffine parts and a subsequent nonlinear decomposition of SW defect interaction to describe material stability based on spatial features in the nonaffine displacement fields.We were able to find evidence for enhancement of tensile strength in certain configurations and the reduction of tensile strength in other configurations.Additionally, we traced the stability relations between configurations hosting one and two defects back to displacement modifications caused by the interaction of the defects.This was achieved by comparing the displacement field of the nonlinear interaction of the two defects with the nonaffine displacement field obtained by the AQS of uniaxial tension applied to the configuration at different time steps close to material failure.The decomposition of nonaffine displacement fields into linear and nonlinear parts appears to be a powerful tool for studying the dynamics of defect interaction.For future research, this concept can be applied to more complex structures to gain deeper insight into the origin and the dynamics of material failure of disordered structures.
displacement field of the configuration with two defects,  •   () is the nonaffine displacement field of the underlying hexagonal structure, and  •   is the reduced nonaffine displacement field of the configuration with only one defect: Figure2shows the results for the relative tensile strength of four selected configurations, replotted from Focks et al.[26], corresponding to Fig.14(a), (c).The dash dotted and the dotted line show the relative tensile strength of the SW1-type and the SW2-type SW defect under uniaxial tension in AC direction, respectively.The triangular marks show the tensile 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.Furthermore, the authors greatly thank M. Falk for helpful comments and fruitful discussions.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-6591