# The Yield-Strain and Shear-Band Direction in Amorphous Solids Under General Loading

Ashwin J., Oleg Gendelman, Itamar Procaccia and Carmel Shor Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel
Faculty of Mechanical Engineering, Technion, Haifa 32000, Israel.
July 12, 2019
###### Abstract

It is well known experimentally that well-quenched amorphous solids exhibit a plastic instability in the form of a catastrophic shear localization at a well defined value of the external strain. The instability may develop to a shear-band that in some cases is followed by a fracture. It is also known that the values of the yield-strain (and yield-stress), as well as the direction of the shear band with respect to the principal stress axis, vary considerably with variations in the external loading conditions. In this paper we present a microscopic theory of these phenomena for 2-dimensional athermal amorphous solids that are strained quasi-statically. We present analytic formulae for the yield-strains for different loading conditions, and well as for the angles of the shear bands. We explain that the external loading conditions determine the eigenvalues of the quadrupolar Eshelby inclusions which model the non-affine displacement field. These inclusions model elementary plastic events and determine both the yield-strain and the direction of the shear-band. We show that the angles of the shear bands with respect to the principal stress axis are limited theoretically between and . Available experimental data conform to this prediction.

## I Introduction

The physics of amorphous solids raises some problems that do not exist in perfect crystals. While it is well known that disorder brings about (Anderson) localization of eigenfunctions 58PWA (), it is commonly assumed that such localization is limited to eigenfunctions associated with high energies. Low-frequency eigenfunctions are believed to be spatially extended. This picture fails in amorphous solids that are subjected to an external strain. Here one finds that the lowest energy eigenfunction can become localized. This localization can take the form of a plastic event that changes the local organization of particles through the creation of essential non-affine displacement. Sometimes the plastic instability can exhibit itself as a system spanning event of shear localization along a line in 2-dimensional systems 12DHP () and a plane in 3-dimensional systems 13DGPS (). These events are crucially important in limiting the toughness of technologically important materials like metallic glasses. These instabilities are the subject of this paper.

To be as precise as possible, we will deal in this paper with 2-dimensional systems (although the 3-dimensional extension is available 13DGPS ()) which are athermal (at temperature ) and strained quasi-statically 06ML (). We will consider systems which are good glass formers, containing particles in a volume , interacting via generic interaction potentials that are sufficiently smooth (with at least two continuous derivatives). The total energy can be written then in terms of the positions of these particles, . The Hessian matrix is defined as the second derivative 99ML ()

 Hij≡∂2U(r1,r2⋯,rN)∂ri∂rj . (1)

The Hessian is real and symmetric, and therefore can be diagonalized. Excluding Goldstone modes whose eigenvalues are zero due to continuous symmetries, all the other eigenvalues are real and positive as long as the system is mechanically stable. In equilibrium, without any mechanical loading, the eigenfunctions associated with the large eigenvalues are localized due to Anderson localization, as it was mentioned above. However, all the eigenfunctions associated with the low eigenvalues (including all the excess modes that are typical to amorphous solids) are spatially extended. It had been a major discovery of the last decade that at small values of the external strain there appear “fundamental plastic events” in which the eigenvalues of Hessian matrix hit zero via a saddle node bifurcation, and simultaneously the associated eigenfunctions get localized 12DKP (). This is a new mechanism for localization, different from the well known Anderson mechanism, and it is the basis for the understanding of plastic instabilities in amorphous solids. It had been shown that at low values of the external stains the spatial ranges of these plastic events are system-size independent 10KLP (), involving a small number of particles, of the order of 100. The displacement fields associated with fundamental plastic events are reliably modeled by Eshelby inclusions (see below for details). It had been a more recent discovery that at higher values of the external strains these fundamental plastic events can organize in highly correlated lines in 2-dimensional system or in planes in 3 dimensional systems 12DHP (); 13DGPS (). These highly correlated densities of Eshelby inclusions are the microscopic manifestation of the shear localization.

The key idea of the recent work on the shear localization is that at with quasi-static loading the preferred spatial organization of the density of Eshelby inclusions could be found by minimizing their total energy in the strained system. In Refs. 12DHP (); 13DGPS () the theory was worked out for simple external shear in both 2 and 3 dimensions. It was found that in the case of a pure external shear the shear-localization is realized by organizing the inclusions on a line (plane) in 2 (3) dimensions that are precisely at with respect to the principal stress axis. An additional important result is an analytic prediction for the yield-strain (where shear localization occurs for the first time) in terms of the properties of the Eshelby solution for the fundamental plastic event (a single inclusion).

It is well known however that for other loading conditions, e.g. uniaxial tension or compression, the value of the yield strain as well as the direction of the shear band can vary considerably 11GWBN (), indicating that the pure shear is a special case. It was reported recently in 13AGPS () that the difference can be related to the more general form of the Eshelby inclusion which models the fundamental plastic instability at different loading conditions. The aim of the present paper is to offer the theory in sufficient detail and to work out analytic predictions of the yield strain. The calculations involved are straightforward but sometimes cumbersome, and we make an effort to offer them in an optimally didactic way in the sequel.

In Sect. II we present results for AQS numerical simulations in uniaxial compression and extension to provide us with data on yield strains and directions of shear bands. We find the well known asymmetry in both measures. In Sect. III we begin to develop the general theory for 2-dimensional Eshelby inclusions that is appropriate for the most general loading conditions. In Sect. IV we compute the all-important energy of inclusions, and minimize it to find the preferred organizations and the resulting angle of the shear band. On the basis of this calculation we offer in Sect. V an analytic prediction for the yield strain. The final section VI provides a short summary and a discussion.

The reader should note that throughout this paper we will reserve the Greek letters for tensor subscripts, while the Latin letters will be used for particles and Eshelby inclusions.

## Ii Numerical Simulations

To prepare data for the present analysis we have performed two-dimensional (2D) Molecular Dynamics simulations on a binary system which is known to be a good glass former and has a quasi-crystalline ground state 87WSS (); 88LB (). Each atom in the system is labeled as either “small”(S) or “large”(L) and all the particles interact via Lennard-Jones (LJ) potential. All distances are normalized by , the distance at which the LJ potential between the two species becomes zero and the energy is normalized by which is the interaction energy between the two species. Temperature was measured in units of where is Boltzmann’s constant. For detailed information on the model potential and its properties, we refer the reader to Ref 87WSS (). The number of particles in our simulations is 10000 at a number density with a particle ratio . The mode coupling temperature for this system is known to be close to 0.325. The mass of all particles is and time is normalized to . For the sake of computational efficiency, the interaction potential is smoothly truncated to zero along with its first two derivatives at a cut-off distance . To prepare the glasses, we first start from a well equilibrated liquid at a high temperature of which is supercooled to at relatively fast quenching rate of . Then, we equilibrated these supercooled liquids for times greater than , where is the time taken for the self intermediate scattering function to approach of its initial value. Lastly, following this equilibration, we quenched these supercooled liquids deep into the glassy regime at a temperature of at a reduced quench rate of . Following this quench we took the glass to mechanical equilibrium (nearest energy minima) by a conjugate gradient energy minimization. After that we start loading our glasses under athermal quasi-static conditions. Two different loading protocols were applied, one is uniaxial compression and the other - uniaxial extension. The simulations were followed until the yield strain and slightly above to ascertain the angle of the shear band formed. The shear band is marked on the sample by coloring particles using their values of 98FL (). We use a binary coloring scheme which means that when a given particle has , it is colored black, else it is colored white.

We note the asymmetry in , close to 5.5% for compression and 3.5% for extension, and the asymmetry in angles, 46 and 54 respectively. The theory outlined in this paper provides analytic formulae for both measures, cf Eqs. (74) and (78).

## Iii Theory of the Fundamental Plastic Event

In the past, phenomenological models were applied to this problem 75RR (); 98RO (); 03ZES (); 08ZL (); 09ZL (), but a microscopic approach was lacking. In this section we present the theory of 2-dimensional Eshelby inclusions for the general loading conditions.

### iii.1 Two Dimensional Circular Inclusion

Consider an elastic solid having a volume and surface area [Fig. 3]. The material will be assumed to be homogeneous with an elastic stiffness tensor given by . Let a sub-volume with surface area undergo a uniform permanent (inelastic) deformation, such as a structural phase transformation. The material inside is referred to as an inclusion and the material outside is called the matrix. If we could remove this inclusion from its surrounding material then it would attain a state of a uniform strain and zero stress. Such a stress free strain is referred to as the eigen-strain . The eigen stress is then given by .

In reality, the inclusion is surrounded by the matrix. Therefore, it is not able to reach the state of zero stress. Instead, both the inclusion and the matrix will deform and experience an elastic stress field. The Eshelby’s transformed inclusion problem 57JDE () is to solve the stress, strain and displacement fields both in the inclusion and in the matrix.

We consider a 2D circular inclusion that has been strained into an ellipse using an eigen-strain and which allows for a volume change (). A general expression for such a tensor can be written in terms of a unit eigenvectors () and corresponding eigenvalues () as follows:

 ϵ∗αβ=ζn^nα^nβ+ζk^kα^kβ (2)

Using orthogonality of the eigen directions: , we get:

 ϵ∗αβ =(ζn−ζk)2(2^nα^nβ−δαβ)+(ζn+ζk)2δαβ =ϵ∗,0αβ+ϵ∗,Tαβ (3)

where the traceless part and the nonzero trace part are given as

 ϵ∗,0αβ =(ζn−ζk)2(2^nα^nβ−δαβ) ϵ∗,Tαβ =(ζn+ζk)2δαβ (4)

We also assume that the system is acted upon by a homogeneous strain that acts globally (which in our case also triggers the local transformation of the inclusion). This strained ellipsoidal inclusion feels a traction exerted by the surrounding elastic medium resulting in a constrained strain in the inclusion and also exerts a traction at the inclusion-matrix interface resulting in the originally unstrained surroundings developing a constrained strain field .

The eigen-strain in the inclusion is related to the constrained strain via the fourth order Eshelby tensor :

 ϵcαβ=Sαβγδϵ∗γδ (5)

Now for an inclusion of arbitrary shape the constrained strain , both the stress and the displacement field inside the inclusion are in general functions of space. For ellipsoidal inclusions, however, it was shown by Eshelby 57JDE (); 59JDE () that the Eshelby tensor and the constrained stress and strain fields inside the inclusion become independent of space. We work here with a circular inclusion which is a special case of an ellipse and hence for such an inclusion, the Eshelby tensor reads 2004WCB ()

 Sαβγδ=(λ−μ)4(λ+2μ)δαβδγδ+(λ+3μ)4(λ+2μ)(δαδδβγ+δαγδβδ) (6)

where are Lamé coefficients. Plugging Eq (6) in Eq (5), we get

 ϵcαβ=(λ+μ)(λ+2μ)ϵ∗,Tαβ+(λ+3μ)2(λ+2μ)ϵ∗,0αβ (7)

The total stress, strain and displacement field inside the circular inclusion are then given by

 ϵIαβ = ϵcαβ+ϵ∞αβ σIαβ = σcαβ−σ∗αβ+σ∞αβ≡Cαβγδ(ϵcγδ−ϵ∗γδ+ϵ∞γδ) uIα = ucα+u∞α=(ϵcαβ+ϵ∞αβ)Xβ (8)

where the superscript indicates the inclusion. The eigen stress is related to the eigen strain as

 σ∗αβ =Cαβγδϵ∗γδ =2μϵ∗αβ+λϵ∗ηηδαβ (9)

where we have used the following definition of the fourth order elastic stiffness tensor for an isotropic elastic medium:

 Cαβγδ=λδαβδγδ+μ(δαγδβδ+δαδδβγ) . (10)

The stress in the inclusion can now be written down in terms of the independent variables using equations (9) as

 σIαβ=2μ(ϵcαβ−ϵ∗αβ+ϵ∞αβ)+λ(ϵcηη−ϵ∗ηη+ϵ∞ηη)δαβ (11)

### iii.2 Constrained Fields in the Matrix

In the surrounding elastic matrix, the stress, strain and displacement fields are all explicit functions of space and can be written as

 ϵmαβ(X) = ϵcαβ(X)+ϵ∞αβ σmαβ(X) = σcαβ(X)+σ∞αβ umαβ(X) = ucαβ(X)+u∞αβ (12)

The displacement field in the isotropic elastic medium will satisfy the Lamé-Navier equation (without any body forces) 70LL ()

 (μ+λ)∂2ucγ∂Xα∂Xγ+μ∂2ucα∂Xβ∂Xβ=0 (13)

The constrained fields in the inclusion will supply the boundary conditions for the displacement field in the matrix at the inclusion boundary. Also as , the constrained displacement field will vanish.

All solutions of Eq. (13) also obey the higher order bi-harmonic equation

 ∂4ucα∂Xβ∂Xβ∂Xλ∂Xλ=∇4ucα=0 (14)

Thus our objective is to construct from the radial solutions of the bi-Laplacian equation Eq. (14) derivatives which also satisfy Eq. (13). Note that the bi-Laplacian equation is only a necessary (but not a sufficient) condition for the solutions and Eq. (13) still needs to be satisfied.

### iii.3 Solution of the Lamé-Navier Equation

From the foregoing section, we note that the constrained displacement field due to the Eshelby solution is given as

 ucα =ϵcαβXβ =[(λ+μ)(λ+2μ)ϵ∗,Tαβ+(λ+3μ)2(λ+2μ)ϵ∗,0αβ]Xβ =(ζn+ζk)(λ+μ)2(λ+2μ)Xα+(λ+3μ)2(λ+2μ)ϵ∗,0αβXβ (15)

where we take the expression for from Eq. (7). We will look for linear combinations of the derivatives of the radial solutions of the bi-harmonic equation (14) which are linear in the eigen-strain and go to zero at large distance. In addition the terms must transform as a vector field. Let and be the solutions to the Lamé-Navier equations. We then have:

 (μ+λ)∂2uc,Tγ∂Xα∂Xγ+μ∂2uc,Tα∂Xβ∂Xβ =0 (μ+λ)∂2uc,0γ∂Xα∂Xγ+μ∂2uc,0α∂Xβ∂Xβ =0 (16)

Using the radial solutions of the Lamé-Navier equation we can construct the following combinations which transform as a vector and also go to zero as .

 uc,Tα=Aϵ∗,Tαβ∂lnr∂Xβ+Bϵ∗,Tβγ∂3lnr∂Xα∂Xβ∂Xγ+Cϵ∗,Tβγ∂3(r2lnr)∂Xα∂Xβ∂Xγ (17)

and

 uc,0α=A′ϵ∗,0αβ∂lnr∂Xβ+B′ϵ∗,0βγ∂3lnr∂Xα∂Xβ∂Xγ+C′ϵ∗,0βγ∂3(r2lnr)∂Xα∂Xβ∂Xγ (18)

Using the identities

 ∂2lnr∂Xβ∂Xβ=0;∂2(r2lnr)∂Xβ∂Xβ=4lnr+4 (19)

we see from Eq (17)

 ∂2uc,Tα∂Xβ∂Xβ=4Cϵ∗,Tηλ∂3lnr∂Xα∂Xη∂Xλ (20)

and similarly

 ∂2uc,Tγ∂Xα∂Xγ=(A+4C)ϵ∗,Tηλ∂3lnr∂Xα∂Xη∂Xλ (21)

Plugging Eq (20), (21) into Eq (16), we get

 C=−A(λ+μ)4(λ+2μ) (22)

We can now rewrite Eq (17) as

 uc,Tα=Aϵ∗,Tαβ∂lnr∂Xβ+Bϵ∗,Tβγ∂3lnr∂Xα∂Xβ∂Xγ−A(λ+μ)4(λ+2μ)ϵ∗,Tβγ∂3(r2% lnr)∂Xα∂Xβ∂Xγ (23)

We now compute the following identities

 ∂lnr∂Xβ =Xβr2 ∂3lnr∂Xα∂Xβ∂Xγ =−2r2(Xαδβγ+Xβδαγ+Xγδαβ)+8XαXβXγr6 ∂3(r2lnr)∂Xα∂Xβ∂Xγ =2r2(Xαδβγ+Xβδαγ+Xγδαβ)−4XαXβXγr4 (24)

Using these identities, we can now write down Eq (23) as

 (25)

and similarly

 uc,0α=A′ϵ∗,0αβXβr2−[2B′r4+A′(λ+μ)2(λ+2μ)r2]ϵ∗,0βγ(Xαδβγ+Xβδαγ+Xγδαβ)+[8B′r6+A′(λ+μ)(λ+2μ)r4]ϵ∗,0βγXαXβXγ (26)

Now using Eq. (4), we can rewrite equations (25) and (26) as

 uc,Tα=(ζn+ζk)μAXα2(λ+2μ)r2 (27)
 (28)

The complete solution for the displacement field in the matrix is then given as

 ucα =uc,Tα+uc,0α =(ζn+ζk)μAXα2(λ+2μ)r2+(A′μ(λ+2μ)r2−4B′r4)ϵ∗,0αβXβ+(8B′r6+A′(λ+μ)(λ+2μ)r4)ϵ∗,0βγXαXβXγ (29)

Now at (the inclusion boundary), the form of solution (29) must match with the Eshelby solution (15). This implies,

 A=a2(λ+μ)μ ,A′=a2 ,B′=−a4(λ+μ)8(λ+2μ) (30)

Plugging equation (30) into Eq. (29), we get:

 ucα=(λ+μ)2(λ+2μ)(a2r2)[(ζn+ζk)Xα+(2μλ+μ+a2r2)ϵ∗,0αβXβ+2(1−a2r2)ϵ∗,0βγXαXβXγr2)] (31)

Noting that

 ϵ∗,0αβXβ=(ζn−ζk)2(2^nα(^n⋅X)−Xα) (32)

and

 ϵ∗,0βγXαXβXγ=(ζn−ζk)2(2(^n⋅X)2−r2)Xα (33)

we find that Eq. (31) finally becomes

 uc(X) +2(1−a2r2)(2(^n⋅^r)2−1)X] (34)

where . The Cartesian components of Eq. (34) used for visualizing the displacement field in space are given below:

 ucx(X) +2(1−a2r2)((x2−y2)cos2ϕ+2xysin2ϕr2)x]
 ucy(X) +2(1−a2r2)((x2−y2)cos2ϕ+2xysin2ϕr2)y] (35)

where the unit vector makes an angle with the x-axis. Taking derivatives of the displacement field

 ∂ucα∂Xβ −4(μλ+μ+a2r2)(2(^n⋅^r)^nα−Xαr)Xβr +(2μλ+μ+a2r2)(2^nα^nβ−δαβ)−4(1−2a2r2)(2(^n⋅^r)2−1)XαXβr2 (36)

and

 ∂ucβ∂Xα −4(μλ+μ+a2r2)(2(^n⋅^r)^nβ−Xβr)Xαr +(2μλ+μ+a2r2)(2^nα^nβ−δαβ)−4(1−2a2r2)(2(^n⋅^r)2−1)XαXβr2 (37)

The constrained strain in the matrix can be written as

 ϵcαβ=12(∂ucα∂Xβ+∂ucβ∂Xα) (38)

Using equations (36) and (37), Eq. (38) becomes

 ϵcαβ(X) +(2μλ+μ+a2r2)(2^nα^nβ−δαβ)−4(1−2a2r2)(2(^n⋅^r)2−1)XαXβr2 +4(1−a2r2){(^n⋅^r)(^nβXαr+^nαXβr)−2(^n⋅^r)2XαXβr2}+2(1−a2r2)(2(^n⋅^r)2−1)δαβ] (39)

It is easy to see that the trace of is given as

 ϵcηη(X)=−(ζn−ζk)μ(λ+2μ)(a2r2)(2(^n⋅^r)2−1) (40)

We can now calculate the constrained stress in the elastic matrix due to the deformed Eshelby inclusion. It follows from Hooke’s law:

 σcαβ(X)=2μϵcαβ(X)+λϵcηη(X)δαβ (41)

Plugging Eq. (39) in Eq. (41), we get the final expression for the constrained stress

 σcαβ(X)=(ζn−ζk)μ(λ+μ)2(λ+2μ)(a2r2)[2(ζn+ζk)(ζn−ζk)(δαβ−2XαXβr2)−4(μλ+μ+a2r2){(^n⋅^r)(^nαXβr+^nβXαr)−XαXβr2} +(2μλ+μ+a2r2)(2^nα^nβ−δαβ)−4(1−2a2r2)(2(^n⋅^r)2−1)XαXβr2 (42) +4(1−a2r2){(^n⋅^r)(^nβXαr+^nαXβr)−2(^n⋅^r)2XαXβr2}+2(1−a2r2)(2(^n⋅^r)2−1)δαβ−2λλ+μ(2(^n⋅^r)2−1)δαβ] .

## Iv Energy of N Eshelby Inclusions Embedded in the Matrix

The energy of the Eshelby inclusions embedded in a linear-elastic matrix is given by the following expression

 E=12N∑i=1∫Vi0σiαβϵiαβdV+12∫V−∑Ni=1Vi0σmαβϵmαβdV (43)

where the superscript denotes the index of the inclusion and denotes the matrix. Eq (43) can be re-written using the definition of the strain , where :

 E = 14N∑i=1∫Vi0σiαβ(uiα,β+uiβ,α)dV (44) + 14∫V−∑Ni=1Vi0σmαβ(umα,β+umβ,α)dV

Using the symmetry of the stress tensor, we obtain

 E=12N∑i=1∫Vi0σiαβuiβ,αdV+12∫V−∑Ni=1Vi0σmαβumβ,αdV (45)

If there are no body forces, we also have the identity

 σαβuβ,α=(σαβuβ),α−σαβ,αuβ=(σαβuβ),α (46)

Thus we can write Eq. (45) as

 E=12N∑i=1∫Vi0(σiαβuiβ),αdV+12∫V−∫Ni=1Vi0(σmαβumβ),αdV . (47)

Using Gauss’s theorem, we convert these volume integrals into area integrals to obtain

 E = 12N∑i=1∫Si0σiαβuiβ^niαdS−N∑i=112∫Si0σmαβumβ^niαdS (48) + 12∫S∞σmαβumβ^n∞αdS ,

where and are unit vectors both pointing outwards from the surfaces bounding the inclusion volume and the matrix boundary respectively. Eq. (48) can be re-written as follows

 E=12∫S∞σmαβu∞βdS+12N∑i=1∫Si0(σiαβuiβ−σmαβumβ)^niαdS =12σ∞αβϵ∞βγ∫S∞Xγ^n∞αdS+12N∑i=1∫Si0(σiαβuiβ−σmαβumβ)^niαdS =12σ∞αβϵ∞βαV+12N∑i=1∫Si0(σiαβuiβ−σmαβumβ)^niαdS (49)

Thus we can re-write Eq. (12) as

 ϵmαβ(X) = N∑i=1ϵc,iαβ(X)+ϵ∞αβ σmαβ(X) = N∑i=1σc,iαβ(X)+σ∞αβ umαβ(X) = N∑i=1uc,iαβ(X)+u∞αβ (50)

where denotes the constrained strain at in the matrix due to the Eshelby inclusion. We also have for locations inside the inclusions

 ϵiαβ(X) = ∑j≠iϵc,jαβ(X)+ϵc,iαβ(X)−ϵ∗,iαβ+ϵ∞αβ σiαβ(X) = ∑j≠iσc,jαβ(X)+σc,iαβ(X)−σ∗,iαβ+σ∞αβ uiα(X) = ∑j≠iuc,jα(X)+uc,iα(X)−ϵ∗,iαβXβ+u∞α (51)

where is the eigen-strain of the Eshelby inclusion and so on. Note that in the expression for the strain in the inclusion given by Eq, (51), we have subtracted the contribution of eigen-strain from the constrained strain leaving only the elastic contribution to calculate correctly the elastic contribution to the energy. Using these expressions, the elastic energy of the system can be written from Eq. (49) as follows:

 E=12σ∞αβϵ∞βαV+12N∑i=1∫Si0(σiαβuiβ−σmαβumβ)^niαdS (52)

Since the traction force has to be continuous at the inclusion boundary (Newton’s III