Incorporating a soft ordered phase into an amorphous configuration enhances its uniform plastic deformation under shear

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][4][5][6][7][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 * koh.kokketsu@shizuoka.ac.jp, gjjgao@gmail.com 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][13][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 heterogenity 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][4][5][6][7][8][9].
Below we elaborate on the 2D model of amorphous materials and the MD simulation methods in section II, followed by quantitative results of comparing the deformation behavior of homogeneous and inhomogeneous configurations under quasistatic shear in section III. We conclude our study in section IV.

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

System geometry
Our system is a mixture of total N = N s + N l circular particles interacting via the finite-range, pairwiseadditive 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.

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][13][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] V n−6 LJ (r ij ) = 4 λ( where is the characteristic energy scale, λ =  6 ] is a constant that guarantees V n−6 LJ (r ij ) → 0 as r ij → r cut , and Θ(x) is the Heaviside step function. We use r cut = 3.2d s , where d s is the diameter of small particles and n = 8 for the soft V 8−6 LJ LJ potential or n = 12 for the well-known stiff V 12−6 LJ LJ potential.

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 d s m s / . 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 N i=1 m i v 2 i 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 2N − 2 or 2N − 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 m s d 2 s . 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][22][23][24] where β j is a local crystalline angle introduced by k=1 e 6iθ jk , i is the imaginary unit √ −1 and θ jk is the angle between the separation vector r j − r k 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 r j − r k ≤ 1.5d jk [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 averageD bonded i 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 ofD bonded i , for small or large particles, respectively. If we plot the probability density  The inset shows the probability density distribution of Di with the same horizontal axis, and the total area below the two curves is one. The vertical dashed line shows a threshold Dt = 0.5, separating ordered and disordered particles in this study. The data are obtained using ten relaxed configurations, where the results of small (S: purple) and large (L: green) particles are plotted separately.
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 itsD bonded i , 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 (color online) The n − 6 LJ interparticle potential used in the 2D model of amorphous materials, shown in (a). The interaction between a soft ordered (green) particle and another soft ordered or stiff disordered (orange) particle is the soft V 8−6 LJ LJ potential. The interaction between two stiff disordered particles or with a boundary particle (black) is the stiff V 12−6 LJ LJ potential. Exemplary snapshots of an inhomogeneous and a homogeneous configurations are shown in (b1) and (b2), respectively. 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.4406N . 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 V 8−6 LJ potential to the interaction between two ordered particles of D i ≤ D t . We assign the same soft V 8−6 LJ potential to the interaction between an ordered particle of D i ≤ D t and a disor- dered particle of D i > D t , happening at the interface between the ordered and amorphous phases. Moreover, the stiff interparticle V 12−6 LJ 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 V 12−6 LJ 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.2d s , and the total number of top or bottom boundary particles is about 3.2d s √ N . 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.01L and −0.01L, respectively.
The system heats up due to the perturbation from the moved boundary particles. 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. 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 equalsized stripes, ∆ is the width of each stripe, andx i is the averaged x positions of particles whose y positions located within y i − ∆ 2 and y i + ∆ 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 averagedL σ 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 and inhomogeneous configurations without softness inhomogeneity (black). In the inset, particles that are more ordered of an inhomogeneous configuration are semi-transparent. (c) Averaged Lσ of the same homogeneous configurations (red) and the same inhomogeneous configurations except with softness inhomogeneity (blue). In the insets, the stiff and soft phases are colored in orange and green, respectively. For each case, the results are obtained using ten initial configurations. 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 V 8−6 (color online) Averaged Lσ as in Fig. 5 (b) and (c), except the results are obtained using fifty initial configurations for both cases.
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, 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, D 2 min , which identi- (color online) Snapshots of an equilibrated (a1) homogeneous and (b1) inhomogeneous configurations of N = 1000 at T = Ts and γ = 0.02, colored by particle types (soft ordered: green; stiff disordered: orange) defined in the initial conditions. c = 0.1 for the inhomogeneous configuration. Particles are colored in shades from blue to red with increasing value of the von Mises stress σ V M in (a2) and (b2) correspondingly. The color scale is relative to the minimal and maximal von Mises stresses σ V M min and σ V M max in each system.
fies local irreversible particle shuffling in unit of d 2 s [26]. We show a comparison of D 2 min between exemplary homogeneous and inhomogeneous configurations during the quasistatic shear strain interval γ = [0.04, 0.46] in Fig. 8. Both systems have similar probability density distribution P (log 10 (D 2 min )) and a similar amount of particles are subject to irreversible shear deformations equal or greater than the same max(D 2 min ) at γ = 0.46. We can see clearly that particles with large D 2 min 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 σ V M i and D 2 min 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.

V. 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 [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.

VI. 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.