Amorphous materials of homogeneous structures usually suffer from nonuniform deformation under shear, which can develop into shear localization and eventually destructive shear band. One approach to tackle this issue is to introduce an inhomogeneous structure containing more than one phase, which can reduce the local nonuniform shear deformation and hinder its percolation throughout the system. Using thermostated molecular dynamics (MD) simulations, we compare the deformation behavior between a homogeneous amorphous mixture of bidisperse disc particles, interacting via an *n* − 6 Lennard-Jones potential of tunable softness, with an inhomogeneous one containing an evenly-distributed ordered phase. We change the population ratio of large to small particles to create a homogeneous or an inhomogeneous mixture, where the softness of a chosen phase can be manually adjusted by specifying *n* of the interparticle potential. Results of applying extensive quasistatic shear on the prepared mixtures reveal that the inhomogeneous amorphous mixture containing a soft ordered phase overall deforms more uniformly than the homogeneous one, which indicates that both the structure inhomogeneity and the inter-phase softness variance play important roles in enhancing the uniformity of the plastic deformation under shear.

## I. INTRODUCTION

Homogeneous amorphous materials such as bulk metallic glasses are known for exhibiting superior mechanical properties than their crystalline siblings. However, a significant disadvantage of amorphous materials is their low ductility, due to nonuniform shear deformation, which causes premature and unpredictable failure and greatly limits their industrial applications.^{1,2} Recently, introducing an ordered phase into an amorphous material to make the engineered structure inhomogeneous have been shown to improve their mechanical properties, for example, embedding an isolated and soft crystal phase containing dendrites in a bulk metallic glass matrix to enhance tensile ductility,^{3–8} and distributing polycrystalline metallic alloys within amorphous shells to raise mechanical strength and restrict shear localization.^{9}

In this study, we propose a mesoscale model of mixing 2D bidisperse disc particles to study the shear deformation behavior of an amorphous configuration containing an ordered phase using molecular dynamics (MD) simulations. In our model of 2D amorphous materials, we prepare a homogeneous or an inhomogeneous amorphous/ordered configuration by mixing 50 − 50 or 90 − 10 small and large 12 − 6 Lennard-Jones (LJ) disc particles, which allows the system to form an amorphous or a partially amorphous structure under thermal equilibrium.^{10,11} We set the interparticle interaction to be Lennard-Jones for the potential has long been used to study amorphous materials including metallic glasses.^{12–14} We then identify particles in the amorphous phase or the ordered phase based on the values of their disorder parameter. We alter the strength of a chosen phase by assigning the particles belonging to it proper interparticle softness, using an *n* − 6 Lennard-Jones (LJ) potential.^{15} We choose *n* = 8 to make particles in the ordered phase softer than those in the amorphous phase. Finally, we apply quasistatic shear on the prepared configurations and calculate the uniformity of their deformation under extensive shear.

We investigated the difference in shear deformation between a homogeneous and an inhomogeneous configurations and demonstrate that an inhomogeneous configuration favors more uniform shear deformation due to the spatial and strength heterogeneity between the two phases. Specifically, our model shows that introducing an ordered phase into an amorphous phase helps the inhomogeneous system deform more uniformly under shear than a homogeneous single-phased system when the applied shear strain is small. For large shear strain, the difference in softness between the two phases also plays an important role to further enhance the uniform plastic deformation. Our MD simulation results offer evidence supporting the improved ductility in mesoscale reported in experiments of amorphous materials.^{3–9}

## II. NUMERICAL SIMULATION METHOD

Our MD method includes two parts: 1) generating a homogeneous or inhomogeneous initial amorphous configuration and 2) applying quasistatic shear on it. To create an initial configuration, we first pick up an equilibrium configuration in a liquid state where particles interacting via the 12 − 6 LJ potential and cool it down to a solid temperature. Here temperature is defined as the total kinetic energy divided by degrees of freedom of the system. Then we assign particles different softness according to their disorder parameter and prepare boundary particles. Finally, we relax the configuration sandwiched by the boundary particles at the same solid temperature. To test the deformation uniformity of the prepared initial configuration, we apply quasistatic shear by moving the boundary particles stepwise followed by thermostated relaxation and calculate the average deviation of uniform deformation of the sheared configuration. For each system setup, we use at least ten independent initial conditions to obtain the averaged results. The details are given below.

### A. Preparation of an initial configuration

#### 1. System geometry

Our system is a mixture of total *N* = *N*_{s} + *N*_{l} circular particles interacting via the finite-range, pairwise-additive LJ potential. Specifically, it contains *N*_{s} small particles of diameter *d*_{s} and *N*_{l} large particles of diameter *d*_{l}, with the diameter ratio *r* = 1.4 to avoid artificial crystallization in 2D. We use *N* = 1000 throughout this study. The masses of the small particles *m*_{s} and the large particles *m*_{l} are identical. The system occupies a square simulation box of size *L* on the *xy*-plane, where *x* is the horizontal axis and *y* is the vertical one. For a given particle number *N* and particle diameters *d*_{s} and *d*_{l}, the box size is determined from the condition that the configuration has a fixed area packing fraction *ϕ* = 0.793 or soft-core packing fraction of LJ potential *ϕ*_{s} = 1.0. where *ϕ* and *ϕ*_{s} are defined as

The value 2^{1/6}*d*_{s} or 2^{1/6}*d*_{l} is where the 12 − 6 LJ potential reaches its minimal, detailed below.

#### 2. Interparticle Lennard-Jones potential

We choose the finite-range, pairwise additive LJ potential to build our amorphous model for it has been widely used to study amorphous systems.^{12–14} To create a softness difference between the amorphous and ordered phases in an inhomogeneous amorphous configuration, we choose a tunable *n* − 6 LJ potential^{15}

where *ϵ* is the characteristic energy scale, $\lambda =32[2(n/6)/(n\u22126)]$, *α* = *n*/[2(*n* − 6)], *r*_{ij} is the separation between particles *i* and *j*, $dij=12(di+dj)$ is their average diameter, $cij=4\u03f5[\lambda (dij/rcut)n\u2212\alpha (dij/rcut)6]$ is a constant that guarantees $VLJn\u22126(rij)\u21920$ as *r*_{ij} → *r*_{cut}, and Θ(*x*) is the Heaviside step function. We use *r*_{cut} = 3.2*d*_{s}, where *d*_{s} is the diameter of small particles and *n* = 8 for the soft $VLJ8\u22126$ LJ potential or *n* = 12 for the well-known stiff $VLJ12\u22126$ LJ potential.

#### 3. Tuning inhomogeneity and softness

The whole process of preparing a homogeneous or inhomogeneous initial configuration for our amorphous model can be summarized as follows. We (a) create a randomly packed initial configuration of circular particles without interparticle overlap; (b) equilibrate the system at a liquid temperature *T*_{l} to relax it; (c) pick up a relaxed liquid configuration and equilibrate it again at a lower solid temperature *T*_{s}; (d) attach boundary particles to the relaxed solid configuration and reassign the interparticle interactions of all particles based on the value of their disorder parameters; (e) equilibrate the sandwiched solid configuration one last time at *T* = *T*_{s}.

The MD simulations in this study use the diameter *d*_{s} and mass *m*_{s} of the small particles and the interparticle potential amplitude *ϵ* as the reference length, mass, and energy scales. Mass *m*_{l} of the large particles is the same as *m*_{s}. As a result, the unit of time *t* is $dsms/\u03f5$. We control the temperature *T* of this system using the Nóse-Hoover thermostat with a thermostat moment of inertia *Q*.^{16,17} Here temperature *T* is defined as the total kinetic energy $\u2211i=1Nmivi2/2$ divided by degrees of freedom of the system and measured in units of *ϵ*/*k*_{B}, where *m*_{i} and *v*_{i} are the mass and velocity of particle *i*, and *k*_{B} is the Boltzmann’s constant. Degrees of freedom of the system is 2*N* − 2 or 2*N* − 1 for periodic boundary conditions in both *x* and *y* directions or only in the *x* direction. To closely compare our results with those in the literature,^{10,11,18,19} we set *T* = *T*_{l} = 2.0 and *T* = *T*_{s} = 0.2 to study the system in a liquid state or a solid state, respectively. The unit of *Q* is $msds2$. The Newtonian equations of motion of a *N* particle system excluding boundary particles if any are integrated using the velocity Verlet algorithm.^{20}

We perform a thermostated relaxation of the system excluding boundary particles after generating a random non-overlapped initial configuration and whenever making changes in the simulation temperature, introducing boundary particles or reassigning properties of particles for the pairwise *n* − 6 LJ interparticle interactions. The relaxation of a random initial configuration at *T* = *T*_{l} also erases any beginning memory from the system. We terminate the relaxation process when the temperature fluctuation decays to within ± 5% of the required temperature.

It is known that by mixing bidisperse 2D particles with the population ratio *c* = *N*_{l}/*N* increasing from 0 to 0.5, we can create structure from a single-crystal, a partially-amorphous (polycrystal) structure to an amorphous structure.^{10,11} To create a inhomogeneous amorphous configuration containing an ordered phase, we choose a partially-amorphous structure using *c* = 0.1 and *N* = 1000.

After generating a random initial configuration (see the Appendix for the implementing details), we equilibrate it at *T* = *T*_{l}. We then bring the temperature down to *T*_{s} within a time interval Δ*t* = 10^{5} to equilibrium the system again with periodic boundary conditions in both *x* and *y* directions and *Q* = 40.0. Using the relaxed configuration at *T* = *T*_{s}, we calculate the disorder parameter *D*_{j} of particle *j* defined as^{10,11,21–24}

where *β*_{j} is a local crystalline angle introduced by $e6i\beta j=\u2211k=1Nbe6i\theta jk/\u2211k=1Nbe6i\theta jk$, *i* is the imaginary unit $\u22121$ and *θ*_{jk} is the angle between the separation vector $\u21c0rj\u2212\u21c0rk$ of particle *j* and its bonded neighboring particle *k* and the horizontal *x* axis. Particle *k* is considered bonded to particle *j* as long as $\u21c0rj\u2212\u21c0rk\u22641.5djk$.^{11} The summation runs over all *N*_{b} bonded neighboring particles *k* of particle *j*. The value of *D*_{j} is zero if particle *j* and its bonded neighbors form an ordered perfect hexagonal structure and increases if the structure becomes more disordered.

Fig. 1(a) shows an equilibrated fully-amorphous configuration of *c* = 0.5 at *T* = *T*_{s}, where particles are colored from blue to red with increasing *log*_{10}(*D*_{i}). Fig. 1(b) shows a weak positive correlation between *D*_{i} of particle *i* and the average $D\xafibonded$ of the same quantity of its bonded neighbors, obtained by averaging ten relaxed initial conditions. We can see that particle *i* can have ordered or disordered neighbors regardless of its *D*_{i}. The figure shows two curves of $D\xafibonded$, for small or large particles, respectively. If we plot the probability density function *P*(*D*_{i}) of particle *i* for small and large particles separately, as shown in the inset of Fig. 1(b), we can see quantitatively that almost all particles are very disordered with *D*_{i} > 0.5, showing the configuration is indeed amorphous.

Similarly, Fig. 2 shows an equilibrated partially-amorphous configuration of *c* = 0.1 at *T* = *T*_{s}. We can observe in Fig. 2(a) that a majority of large particles are disordered and form the amorphous phase separating ordered islands made of small particles. The most ordered small particles form the cores of the ordered islands surrounded by less ordered particles. The degree of disorder of particles increases monotonically as a function of the distance measured from the ordered cores. Contrary to Fig. 1(b), Fig. 2(b) shows a strong positive correlation between *D*_{i} of particle *i* and its $D\xafibonded$, also using ten relaxed initial conditions. We can clearly see that very ordered particles tend to have ordered neighbors, but highly disordered particles tend to have equally disordered neighbors, which gives the two curves positive slopes close to unity. Moreover, in the inset of Fig. 2(b), we can see quantitatively that almost all large particles are very disordered with *D*_{i} > 0.5, while only about half small particles are in a similar disordered status.

The insets of Fig. 1 and Fig. 2 disclose the basic principle of how to build an homogeneous or inhomogeneous amorphous structure out of a mixture of small and large particles in this study: in a fully-amorphous configuration, all most all particles have their *D*_{i} > *D*_{t}, where *D*_{t} = 0.5 is a threshold disorder value used in this study. On the other hand, in a partially-amorphous configuration, we can use *D*_{i} < *D*_{t} to identify ordered particles, which form an ordered phase, and the rest particles form an amorphous phase. Of the ten initial configurations used in this study, the average number of ordered particles is 0.4406*N*. The amorphous phases in homogeneous and inhomogeneous configurations are similar in terms of their disorder parameter distributions.

For an inhomogeneous amorphous configuration, we assign the soft interparticle $VLJ8\u22126$ potential to the interaction between two ordered particles of *D*_{i} ≤ *D*_{t}. We assign the same soft $VLJ8\u22126$ potential to the interaction between an ordered particle of *D*_{i} ≤ *D*_{t} and a disordered particle of *D*_{i} > *D*_{t}, happening at the interface between the ordered and amorphous phases. Moreover, the stiff interparticle $VLJ12\u22126$ potential governs the interaction between two disordered particles of *D*_{i} > *D*_{t}. Finally, the interaction between an ordered or disordered particle and a boundary particle is controlled by the same stiff $VLJ12\u22126$ potential. The boundary particles are image particles created using the periodic boundary condition in the vertical *y* axis. We keep enough image particles so that the top and bottom boundaries have thickness equals *r*_{cut} = 3.2*d*_{s}, and the total number of top or bottom boundary particles is about $3.2dsN$. Boundary particles are not thermostated, and we relax the prepared system sandwiched by boundary particles at *T* = *T*_{s} again with periodic boundary condition in the horizontal *x* direction and *Q* = 0.01 before using it for the quasistatic shear tests. We observe only slight local position variation but no large-scale rearrangement of particles in the amorphous or ordered phase of the relaxed configuration. Two equilibrated snapshots of exemplary inhomogeneous and homogeneous amorphous configurations and the rules for interparticle interactions are shown in Fig. 3, where particles in soft-ordered and stiff-amorphous phases are colored in green and orange, respectively.

### B. Quasistatic shear deformation

We apply quasistatic shear strain on the prepared homogeneous and inhomogeneous amorphous configurations to test their plastic deformation behavior. To do this, at each step of the quasistatic shear, we shift the *x* position of each top and bottom boundary particle by a small amount of 0.01*L* and −0.01*L*, respectively. The system heats up due to the work from an increment of shear strain *γ* of the system by 0.02 in the quasistatic shear deformation, which is about two-third of the shear strain of shifting a single boundary particle by its own size. Using a dimensionless MD time step *dt* = 0.001 and the periodic boundary condition in the horizontal *x* direction, we relax the system at temperature *T*_{s} with the Nóse-Hoover thermostat of *Q* = 0.01 integrated by the velocity Verlet algorithm until the standard deviation of temperature *T*, periodically calculated within a time interval of 2, 500 MD steps, decays to smaller than 0.03, which corresponds to a temperature fluctuation within ±5% of the assigned temperature. Boundary particles are not thermostated and their positions stay fixed during the relaxation process. The thermostated system quickly dissipates the sudden temperature increase after each quasistatic shear deformation. We have checked our simulation results carefully and made sure there is no significant particle rearrangement that can change the configurational structure of the system substantially per quasistatic shear deformation. Fig. 4 demonstrates the decay of temperature fluctuation while applying this process on an exemplary inhomogeneous configuration. We repeat the two steps of shifting boundary particles and relaxation of the system at *T*_{s} until a prescribed value of shear strain *γ* is reached.

## III. ANALYSIS OF THE UNIFORMITY OF DEFORMATION UNDER QUASISTATIC SHEAR

As mentioned in the Introduction, we expect that an inhomogeneous amorphous configuration allows the system to shear more uniformly than a homogeneous one for the latter lack of the inhomogeneity of two phases with different softness to deter the emergence of shear localization. To measure the uniformity of shear deformation quantitatively, we calculate the deviation of uniform shear deformation *L*_{σ} that measures how far a sheared configuration at a given strain *γ* is away from its completely uniformly deformed counterpart, defined as

where *x*_{i} and *y*_{i} are calculated positions of the idealized uniform deformation profile, evenly divided into *m* equal-sized stripes, Δ is the width of each stripe, and $x\xafi$ is the averaged *x* positions of particles whose *y* positions located within $yi\u2212\Delta 2$ and $yi+\Delta 2$ of stripe *i*. To make sure that the choice of *m* has no influence on our conclusion, we have tried *m* = 5, 10, 20, 30, 40 and 50 and found that the values of *L*_{σ} stay unchanged within negligible fluctuations. We therefore use *m* = 10 throughout our analysis and ten trials with different initial conditions to calculate an averaged $L\xaf\sigma $ as a function of *γ*, as shown in Fig. 5(a).

Our amorphous model with tunable inhomogeneity allows us to compare the effects of introducing an inhomogeneous configuration and creating a softness difference between the two phases in the inhomogeneous configuration on the uniformity of shear deformation independently. We proceed with our investigation by two stages. First, to test the effect of introducing the inhomogeneous configuration, we compare the shear deformation behavior between an amorphous homogeneous configuration and a partially-amorphous inhomogeneous configuration, where the interparticle interactions are $VLJ12\u22126$ LJ potential in both cases. The results of the comparison is shown in Fig. 5(b). We can see clearly that introducing a configurational inhomogeneity effectively reduces nonuniform deformation when strain *γ* is smaller than about 0.2. When sheared further, the inhomogeneous configuration gradually loses to the homogeneous one.

Second, we incorporate the softness difference between the two phases of the inhomogeneous configuration, where now ordered particles interacting with another ordered or disordered particles via the soft $VLJ8\u22126$ LJ potential, and compare it with the homogeneous one as in the first stage. The results are shown in Fig. 5(c). Strikingly, the tuned inhomogeneous configuration deforms even more uniformly (smaller averaged *L*_{σ}) than the homogeneous one until the shear strain *γ* reaches about 0.5.

Due to the large error bars in Fig. 5(b) and (c), we further verify the results using five times more initial configurations to verify the reliability of the observed trends. We confirm the trends stay the same and therefore saturate the data, as shown in Fig. 6.

To get a better insight into the deformation mechanism under shear of the homogeneous and inhomogeneous configurations, we calculate their per-particle von Mises stress, which is the deviatoric part of the stress tensor, defined as $\sigma iVM=(\sigma ixx)2+(\sigma iyy)2\u2212\sigma ixx\sigma iyy+3(\sigma ixy)2$ in 2D, where $\sigma iab=\u2212[miviavib+\u2211j=1Np12rijafijb]$ for particle *i* influenced by *N*_{p} neighbors within *r*_{cut} via the *n* − 6 Lennard-Jones potential, *f*_{ij} is the force acting on particle *i* from its neighbor *j*, and *a*, *b* ∈ [*x*, *y*].^{25} The results are shown in Fig. 7. We observe that in a inhomogeneous configuration, particles subject to high *σ*^{V M} shear stress mostly distribute within the amorphous phase, and the isolated soft-ordered phase disperses them so that they cannot coordinate to percolate through the whole system easily.

Furthermore, we also calculate the local deviation from affine deformation, $Dmin2$, which identifies local irreversible particle shuffling in unit of $ds2$.^{26} We demonstrate a comparison of $Dmin2$ between exemplary homogeneous and inhomogeneous configurations during the quasistatic shear strain interval *γ* = [0.04, 0.46], as shown in Fig. 8. Both systems have similar probability density distribution $P(log10(Dmin2))$ and a similar amount of particles are subject to irreversible shear deformations equal or greater than the same $max(Dmin2)$ at *γ* = 0.46. In the figure, we can see clearly that particles with large $Dmin2$ concentrate on one side in the homogeneous system, responsible for the higher nonuniform deformation. On the other hand, similar particles are mostly distributed evenly within the amorphous phase, and the intermediate ordered phase hinders their percolation, which enhances the uniform deformation. Our calculations of $\sigma iVM$ and $Dmin2$ offer clear evidence that an inhomogeneous amorphous configuration can effectively deter nonuniform shear deformation than a homogeneous one.

## IV. CONCLUSIONS

In this study, we propose a 2D mixture of bidisperse particles to study the deformation behavior of homogeneous and inhomogeneous amorphous materials. A configuration modeling homogeneous amorphous materials is made of 50 − 50 small and large particles. On the other hand, Our mesoscale model of inhomogeneous amorphous materials distributes large circular particles in a sea of small particles with a 1:9 population ratio. In the inhomogeneous model, about half (≈44%) small particles form an isolated and ordered phase. On the other hand, most large particles and the remaining half small particles together form an amorphous phase filling the space in between the ordered phase. To give the ordered or the amorphous phase proper mechanical softness as proposed in experiments on amorphous materials to improve their ductility, we introduce a tunable *n* − 6 LJ potential. A particle in the soft-ordered phase interacts with another particle via the 8 − 6 LJ potential. Two disordered particles interact with each other via the stiff 12 − 6 LJ potential. We apply quasistatic shear on the homogeneous and inhomogeneous configurations by repeatedly moving the boundary particles with a stepwise shear strain of 0.02, followed by thermostated relaxation.

Our simulation results show that when the applied shear strain *γ* is small, the configurational inhomogeneity between the two phases in an inhomogeneous structure alone works well to improve uniform deformation. However, with increasing the shear deformation further, the difference in softness between the two phases plays an essential role to enhance the uniform shear behavior. In general, the inhomogeneous configurations deform more uniformly than the homogeneous ones. Our simplified model offers clear numerical evidence supporting the experimental results and could open a new quantitative approach to systematically improving the ductility of inhomogeneous amorphous materials including an ordered phase. For the future work, we will examine a hard ordered phase embedded in a soft disordered matrix and explore the optimal ordered/disordered area ratio that gives the best deformation uniformity, which will give a more complete picture of this study.

## ACKNOWLEDGMENTS

GJG acknowledges financial support from Shizuoka University startup funding. YJW acknowledges financial support from National Natural Science Foundation of China (NSFC No. 11672299). We also thank Corey S. O’Hern for insightful comments.

The authors declare no conflict of interest.

### APPENDIX: GENERATING A RANDOMLY PACKED CONFIGURATION OF PARTICLES

When conducting MD simulation in a liquid state, we start with an initial configuration at *ϕ* = 0.793 without overlap between particles. Practically, it is very difficult to generate such random initial configuration using a completely random process which places particles one by one. To avoid this issue, we first generate a mechanically stable (MS) packing of frictionless particles interacting via the finite-range pairwise-additive repulsive spring potential at area packing fraction *ϕ* = 0.84 (*ϕ*_{s} = 1.06), close to random-close packing density in 2D.

The MS packing of particles is generated using the procedure detailed in Ref. 27. We start with a sparse initial configuration, and the procedure increases the sizes of the particles followed by energy minimization to remove overlap between particles. Periodic boundary conditions are implemented in both *x* and *y* directions. Occasionally, particles have to be shrunk if the energy minimization procedure fails to remove interparticle overlap. We repeat the two steps of particle size perturbation and energy minimization until all particles are force-balanced with their neighbors, and any attempt to increase particle size will result in an increase of the total energy of the system.

We then decrease *ϕ* from about 0.84 (*ϕ*_{s} = 1.06) to 0.793 (*ϕ*_{s} = 1.0) manually and run the thermostated MD simulation at *T* = *T*_{l}. Since the system is in a liquid state, it will soon forget it was from an MS state as long as it is fully relaxed.