Shear Banding in a Contact Problem between Metallic Glasses

: The case of a frictionless contact between a spherical body and a ﬂat metallic glass is studied using a mesoscopic description of plasticity combined with a semi-analytical description of the elastic deformation in a contact geometry (code ISAAC). Plasticity is described by irreversible strain rearrangements in the maximum deviatoric strain direction, above some random strain threshold. In the absence of adhesion or friction, the plastic deformation is initiated below the surface. To represent the singularities due to adhesion, initial rearrangements are forced at the boundary of the contact. Then, the structural disorder is introduced in two different levels: either in the local strain thresholds for plasticity or in the residual plastic strains. It is shown that the spatial organization of plastic rearrangements is not universal, but it is very dependent on the choice of disorder and external loading conditions. Spatial curved shear bands may appear below the contact but only for a very speciﬁc set of parameters, especially those characterizing the random thresholds compared to externally induced strain gradients.


Introduction
Bulk metallic glasses are made of a disordered assembly of metallic atoms. Thus, they exhibit an atomic structure that is very close to the liquid one, but below the glass transition temperature, they are trapped into a metastable mechanical equilibrium identified by the absence of significant viscous flow at rest (viscosity larger than 10 13 poise). This similarity to the liquid structure and the absence of both long-range order and characteristic length, despite the interatomic distance, is characteristic of amorphous materials. The original microstructure of amorphous materials bestows them very specific properties, such as a very smooth and moldable surface, low rigidity, but very high hardness, and a strong capacity to withstand very high loads before failure [1]. Indeed, the atomistic disorder precludes the formation of extended linear defects such as dislocations and thus forces to maintain irreversible deformation restricted to point-like defects. This makes them highly interesting for their ability to withstand heavy loads while maintaining their elasticity [2]. However, these properties are very dependent on the temperature: they can change drastically at the glass transition temperature, which is between 100 and 1000 K depending on the composition [3]. At finite temperatures [4] or under high-frequency acoustic waves [5], some metallic glasses can even exhibit superplasticity, with up to 160% irreversible plastic flow strain without breaking, at centimeter scales. Moreover, at low temperature or strain rate, irreversible (plastic) deformation is observed at small scale, and it can organize along macroscopic shear bands, which break the apparent homogeneity of the deformation and are the precursors of catastrophic failures [6][7][8]. Thus, it is important to understand the origin of the plastic deformation (irreversible non-viscous deformation) and its spatial organization.
The microscopic origin of plastic deformation in metallic glasses was initially related to the emergence and diffusion of local voids described successfully within the framework 2 of 16 of the free volume theory [9,10] and then to isochoric shear transformation zones analogous to dislocation loops [11]. Recently, it was proved that the constitutive laws of plastic flow results from the accumulation of plastic deformations described as Eshelby inclusions [12]. These Eshelby inclusions result from a loss of material stability [13][14][15]. Thus, they unfold in the direction of instability. The latter depends on the glass composition. In metallic glasses, it is mainly due to the local shear. It is generally given by the direction of the maximum deviatoric strain [16].
The occurrence of shear bands may have two origins: either the unfolding of the initial instability along a plane (elementary shear bands [17][18][19]), or the progressive accumulation of damage encouraged by the self-sustainability of the shear bands [20][21][22]. The accumulation of plasticity without self-sustainability can also give rise to the sparse distribution of shear bands, allowing the occurrence of a yielding transition appearing as a universal critical phenomenon [23,24] but without permanent shear banding. Finally, the nucleation of the shear bands may be intrinsic (with a volumetric origin) or extrinsic (from surface defects at the free surfaces, for example) [20], the later favoring permanent shear banding. Thus, under the same loading conditions, the presence of interfaces and specific boundary conditions may induce shear bands that would be absent elsewhere [25,26].
In this paper, we focus on the nucleation of shear bands under a nanoindenter. It was shown that under such a loading condition, glasses with low free volume, such as metallic glasses, and contrary to pure silica for example, have a tendency to develop shear bands [7,8]. In Zr-based metallic glasses, Su and Anand [27] have shown a nice and regular pattern of curved shear bands that was interpreted with the help of a specific constitutive equation. This continuum constitutive law was supposed to relate the local shearing rate to the local shear stresses, weighted by the local pressure, and also to give rise to a kinetic increase of the free volume and consequently to post yield-strain softening [27]. The crossing between the shear bands was characterized by an angle of 84 • . Budrikis et al. [28] have observed similar features in a mesoscopic modeling of plasticity under a nanoindenter. They related them to a universal critical phenomenon for shear bands nucleation at the yield strain for macroscopic plastic flow. This critical phenomenon would be similar to what happens under more simple loadings such as pure shear or hydrostatic compression [23,24]. However, Shi and Falk [29] have shown using molecular dynamics simulations that this pattern depends indeed on the processing history of the glass (cooling rate) and also on the strain rate. Moreover, Amon and Crassous [30] have observed in granular materials a progressive transition from initially curved shear band patterns to straight bands across the sample. Thus, it is not clear whether the shape of such a shear band is universal or not, and how it is nucleated. To understand the necessary ingredients for curved shear band nucleation, we have run a mesoscopic model of local plasticity based on microscopic known processes and a realistic description of long-range elastic couplings.
In Section 2, we will describe the model used and the parameters involved. In Section 3, we will discuss the results, more specifically the role of disorder in the plastic transformation compared to the role of disorder in the local plastic threshold. Section 4 will discuss the conditions for getting shear bands in a contact problem. Section 5 summarizes the results.

Materials and Methods
The material studied is a typical bulk metallic glass. The material is not studied at the atomic scale, nor at a continuous scale, but the numerical model (combining the software ISAAC with local fortran programs, following the steps described in Figure 1 (right) used in this article is inspired from a set of mesoscopic models proposed in the last twenty years to model plasticity in disordered materials at intermediate length-scales (nm-µm) [31]. These models keep the necessary ingredients for local plasticity in solids, inspired from continuum mechanics with local criteria for the nucleation of plastic rearrangements [32], or more recently fed by molecular dynamics data concerning the details involved in the local criteria when they exist [33]. These models can indeed be shown as a good description of the atomistic modeling, but at a slightly larger coarse-graining length scale [34]. In these models, the effect of disorder is mainly located near the plastic rearrangements, and the elastic couplings are described within the framework of continuous linear isotropic and homogeneous elasticity. We will describe in this part the local plastic criteria and plastic events, then the numerical tools used for long-range elasticity in a contact mechanics problem, and finally the details of the contact conditions.

Local Plastic Criteria and Plastic Transformations
Plasticity is clearly related to an instability process [35], that is, a loss of mechanical stability due to the crossing of an energy barrier, and revealed by the loss of positive definiteness of the dynamical matrix of the full system (3N × 3N matrix of the second-order derivatives of the potential energy with N the number of atoms). In amorphous materials, such a loss of stability manifests itself by the nucleation of a high local strain concentration [13] that becomes unstable, and whose unfolding gives rise to a single Eshelby inclusion, or to a bidimensional alignment of Eshelby inclusions along elementary (embryonic) shear bands [12]. Although the instability process is a collective one (occurrence of a negative eigenvalue in the dynamical matrix of the full system), in the presence of disorder, the unstable mode can be described as a set of Eshelby inclusions, only one of them being the initiator of the mechanical instability. Thus, the mesoscopic models [31,32,36,37] are based on the nucleation of local Eshelby inclusions, interacting through long-range isotropic elasticity characteristic of the mechanical response of amorphous solids at large scale [38]. The nucleation criteria are still debated [33], since the instability process is indeed not local but global. However, a pragmatic analysis of the numerical data at the atomic scale shows that the nucleation of the local Eshelby inclusion may be related to large local strains, or low local elastic modulus. In the case of bulk metallic glasses, the plastic rearrangement has mainly a shear component. Thus, the nucleation criterion is related to the maximum local deviatoric strain, and naturally, the direction of the instability unfolding is the eigendirection related to this maximum deviatoric strain in the deviatoric part of the local strain tensor.
The implementation of the model is the following. The first two materials of Hertzian geometry are put in contact with a vertical load F z giving rise to a spherical contact area with a radius a (see Section 2.3). Starting from this initial configuration, a threshold ε d thresh is chosen as detailed later. If the maximum local deviatoric strain becomes larger than ε d thresh , then an Eshelby inclusion is nucleated at this place. ε d thresh is random with where ε 0 thresh is a homogeneous strain threshold, and r(i) is a Gaussian random number with zero average and variance unity, which is characteristic of the site i before rearrangement occurs. Q is a parameter: it is the amplitude of the random contribution to the local strain threshold. In the paper, we varied Q between 0 and 10 −1 (the maximum value is one order of magnitude larger than the maximum deviatoric strain induced by the Hertzian contact for the chosen load F z ). ε 0 thresh is chosen in order to allow for a progressive nucleation of plastic events, with only a limited number of sites. It is adjusted to the external load in order to get a given number of sites (5, 10, or 100) with a deviatoric strain larger than ε d thresh , that is to nucleate a fixed number of plastic events (5, 10, or 100 respectively) at each step. This allows imposing an average constant plastic strain rate. In practice, the sample is discretized into sites with lateral size r 0 corresponding to the size of the local rearrangements. r 0 is chosen to be equal to 10 −3 times the radius of curvature R at the contact (that is 1 nm for R = 1 µm). At each deformation step, the deviatoric strain will be computed on every site in the sample. The random value (-Q.r (i)) will be added at each site i to quantify the departure from the local plastic threshold. A fixed number of sites (5, 10, or 100) with the largest value of (ε d i − Q.r(i)) is identified, and these sites undergo a plastic rearrangement. After a plastic rearrangement occurred, a new value of r(i) is chosen at this site. The unfolding of the plastic rearrangements thus nucleated gives rise to a spherical Eshelby inclusion with radius r 0 . After the nucleation sites have been identified, an eigenstrain ε*(i) is attributed to each site i. The strain fields are computed again, taking into account the presence of the (5, 10, or 100) newly nucleated Eshelby inclusions, and the procedure is repeated with the calculation of the largest deviatoric strains-Q.r (i) in the new configuration. The method is summarized in Figure 1 (right).
The residual plastic strain inside the inclusion is the eigenstrain ε*. In the absence of disorder, the residual plastic strain is proportional to (ε d i − ε d thresh ) (see Figure 1 (left)) with ε d i the deviatoric strain measured at site i. In the presence of disorder, the eigenstrain ε* is chosen as with r*(i) another gaussian random number with zero average and variance unity, and S and Q* are two parameters of the model. S quantifies the exceeding of the local threshold. It is chosen in the absence of disorder in order to get similar values in the maximum deviatoric strain after the first nucleations of plastic rearrangements (too small values will stop the plastic deformation in the absence of disorder, and too large values will drive too much plastic rearrangement in the next step). For example, in the case of the contact problem described in Section 2.3, the values of the maximum deviatoric strain obtained after a given number n sites of plastic events have been nucleated simultaneously in the first step, as summarized in Table 1. It appears in Table 1 that the value S = 0.10 allows staying at the same values of deviatoric strains as in usual mesoscopic models [28] and will thus be kept in the following. Thus, Q and Q* are the only remaining parameters of the model. They will allow discussing the role of disorder in the organization of plastic events in our system. plastic rearrangement. After a plastic rearrangement occurred, a new value of r(i) is chosen at this site. The unfolding of the plastic rearrangements thus nucleated gives rise to a spherical Eshelby inclusion with radius r0. After the nucleation sites have been identified, an eigenstrain ε*(i) is attributed to each site i. The strain fields are computed again, taking into account the presence of the (5, 10, or 100) newly nucleated Eshelby inclusions, and the procedure is repeated with the calculation of the largest deviatoric strains-Q.r (i) in the new configuration. The method is summarized in Figure 1 (right).
The residual plastic strain inside the inclusion is the eigenstrain ε*. In the absence of disorder, the residual plastic strain is proportional to ( − ) (see Figure 1 (left)) with the deviatoric strain measured at site i. In the presence of disorder, the eigenstrain ε* is chosen as with r*(i) another gaussian random number with zero average and variance unity, and S and Q* are two parameters of the model. S quantifies the exceeding of the local threshold. It is chosen in the absence of disorder in order to get similar values in the maximum deviatoric strain after the first nucleations of plastic rearrangements (too small values will stop the plastic deformation in the absence of disorder, and too large values will drive too much plastic rearrangement in the next step). For example, in the case of the contact problem described in Section 2.3, the values of the maximum deviatoric strain obtained after a given number nsites of plastic events have been nucleated simultaneously in the first step, as summarized in Table 1. Table 1. Values of the maximum local deviatoric strain ε d 1 in a contact problem without disorder after the first step, compared to the initial value of the maximum deviatoric strain ε d 0. The contact problem is described in Section 2.3, here for different values of the parameter S and different numbers nsites of simultaneously nucleated sites. nsites S = 0.01 S = 0.10 S = 1.00 It appears in Table 1 that the value S = 0.10 allows staying at the same values of deviatoric strains as in usual mesoscopic models [28] and will thus be kept in the following. Thus, Q and Q* are the only remaining parameters of the model. They will allow discussing the role of disorder in the organization of plastic events in our system.

A Semi-Analytical Method for Long-Range Elasticity
We face to solve here the calculation of elastic strain in a contact problem and in the presence of Eshelby inclusions. The contact problem is first considered as a purely elastic contact between two linear elastic bodies with neither friction nor adhesion. One surface is flat, while the other has a radius of curvature R: it is a Hertz contact [39]. The calculation of the strain field beneath the contact can be performed with the help of superposition theorem and green functions for any given or computed pressure applied at the surface (similar to a Hertzian pressure) [40]. The Eshelby inclusions [41,42] can be computed thanks to the use of images symmetric to the interface, in order to still fulfill the limit boundary conditions [43].
The method was implemented in the software ISAAC (Version 2020, LaMCoS, Villeurbanne, France) [44] developed by D. Nélias et al. and adapted by T. Chaise, among others, for the calculation of residual plastic strain [45]. The semi-analytical method proposed to solve the contact problem is based on Eshelby's formalism but uses 3D Fast Fourier Transforms to speed up the computation. Thus, the time and memory necessary are greatly reduced in comparison with the classical finite element method. It allows a realistic computation of the stress and strain fields, taking into account the exact boundary conditions. Here, only an elastic contact is studied, but the same method can be extended to frictional contacts [46] and to adhesive contacts [47]. The only approximation used is that the only interface that matters is the contact surface: the bodies are considered as semi-infinite and treated within the framework of linear elasticity. The Eshelby inclusions and the corresponding heterogeneous solution are described using enrichment fields that are superimposed to the homogeneous problem [43]. To perform the simulations, a standard Young modulus E = 115 GPa and Poisson's ratio ν = 0.3 is used for both materials. These values correspond to the elastic moduli of Ti-based and Fe-based metallic glasses, the range of Young moduli spanning generally from 40 (La-based glasses) to 150 GPa approximately [1], and the Poisson's ratio being between 0.3 and 0.4 in bulk metallic glasses. The discretization of space is made through N = 65,559 nodes, which are separated by a distance dx = dy = 2r 0 horizontally and dz = r 0 vertically. Note that our goal here is not to describe precisely a given material but to highlight the disorder-induced mechanisms of localization of plastic deformation in a contact geometry.

Contact Problem
The initial contact problem is represented in Figure 2. It is made of a flat surface with Young modulus E and Poisson ration ν in contact with a spherical surface having a radius of curvature R and the same elastic properties. The bodies are pressed together using a normal force F z . Only the strain field inside the flat surface is studied here. The scaling relations in the Hertz contact are the following [39]: The radius of the spherical contact surface is The vertical displacement in the center is The maximum shear stress appears below the surface, at a distance z max ≈ 0.48 a for ν = 0.3, and its value is corresponding to a maximum deviatoric strain In our case, the initial force is imposed, thus fixing the initial contact radius , and then the deformation is induced by the succession of plastic rearrangements.
In the case of a purely elastic contact, the initiation of plasticity cannot appear at the surface, because the maximum shear stress is below the surface. In the absence of structural heterogeneities, the nucleation of shear bands is necessarily homogeneous in this case and located close to the depth . However, as soon as some adhesion or a sufficiently high static friction coefficient (µ ≥ 0.3) is present, then the shear stress becomes maximum at the surface [39,48]. Thus, the nucleation of plasticity and of shear banding is heterogeneous (i.e., from the surface) in these cases.
In order to allow heterogeneous nucleation or to modelize the effect of structural disorder, we have compared two cases: the case of isotropic and homogeneous linear elastic solids, and the case where two Eshelby inclusions are already located, close to the boundary of the contact, that is at a distance from the center of the contact zone. In this latter case, in order to get a real competition between the nucleation in the center and nucleation at the surface, the Eshelby inclusions have a shear eigenstrain slightly larger than .

Results
We will now compare the shape in the organization of the successive plastic events taking place in the flat solid for the different parameters studied.

Homogeneous vs. Heterogeneous Nucleation
The deviatoric stress induced by the presence of pre-existing inclusions competes with the deviatoric stress due to the Hertz contact. Figure 3 gives an idea of the spatial distribution of deviatoric strain (that is, the proximity to the local plastic yield strain assumed constant here as a first step) for different strengths of the two pre-existing defects. From left to right, the role of the two defects on the strain field increases. In the last case, In our case, the initial force is imposed, thus fixing the initial contact radius a, and then the deformation is induced by the succession of plastic rearrangements.
In the case of a purely elastic contact, the initiation of plasticity cannot appear at the surface, because the maximum shear stress τ max is below the surface. In the absence of structural heterogeneities, the nucleation of shear bands is necessarily homogeneous in this case and located close to the depth z max . However, as soon as some adhesion or a sufficiently high static friction coefficient (µ ≥ 0.3) is present, then the shear stress becomes maximum at the surface [39,48]. Thus, the nucleation of plasticity and of shear banding is heterogeneous (i.e., from the surface) in these cases.
In order to allow heterogeneous nucleation or to modelize the effect of structural disorder, we have compared two cases: the case of isotropic and homogeneous linear elastic solids, and the case where two Eshelby inclusions are already located, close to the boundary of the contact, that is at a distance a from the center of the contact zone. In this latter case, in order to get a real competition between the nucleation in the center and nucleation at the surface, the Eshelby inclusions have a shear eigenstrain slightly larger than ε d max .

Results
We will now compare the shape in the organization of the successive plastic events taking place in the flat solid for the different parameters studied.

Homogeneous vs. Heterogeneous Nucleation
The deviatoric stress induced by the presence of pre-existing inclusions competes with the deviatoric stress due to the Hertz contact. Figure 3 gives an idea of the spatial distribution of deviatoric strain (that is, the proximity to the local plastic yield strain assumed constant here as a first step) for different strengths of the two pre-existing defects. From left to right, the role of the two defects on the strain field increases. In the last case, where the deviatoric strain field at the defects is 10 times stronger than ε d max , the plastic deformation should localize on the defects, but in the two first cases with lower strain field, a connection is expected between the pre-existing defects at the surface and the subsurface zone of high shear strain due to the strain gradients induced by the Hertz contact. Figure 4 shows the position of the first plastic events in three cases: without defects, with defects comparable to the Hertz fields, with defects far stronger (10.1 times) than the Hertz fields. In the first two cases, the initial plastic events tend to localize all near z max ≈ 0.48 a, and there is no reason for the occurrence of shear bands. However, in the case of very strong defects, the simultaneous accumulation of plastic events near the defects and at z max suggests a future connection between them, for example along plastic bands. Clearly, from Figure 4, it appears that when the defects are not sufficiently strong, the plastic activity remains grouped below the surface at z max , while the presence of strong defects favors shear banding between z max and the surface. Thus, in the following, we will always keep the pre-existing defects at their highest value (with ε xy = 10.1 ε d max ).
Here, we have not taken into account the effect of disorder neither in the local yield strain nor on the eigenstrain characterizing the core of the unfolded inclusion. We will now discuss the effect of disorder on the strain fields and consequently on the nucleation of plastic events in the metallic glasses. where the deviatoric strain field at the defects is 10 times stronger than , the plastic deformation should localize on the defects, but in the two first cases with lower strain field, a connection is expected between the pre-existing defects at the surface and the subsurface zone of high shear strain due to the strain gradients induced by the Hertz contact. Figure 4 shows the position of the first plastic events in three cases: without defects, with defects comparable to the Hertz fields, with defects far stronger (10.1 times) than the Hertz fields. In the first two cases, the initial plastic events tend to localize all near ≈ 0.48 , and there is no reason for the occurrence of shear bands. However, in the case of very strong defects, the simultaneous accumulation of plastic events near the defects and at zmax suggests a future connection between them, for example along plastic bands. Clearly, from Figure 4, it appears that when the defects are not sufficiently strong, the plastic activity remains grouped below the surface at , while the presence of strong defects favors shear banding between zmax and the surface. Thus, in the following, we will always keep the pre-existing defects at their highest value (with = 10.1 ). Here, we have not taken into account the effect of disorder neither in the local yield strain nor on the eigenstrain characterizing the core of the unfolded inclusion. We will now discuss the effect of disorder on the strain fields and consequently on the nucleation of plastic events in the metallic glasses.   where the deviatoric strain field at the defects is 10 times stronger than , the plastic deformation should localize on the defects, but in the two first cases with lower strain field, a connection is expected between the pre-existing defects at the surface and the subsurface zone of high shear strain due to the strain gradients induced by the Hertz contact. Figure 4 shows the position of the first plastic events in three cases: without defects, with defects comparable to the Hertz fields, with defects far stronger (10.1 times) than the Hertz fields. In the first two cases, the initial plastic events tend to localize all near ≈ 0.48 , and there is no reason for the occurrence of shear bands. However, in the case of very strong defects, the simultaneous accumulation of plastic events near the defects and at zmax suggests a future connection between them, for example along plastic bands. Clearly, from Figure 4, it appears that when the defects are not sufficiently strong, the plastic activity remains grouped below the surface at , while the presence of strong defects favors shear banding between zmax and the surface. Thus, in the following, we will always keep the pre-existing defects at their highest value (with = 10.1 ). Here, we have not taken into account the effect of disorder neither in the local yield strain nor on the eigenstrain characterizing the core of the unfolded inclusion. We will now discuss the effect of disorder on the strain fields and consequently on the nucleation of plastic events in the metallic glasses.    Figure 5 shows the effect of the two parameters S and Q* on the position of plastic events after 250 plastic rearrangements took place (25 steps with 10 rearrangements at each step). In these figures, there is no disorder in the yield strain (Q = 0), and the pre-existing defects are 10 times stronger than the Hertz field, as detailed in the previous part. It appears clearly from this figure that the parameter S, which is used to counterbalance the lack of precision on the determination of the local yield strain and allow for a better definition of plastic flow, induces spatial correlations in the location of plastic events: the deviatoric strain is higher close to the Eshelby inclusions, thus giving rise to higher eigenstrains that will contribute to self-sustain the process. The location of the plastic events shows a better alignment close to the pre-existing defects, and it looks more sparse close to z max . As discussed in Section 2.1 for the disorder-free case, the value of S = 0.1 will be kept in the following.

Role of Disorder in the Residual Plastic Eigenstrain
Let's look now at the effect of the disorder through the parameter Q*. In Figure 6, the distribution of deviatoric strains in the sample is shown for different values of Q*. Increasing Q* not only contributes to enlarging the distribution of deviatoric strains as compared to the disorder-free Hertz case, but it can even give rise to very large cumulated deviatoric strains, up to 10% strain for Q* = 0.01 at the 25th step (Figure 6 bottom). Let us remind that Q* = 0.01 corresponds to a Gaussian distribution of eigenstrains with variance Q* = 1%, owing in our case at the 25th step extreme eigenstrains of about 3.5% only (Figure 6 top). In Figure 5, increasing Q* up to the 1% strain clearly contributes to counterbalance the role of the Hertz fields and of their spatial gradients. The plastic rearrangements localize close to the pre-existing defects. The resulting shear aligns along large shear bands, propagating straightly from the surface defects. While the role of defects is clearly enhanced in this case, it does not gives rise to curved shear bands, but only to straight bands initiated from the surface, even after a large number of events.
surface shown as Eshebly inclusions with the shear eigenstrain = 10.1 (right). Bottom: at step 10 of the deformation (nsites = 10) with pre-existing defects with the shear eigenstrain = 1.01 (left) and = 10.1 (right). Figure 5 shows the effect of the two parameters S and Q* on the position of plasti events after 250 plastic rearrangements took place (25 steps with 10 rearrangements a each step). In these figures, there is no disorder in the yield strain (Q = 0), and the pre existing defects are 10 times stronger than the Hertz field, as detailed in the previous par It appears clearly from this figure that the parameter S, which is used to counterbalanc the lack of precision on the determination of the local yield strain and allow for a bette definition of plastic flow, induces spatial correlations in the location of plastic events: th deviatoric strain is higher close to the Eshelby inclusions, thus giving rise to higher e genstrains that will contribute to self-sustain the process. The location of the plastic event shows a better alignment close to the pre-existing defects, and it looks more sparse clos to zmax. As discussed in Section 2.1 for the disorder-free case, the value of S = 0.1 will b kept in the following.

Role of Disorder in the Residual Plastic Eigenstrain
Let's look now at the effect of the disorder through the parameter Q*. In Figure 6, th distribution of deviatoric strains in the sample is shown for different values of Q*. Increas ing Q* not only contributes to enlarging the distribution of deviatoric strains as compare to the disorder-free Hertz case, but it can even give rise to very large cumulated deviatori strains, up to 10% strain for Q* = 0.01 at the 25th step (Figure 6 bottom). Let us remind tha Q* = 0.01 corresponds to a Gaussian distribution of eigenstrains with variance Q* = 1% owing in our case at the 25th step extreme eigenstrains of about 3.5% only ( Figure 6 top In Figure 5, increasing Q* up to the 1% strain clearly contributes to counterbalance the rol of the Hertz fields and of their spatial gradients. The plastic rearrangements localize clos to the pre-existing defects. The resulting shear aligns along large shear bands, propagatin straightly from the surface defects. While the role of defects is clearly enhanced in thi case, it does not gives rise to curved shear bands, but only to straight bands initiated from the surface, even after a large number of events.

Role of Disorder in the Local Plastic Thresholds
Let us now look at the effect of disorder in the yield strains. For low disorder Q = 10 −3 , large shear bands are initiated from the surface defects discussed before in Section 3.2. When Q increases, the nucleation of plastic events is focused on the surface defects, and some of them become sensitive to the Hertz fields nucleate preferentially close to zmax. When the disorder is very strong and far larger t the maximum deviatoric strain induced by the Hertz contact, the events are sp and can occur everywhere. However, the Hertz fields will introduce a bias that will fa a nucleation close to zmax. In this latter case, the plastic events do not organize clearly al shear bands, and more importantly, they do not show any sensitivity to the surface appears completely ignored in this organization.
Interestingly, between these two extreme cases, for a contribution of the disorde the local yield strain that is comparable to , a connection is made between ev nucleated at the surface, and events deeper below the surface. For Q ≈ 2 ≈ 10 , characteristic shape of curved shear bands initiated from the strongest defects at the face is recognizable.

Role of Disorder in the Local Plastic Thresholds
Let us now look at the effect of disorder in the yield strains. Figure 7 summarizes the effect of the disorder induced by the increase in the parameter Q in the random part of the local nucleation threshold of Eshelby inclusions. The colored map shows the local deviatoric strain -Q.r (i). Thus, the plastic events always take place at the brightest colors (maximum values) in this map. While for low values of Q, the gradients of the Hertz field are visible, they progressively disappear for large disorder values. The previous values of S = 0.10 and Q* = 0.01 are kept in this case.
For low disorder Q = 10 −3 , large shear bands are initiated from the surface defects, as discussed before in Section 3.2. When Q increases, the nucleation of plastic events is less focused on the surface defects, and some of them become sensitive to the Hertz fields and nucleate preferentially close to z max . When the disorder is very strong and far larger than the maximum deviatoric strain ε d max induced by the Hertz contact, the events are sparse and can occur everywhere. However, the Hertz fields will introduce a bias that will favor a nucleation close to z max . In this latter case, the plastic events do not organize clearly along shear bands, and more importantly, they do not show any sensitivity to the surface that appears completely ignored in this organization.
Interestingly, between these two extreme cases, for a contribution of the disorder to the local yield strain that is comparable to ε d max , a connection is made between events nucleated at the surface, and events deeper below the surface. For Q ≈ 2ε d max ≈ 10 −2 , the characteristic shape of curved shear bands initiated from the strongest defects at the surface is recognizable. Metals 2021, 10, x FOR PEER REVIEW 10 of 16 Figure 7. Effect of Q characterizing the disorder in the yield strain from Equation (1), on the localization of plastic rearrangements, for different number of steps, with Q* = 0.01 and S = 0.100 and initial defects with = 10.1 . With our choice of loading Fz, the maximum deviatoric strain due to Hertz contact in the absence of disorder is equal to 0.65 × 10 −2 here.

Discussion
We have studied here with the help of a mesoscopic model the spatial organization of plastic events inside a disordered material in a contact geometry. The results obtained in this context underline the fact that the sliding conditions at the surface play an important role on the plasticity pattern when the structural disorder is sufficiently low. The role of adhesion or friction at the interface is modeled with the help of pre-existing inclusions located near the free surface of the solid and at the boundary of the contact zone. Indeed, adhesion or friction with small tangential load induces very strong deviatoric strains at the boundary of the contact [39]. It has been shown ( Figure 4) that the presence of initial defects-or equivalently non sliding conditions-is crucial to force the emergence of shear bands at the surface, as observed experimentally in many cases [27][28][29]49]. These pre-existing defects are also representative of the possible role played by a very inhomogeneous structural disorder or by surface defects. This is in agreement with the possibility to generate extrinsic nucleation of shear bands [20], and this enhanced surface sensitivity can be used for specific materials design [26], for example to prevent subsurface damage. We have discussed in the present paper the role of the eigenstrains included  (1), on the localization of plastic rearrangements, for different number of steps, with Q* = 0.01 and S = 0.100 and initial defects with ε xy = 10.1 ε d max . With our choice of loading F z , the maximum deviatoric strain ε d max due to Hertz contact in the absence of disorder is equal to 0.65 × 10 −2 here.

Discussion
We have studied here with the help of a mesoscopic model the spatial organization of plastic events inside a disordered material in a contact geometry. The results obtained in this context underline the fact that the sliding conditions at the surface play an important role on the plasticity pattern when the structural disorder is sufficiently low. The role of adhesion or friction at the interface is modeled with the help of pre-existing inclusions located near the free surface of the solid and at the boundary of the contact zone. Indeed, adhesion or friction with small tangential load induces very strong deviatoric strains at the boundary of the contact [39]. It has been shown ( Figure 4) that the presence of initial defects-or equivalently non sliding conditions-is crucial to force the emergence of shear bands at the surface, as observed experimentally in many cases [27][28][29]49]. These pre-existing defects are also representative of the possible role played by a very inhomogeneous structural disorder or by surface defects. This is in agreement with the possibility to generate extrinsic nucleation of shear bands [20], and this enhanced surface sensitivity can be used for specific materials design [26], for example to prevent subsurface damage. We have discussed in the present paper the role of the eigenstrains included in these pre-existing defects, and we focused mainly on the case of defects whose eigenstrain is about 10 times larger than the maximum deviatoric strain induced by the Hertz contact beneath the surface. Indeed, this high value appears necessary to avoid a concentration of plasticity only below the surface. Moreover, we have shown ( Figure 5) that when the same amount of disorder appears as well in the local strain rearrangements (value of Q* in the eigenstrains), then plasticity localizes along straight lines initiated from the initial defects. However, this type of disorder is not sufficient to induce curved shear bands, as sometimes observed [27][28][29]49,50].
The occurrence of curved shear bands is shown here to result only from a specific set of parameters in the random local strain thresholds for plasticity. To get well marked permanent shear bands (Figure 7), it is necessary to have an amount of disorder in the local strain threshold that is comparable to the maximum deviatoric strain due to the Hertz contact between curved surfaces. However, this one depends not only on the contact geometry (radius of curvature R) but also on the external load (Equation (6)). Thus, the occurrence of curved shear bands is not universal as initially affirmed [28] but depends on the comparison between intrinsic heterogeneities (local strain thresholds for plasticity) and extrinsic parameters (geometry, external load). It is very important here to note that the structural disorder in itself is not sufficient to predict the shape of shear bands but must be put in perspective of the loading conditions: size of the contact area, radius of curvature, and external load, as detailed in Equation (6). The respective role of the disorder in the eigenstrains of the unfolded defects and in the strain thresholds may also explain why the curved shape of initial shear bands is not maintained for higher loads in easily rearranging materials [30]. The disorder in the eigenstrain or some kind of self-sustainability [20][21][22] may finally dominate and move away from the initial range of parameters. This may also explain the variety of the plasticity patterns observed in metallic glasses under a nanoindentor [29].
These effects were not discussed in [28], aiming at proving the universal character of plastic yield in amorphous solids for any kind of load. However, indeed, many assumptions were made in [28] concerning the relative values of disorder. First, the yield strains and the eigenstrains are not chosen independently, but the eigenstrain is assumed to be proportional to the yield strain. Moreover, non-sliding contact, meaning very strong solid friction, has been assumed that favors nucleation at the surface and finite size effects, as discussed in our article.
The disorder in the strain thresholds is usually related to the structural disorder characteristic of amorphous materials [33,51,52]. Experimentally, the amount of disorder depends on the preparation protocol, such as the quenching rate or the annealing for example [53,54], the free volume [9,55], or even defects induced by neutron irradiation [56]. The amount of disorder is known to affect the elastic and the plastic properties of the glass [57]. In all cases, increasing the effective temperature of the glass, or the amount of soft-equivalently less relaxed-zones [58] will make the glass more ductile, less fragile, and contribute to homogenize the deformation. In case of very strong disorder, the role of surface defects becomes negligible (Figure 7). Thus, the extreme values of strains are sensitive again to the Hertzian field. Consequently, the deformation, while irregular, concentrates below the surface, at the maximum deviatoric strain due to the Hertz contact. Such a kind of deformation without shear bands was already observed experimentally with in situ X-ray diffraction in Zr-based metallic glass samples [59]. However, the disorder was unfortunately not characterized in these samples. Thus, the absence of shear bands may also be simply due to perfectly sliding conditions at the surface.
In general, very different kinds of plasticity patterns have been observed below an indentor. Perepezko et al. [60] have identified experimentally on the same sample at least two kinds of events: at the surface and below the surface. Luo et al. [6] have observed different plasticity patterns depending on the amount of free volume induced by chemical strengthening in silicate glasses. The decrease in the free volume always favors localization along shear banding; this may result from a lower quenching rate [29] or the increase in the relaxation induced by longtime annealing, for example [51,53]. This effect also appears in our analysis of the role of the amplitude of disorder in the local strain thresholds: a strong disorder will make every site equivalent, thus masking the role of surface or isolated defects and preventing the nucleation of shear bands. Then, the average mechanical response becomes that of a homogeneous sample.
Note that the size of the samples is shown in experiments to play a role on the ductile behavior of bulk metallic glass samples [61,62]. In general, diminishing the size of the sample, as in films or in micropillars, will induce a more homogeneous repartition of plastic deformation. In this case, the increased role of surface defects (due to the large surface over volume ratio of the samples) may be counterbalanced by the lack of disorder or low effective temperature of the sample, due to its small size. However, on the contrary, the enhanced ductile behavior could also be due to an increase in the effective temperature due to the pressure increase and the confinement. As suggested in this paper, contact force experiments with different friction coefficients, or with different radius of curvature, could help disentangle these explanations by looking at the shear band patterns.
In our mesoscopic model, the formation of shear bands results only from the succession of discrete Eshelby events. Instantaneous shear bands (also called elementary shear bands [17] or embryonic shear bands [20]) are not properly described here, as it is the case in all mesoscopic models [31]. Elementary shear bands result indeed from a specific outof-equilibrium unfolding process with a collective kinetic alignment of events. However, mesoscopic models are able to reproduce some characteristics of these bands: the resulting spatial organization of plastic deformation is finally not very far from quasi-static atomistic simulations [51]. This may be due to the fact that the driving force in the unfolding of the instability is the elastic force, as in the mesoscopic models. However, kinetic effects may play a role. For example, inertial effects could affect the results, as already discussed in [58]. The goal of the present paper is to focus on the primary role of disorder. To compare different shear rates, we have imposed a different number of simultaneous plastic events at each step in the model before redistributing the elastic energy. This does not induce any significant difference in the results, at least in the low disorder case. Thus, we have chosen to keep, in the majority of the simulations presented here, an intermediate number of simultaneous plastic events corresponding to a low shear rate. However, the role of the shear rate should be deepened, since it is known to induce different spatial organizations of events [63]. The role of inertia and of shear rates should deserve further studies if needed.

Conclusions
The choice of the approximate mesoscopic model used in this study allows keeping the main ingredients for the nucleation and succession of plastic events in an amorphous material: amorphous materials are linear elastic at large scale, with homogeneous and isotropic elasticity [38], and plastic rearrangements are described as Eshelby inclusions in agreement with recent molecular dynamics simulations [12].
In this work, two important effects were evidenced that have never been discussed in previous studies: the first one is the role of initial defects, and the second one is the relative role of disorder vs. Hertz mechanical inhomogeneities. Concerning the role of initial defects, it is clear that if these defects are absent, then the plastic deformation will be self-sustained below the contact at the maximum deviatoric strain location due to the Hertzian contact, as in homogeneous solids. When the deviatoric strain supported by the pre-existing defects becomes comparable or larger to the maximum Hertzian deviatoric strain, then the extrinsic nucleation of shear bands occurs, and new processes take place. We saw in Section 3.2 that the role of initial defects is reinforced when the residual plastic strain (eigenstrain of the nucleated Eshelby inclusions) has a random component. However, this type of disorder is not sufficient to induce well-defined curved shear bands inside the sample. It only contributes to enlarging the plastic zone around an initial one. We saw also in the simple model used here that the sensitivity of the residual plastic strain on the amount of deviatoric strain at the nucleation of the plastic event induces a small memory effect, tending to localize the plastic activity around its initial place. Concerning the role of disorder, we compared different sources of disorder in the amorphous sample.
Interestingly, we saw that the spatial organization of the plastic activity and the connection between different places inside the sample (that would correspond to extended shear bands) depend mainly on the amount of disorder in the local yield strain only. This can be understood by the fact that the high strain disorder due to structural origins should counterbalance the strain inhomogeneities induced by the external load. For very high disorder, the differences in extreme values responsible for the nucleation of plastic events will be reinforced by Hertzian inhomogeneities and become certainly insensitive to surface effects. Two simple cases have been easily identified: the case of low disorder where pre-existing defects will control the nucleation of shear bands (initiated at the surface when these defects represent adhesive or frictional contact). In this case, the shear band extends straight inside the sample. The second limit case is the case of very high disorder where plasticity is sparse and insensitive to surface defects. At intermediate values of disorder, more precisely when the contribution of disorder in the yield strain becomes comparable to the maximum Hertzian deviatoric strain, well-defined curved shear bands connecting the surface and the subsurface are shown. This latter case occurs for different amounts of disorder, depending on the geometry and on the external load: as shown in Equation (3), increasing the load or decreasing the radius of curvature at the contact increase the deviatoric strain to which the amount of disorder in the strain threshold must be compared. However, the disorder that is needed to nucleate shear bands and inhomogeneous deformation has less influence on strain when the contact radius of curvature is smaller. Note that in the absence of surface defects, the plastic activity also remains located at the place of maximum Hertzian deviatoric strain without extended shear banding.
These different cases were all observed in experiments: less relaxed, fast quenched, or irradiated samples giving rise to more homogeneous plastic deformation, as observed here in the strong disorder case. However, it is difficult to have a more quantitative estimation of the degree of structural disorder needed. Indeed, in bulk metallic glasses, species are mixed, and it is quite difficult in general to characterize the disorder [53], unless in very specific cases [64]. In sodo-silicate glasses, at least the nanoporous silicon skeleton and channels of mobile ions give some keys of spatial characterization of the mattering structural disorder [65][66][67].
To conclude, boundary conditions, and also the relative amount of intrinsic strain disorder (to some extend in the residual plastic strain, but more importantly in the local yield strain), compared to the maximum Hertzian deviatoric subsurface strain, have a significant impact on the spatial organization of plasticity in amorphous materials in a contact geometry. The boundary conditions as well as the Hertzian strain depend on the external load, making the mechanical behavior of glasses non-universal in a contact geometry.