O U Salman

R Baggio

B Bacroix

G Zanzotto
email: zanzotto@dmsa.unipd.it

N Gorbushin
email: nikolai.gorbushin@espci.fr

L Truskinovsky
email: lev.truskinovsky@espci.fr

Discontinuous yielding of pristine micro-crystals

Keywords: . plasticity, dislocations, Landau theory, nucleation, pattern formation, brittleness. Mathematical subject classification (2010). 00X99

We study the mechanical response of a dislocation-free 2D crystal under homogenous shear using a new mesoscopic approach to crystal plasticity, a Landau-type theory, accounting for the global invariance of the energy in the space of strain tensors while operating with an infinite number of equivalent energy wells. The advantage of this approach is that it eliminates arbitrariness in dealing with topological transitions involved, for instance, in nucleation and annihilation of dislocations. We use discontinuous yielding of pristine micro-crystals as a benchmark problem for the new theory and show that the nature of the catastrophic instability, which in this setting inevitably follows the standard affine response, depends not only on lattice symmetry but also on the orientation of the crystal in the loading device. The ensuing dislocation avalanche involves cooperative dislocation nucleation, resulting in the formation of complex microstructures controlled by a nontrivial selfinduced coupling between different plastic mechanisms.

Introduction

With the advance of nanotechnology and broad fabrication of nano-scale structures, the focus in the study of plastic deformation has shifted to atomic dimensions. The emerging science of nanomaterials deals, for instance, with machine parts printed by chemical vapor deposition or made of nano-grained metals. At these manufacturing scales, external and internal (microstructural) lengths become comparable, and the dislocation-based description of plasticity comes to the forefront in providing design guidelines for miniaturized mechanical devices [START_REF] Lee | Uniaxial compression of fcc Au nanopillars on an MgO substrate: The effects of prestraining and annealing[END_REF][START_REF] Li | A promising structure for fabricating high strength and high electrical conductivity copper alloys[END_REF][START_REF] Lu | Stabilizing nanostructures in metals using grain and twin boundary architectures[END_REF]. It was found, for instance, that sub-micron size objects, serving as components of such systems, are characterized by remarkably high strength. This opened a way to a broad range of novel engineering applications, including nano-metric machining and hierarchical steels, e.g. [START_REF] Schaedler | Ultralight metallic microlattices[END_REF].

Considerable research efforts have been focused on the study of dislocation plasticity in submicron crystals in an attempt to assess the extent of the failure of traditional inelastic constitutive models concerning phenomena at these scales [START_REF] Mordehai | Size effect in compression of single-crystal gold microparticles[END_REF][START_REF] Maaß | Ultrahigh strength of dislocation-free Ni3Al nanocubes[END_REF][START_REF] Han | From "Smaller is Stronger" to "Size-Independent Strength Plateau": Towards Measuring the Ideal Strength of Iron[END_REF][START_REF] Papanikolaou | Avalanches and plastic flow in crystal plasticity: an overview[END_REF][START_REF] Maaß | Micro-plasticity and recent insights from intermittent and small-scale plasticity[END_REF]. As a result of these efforts, it has become clear that the deformation mechanisms, which we habitually associate with plastic flows, change dramatically once the sample size is reduced below the micrometer range. It was found that the strength of such crystals reaches theoretical (ideal) limit [START_REF] Nix | Indentation size effects in crystalline materials: A law for strain gradient plasticity[END_REF][START_REF] Uchic | Sample dimensions influence strength and crystal plasticity[END_REF][START_REF] Greer | Size dependence of mechanical properties of gold at the micron scale in the absence of strain gradients[END_REF][START_REF] Dimiduk | Scale-free intermittent flow in crystal plasticity[END_REF] and that the plastic flow proceeds through stress drops reminiscent of brittle fracture [START_REF] Bei | Effects of pre-strain on the compressive stress-strain response of Mo-alloy single-crystal micropillars[END_REF][START_REF] Chrobak | Deconfinement leads to changes in the nanoscale plasticity of silicon[END_REF][START_REF] Wang | Pristine-to-pristine regime of plastic deformation in submicron-sized single crystal gold particles[END_REF][START_REF] Cui | Influence of loading control on strain bursts and dislocation avalanches at the nanometer and micrometer scale[END_REF]. The attendant intermittency compromises the reliable functioning of ultra-small machinery and jeopardizes our ability to ensure the predictable performance of MEMS and other similar systems [START_REF] Csikor | Dislocation Avalanches, Strain Bursts, and the Problem of Plastic Forming at the Micrometer Scale[END_REF][START_REF] Benzerga | Micro-pillar plasticity: 2.5D mesoscopic simulations[END_REF][START_REF]Plasticity of Micrometer-Scale Single Crystals in Compression[END_REF][START_REF] Argon | Strain avalanches in plasticity[END_REF][START_REF] Zhang | Taming intermittent plasticity at small scales[END_REF].

Of primary interest to our study will be the phenomenon of discontinuous yielding in sub-micron (initially) dislocation free volumes. According to the classical continuum paradigm, the elastic deformation in strained crystals must be followed by either an abrupt brittle-like failure or gradual plastic deformation. Instead, loading of sub-micron crystals revealed a mixed behavior. Microcrystals, exhibiting smooth classical yield at macro-scales, were shown to fail catastrophically, with conventional work hardening replaced by an abrupt system spanning plastic event [START_REF] Greer | Plasticity in small-sized metallic systems: Intrinsic versus extrinsic size effect[END_REF].

The precipitous brittle yield, revealing a high level of dislocation correlations and manifesting plastic collapse, was first discovered in metal whiskers [START_REF] Brenner | Tensile Strength of Whiskers[END_REF][START_REF]Growth and properties of "whiskers[END_REF]. Brittleness of this type is now routinely observed in nanoparticles, which "break plastically" while generating a large number of dislocations [START_REF] Sharma | Nickel nanoparticles set a new record of strength[END_REF][START_REF] Mordehai | Nucleation-Controlled Plasticity of Metallic Nanowires and Nanoparticles[END_REF]; similar global plastic instabilities have also been recorded during nanoindentation [START_REF] Lilleodden | Microstructural length-scale effects in the nanoindentation behavior of thin gold films[END_REF][START_REF] Corcoran | Anomalous plastic deformation at surfaces: Nanoindentation of gold single crystals[END_REF]. The implied system size plastic avalanches result from a large number of highly cooperative individual dislocation nucleation and dislocation glide events, taking place almost simultaneously [START_REF] Dimiduk | Scale-free intermittent flow in crystal plasticity[END_REF][START_REF] Csikor | Dislocation Avalanches, Strain Bursts, and the Problem of Plastic Forming at the Micrometer Scale[END_REF][START_REF] Maaß | Small-scale plasticity: Insights into dislocation avalanche velocities[END_REF][START_REF] Lee | In-situ observation of the initiation of plasticity by nucleation of prismatic dislocation loops[END_REF].

The distinctive features of crystal plasticity at sub-micron scales can be attributed to the high degree of structural perfection of such samples, which are usually almost pristine [START_REF] He | In situ observation of shear-driven amorphization in silicon crystals[END_REF][START_REF] Merabet | Low-temperature intrinsic plasticity in silicon at small scales[END_REF]. In particular, the catastrophic yield depends critically on the absence of dislocations before straining. Despite the global, system-spanning scale of experienced plastic deformation, such samples produce nearly pristine postmortem microstructures. Almost all dislocations mediating the explosive stress drop manage to escape to the boundaries, sometimes even producing pristine-to-pristine transition [START_REF] Wang | Pristine-to-pristine regime of plastic deformation in submicron-sized single crystal gold particles[END_REF][START_REF] Chisholm | Dislocation starvation and exhaustion hardening in Mo alloy nanofibers[END_REF]. It was also observed that as the deformation volume of a material decreases, the effect of crystal orientation on the operative deformation mechanisms increases; in particular, the mechanical response of sub-micron defect-free nanopillars is different when they are deformed along high symmetry orientations and low symmetry orientations [START_REF] Ziegenhain | Influence of crystal anisotropy on elastic deformation and onset of plasticity in nanoindentation: A simulational study[END_REF][START_REF]Effect of crystal orientation on the size effects of nano-scale fcc metals[END_REF]. The compression tests on confined micropillars also revealed higher (and therefore closer to theoretical limit) stresses of massive dislocation nucleation than in the presence of free surfaces [START_REF]Effect of crystal orientation on the size effects of nano-scale fcc metals[END_REF][START_REF]The effect of crystal anisotropy and pre-existing defects on the incipient plasticity of FCC single crystals during nanoindentation[END_REF].

In this paper, addressing the mechanical response of dislocation-free 2D crystals under homogenous shear, we attempt to rationalize the above observations using the mesoscopic model of crystal plasticity [START_REF] Salman | Minimal integer automaton behind crystal plasticity[END_REF][START_REF] Salman | On the critical nature of plastic flow: One and two dimensional models[END_REF][START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF]. We use discontinuous yielding of pristine micro-crystals as a benchmark problem for this new theory and study in detail how the nature of the catastrophic instability depends on the lattice symmetry and the orientation of the crystal in the loading device. The initial plasticity in pristine nano-scale volumes is of particular theoretical importance as the process dominated by cooperative dislocation nucleation vs more conventional flow mechanisms involving multiplication and glide of pre-existing dislocations.

The rest of the paper is organized as follows. In Section 2, we further motivate our study reviewing the problem from the viewpoint of materials science of sub-micron crystals. We then briefly discuss the existing computational tools. The new mesoscopic model is introduced in Section 3, where we recall the geometrically nonlinear kinematics of crystal lattices and discuss the construction of strain-C. R. Physique, 0000, 1, n o 0, 000-000 energy density respecting the finite strain symmetry of the crystal lattice. In Section 4, we present the macroscopic picture of the discontinuous yield revealed by our numerical experiments and, to justify these results, compute the theoretical (ideal) shear strength for perfect crystals with different point groups. In Section 5, we study post-bifurcation behavior of pristine crystals and compare the representations of post avalanche dislocation patterns. Along with their representation in the physical space, we also visualize such patterns in the configurational space of lattice strains. In Section 6, we briefly discuss the mechanical response of the simulated crystals beyond the catastrophic brittle event. Finally, in Section 7 we present our conclusions.

Some background

Plastic deformation in crystalline solids involves discontinuous changes in the configuration of nearest neighbors as some stress/strain thresholds are exceeded upon external loading. The associated relative movement of atoms can be interpreted as the creation and propagation of topological lattice defects. The most prominent among them are dislocations, moving along crystallographic planes and leaving behind quantized lattice slip [START_REF] Read | Dislocations in Crystals[END_REF][START_REF] Cottrell | Commentary. A brief view of work hardening[END_REF][START_REF] Kubin | Dislocations, mesoscale simulations and plastic flow[END_REF]. Understanding the collective motion of dislocations is the key to control the ductile failure of crystals as these linear defects are crucially involved in shear band formation, fatigue, and even fracture [START_REF] Differt | A model of extrusions and intrusions in fatigued metals II. Surface roughening by random irreversible slip[END_REF][START_REF] Antolovich | Plastic strain localization in metals: origins and consequences[END_REF][START_REF] Weiss | Plastic intermittency during cyclic loading: From dislocation patterning to microcrack initiation[END_REF].

Dislocation motion is inherently complex due to long-range elastic interactions of dislocation lines and strongly nonlinear, threshold type short-range interactions of dislocation cores [START_REF] Madec | From dislocation junctions to forest hardening[END_REF][START_REF] Sethna | Deformation of Crystals: Connections with Statistical Physics[END_REF]. Facilitated by lattice trapping, transient and sessile dislocation patterns cover a broad range of scales from microscopic (junctions) to macroscopic (grains) [START_REF] Gómez-García | Dislocation Patterns and the Similitude Principle: 2.5D Mesoscale Simulations[END_REF][START_REF] Chen | Bending crystals: emergence of fractal dislocation structures[END_REF][START_REF] Li | Fundamental factors on formation mechanism of dislocation arrangements in cyclically deformed fcc single crystals[END_REF].

Despite the presence of microscopic heterogeneities, plastic flows appear at the scale of bulk materials as smooth phenomena amenable to continuous description [START_REF] Takeuchi | Work Hardening of Copper Single Crystals with Multiple Glide Orientations[END_REF]. Therefore, macroscopic crystal plasticity is usually formulated within classical continuum mechanics and involves macroscale constitutive relations. The inelastic component of strain tensor is parametrized by a finite number of order parameters representing amplitudes of pre-designed plastic mechanisms; in rate independent limit each of them is assumed to be governed by dry friction dynamics [START_REF] Han | Mechanism-based strain gradient crystal plasticity -I. Theory[END_REF][START_REF] Gurtin | Thermodynamics applied to gradient theories involving the accumulated plastic strain: The theories of Aifantis and Fleck and Hutchinson and their generalization[END_REF][START_REF] Roters | Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications[END_REF][START_REF]Crystal Plasticity and Evolution of Polycrystalline Microstructure[END_REF][START_REF] Forest | Micromorphic Crystal Plasticity[END_REF][START_REF] Salman | Tempering the mechanical response of FCC micro-pillars: an Eulerian plasticity approach[END_REF]. The fact that fluctuations are effectively averaged out 1 , opens for continuum plasticity (CP) access to macroscopic time and length scales and allows one to model realistic 3D structures with complex geometries while accounting phenomenologically for such complex effects as hardening, rate-dependence and even polycrystallinity [START_REF] Franciosi | The concepts of latent hardening and strain hardening in metallic single crystals[END_REF][START_REF] Zhang | Microstructure-based crystal plasticity modeling of cyclic deformation of Ti-6Al-4V[END_REF][START_REF] Roters | Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications[END_REF][START_REF] Forest | Generalized Continua and Phase-Field Models: Application to Crystal Plasticity[END_REF][START_REF] Marano | Intragranular localization induced by softening crystal plasticity: Analysis of slip and kink bands localization modes from high resolution FFT-simulations results[END_REF].

However, it has been recently realized that plastic flows in ultra-small samples are beyond the reach of such theories. A description of the cross-over from a 'mild' plastic flow in bulk materials to a 'wild' scale-free intermittent plastic response at the sub-micron scale requires significant paradigm change [START_REF] Weiss | From mild to wild fluctuations in crystal plasticity[END_REF]. In particular, the current CP theory fails to resolve intermittent stress-drops (or strain-jumps) at submicron scales [START_REF] Uchic | Sample dimensions influence strength and crystal plasticity[END_REF][START_REF] Dimiduk | Scale-free intermittent flow in crystal plasticity[END_REF][START_REF] Zhang | Taming intermittent plasticity at small scales[END_REF] and is also unable to explain catastrophic events accompanying plastic flows in nanoparticles [START_REF] Wang | Pristine-to-pristine regime of plastic deformation in submicron-sized single crystal gold particles[END_REF], nanowires [START_REF] Lu | Surface dislocation nucleation mediated deformation and ultrahigh strength in sub-10-nm gold nanowires[END_REF] and nanopillars [START_REF] Bei | Effects of pre-strain on the compressive stress-strain response of Mo-alloy single-crystal micropillars[END_REF][START_REF] Issa | In situ investigation of MgO nanocube deformation at room temperature[END_REF][START_REF] Hu | Dislocation avalanche mechanism in slowly compressed high entropy alloy nanopillars[END_REF][START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF]. None of these phenomena can be rationalized without a direct reference to dislocation motion.

If in bulk materials dislocation motion is largely uncorrelated and can be indeed averaged out, plastic flows in sub-micron crystal involve highly cooperative dislocational rearrangements. The most striking effect of such cooperativity is provided by intermittent system size dislocation avalanches, which defy self-averaging and challenge any attempt of continuous description. More generally, the scale-free CP fails at sample sizes comparable to the emerging microstructural scale of defect patterning, simply because the latter is assumed to be zero in the classical theory. The attempts to regularize the CP and link the internal length scale with the presence of 'geometrically necessary' strain gradients have been so far only partially successful as the size effect was also observed in the absence of strain gradients, e.g. [START_REF] Bittencourt | Interpretation of the size effects in micropillar compression by a strain gradient crystal plasticity theory[END_REF].

An alternative approach to the rationalization of the size effect was to shift the attention from the role of gradients in dislocation arrangements to the scarcity of dislocation sources in small crystals. It was noticed that when sufficiently small crystals were strained to the bulk yield stresses, the rate of dislocation escape to the surfaces grew to become larger than the rate of dislocation multiplication, and therefore, plasticity could proceed only by the nucleation of new dislocations at considerably higher stresses [START_REF] Shan | Mechanical annealing and source-limited deformation in submicrometre-diameter Ni crystals[END_REF]. The implied 'dislocation starvation' was therefore linked to the fact the 'breeding' distance for dislocation multiplication becomes larger than the system size [START_REF] Rao | Athermal mechanisms of sizedependent crystal flow gleaned from three-dimensional discrete dislocation simulations[END_REF][START_REF] Weinberger | Surface-controlled dislocation multiplication in metal micropillars[END_REF][START_REF] Greer | Plasticity in small-sized metallic systems: Intrinsic versus extrinsic size effect[END_REF][START_REF] Bagheripoor | Length Scale Plasticity: A Review from the Perspective of Dislocation Nucleation[END_REF][START_REF] Shan | Mechanical annealing and source-limited deformation in submicrometre-diameter Ni crystals[END_REF][START_REF] Mordehai | Nanoindentation size effect in single-crystal nanoparticles and thin films: A comparative experimental and simulation study[END_REF][START_REF] Shan | Mechanical annealing and source-limited deformation in submicrometre-diameter Ni crystals[END_REF].

Following this logic, the 'brittle' response of pristine ultra-small samples after they reached the level of theoretical (ideal) strength, can be explained by massive homogeneous nucleation of dislocations [START_REF] Plans | Homogeneous nucleation of dislocations as bifurcations in a periodized discrete elasticity model[END_REF][START_REF] Miller | On the nonlocal nature of dislocation nucleation during nanoindentation[END_REF][START_REF] Garg | Universal scaling laws for homogeneous dislocation nucleation during nano-indentation[END_REF]. In other words, in such 'starved' samples yielding can be expected to proceed as a dislocation nucleation avalanche which can reach the size of the system [START_REF] Kiener | Source-controlled yield and hardening of Cu(100) studied by in situ transmission electron microscopy[END_REF][START_REF] Oh | In situ observation of dislocation nucleation and escape in a submicrometre aluminium single crystal[END_REF][START_REF] Benzerga | Micro-pillar plasticity: 2.5D mesoscopic simulations[END_REF][START_REF] Lee | Dislocation junctions and jogs in a free-standing FCC thin film[END_REF][START_REF] Ng | Effects of trapping dislocations within small crystals on their deformation behavior[END_REF]. The implied cooperative response becomes possible because the defect-free environment allows the nucleated dislocations to dynamically self-organize, taking full advantage of the un-screened long-range elastic interactions [START_REF] Miguel | Intermittent dislocation flow in viscoplastic deformation[END_REF][START_REF] Koslowski | Avalanches and scaling in plastic deformation[END_REF][START_REF] Uchic | Micro-compression testing of fcc metals: A selected overview of experiments and simulations[END_REF]. The emerging mechanical instability ultimately originates from strain softening whose origin will be revealed as a part of this study.

We note that sharp peak stresses at yielding have been also observed in some macroscopic samples where the large yield drop is usually followed by a Luders-like flow instability, e.g. [START_REF] Li | Flow stress behavior and deformation characteristics of Ti-3Al-5V-5Mo compressed at elevated temperatures[END_REF][START_REF] Ruan | Plastic yielding and tensile strength of near-micrometer grain size pure iron[END_REF]. In such cases the initial mobile dislocation density is usually low either because dislocations were annihilated during the annealing process or because they were blocked by solute atoms; in both cases the discontinuous yielding is associated with either collective depinning from Cottrell atmosphere or explosive nucleation from nano-sized grain boundary sources. A theoretical explanation of discontinuous yielding in the prototypical context of transformational plasticity was proposed in [START_REF] Truskinovsky | The origin of nucleation peak in transformational plasticity[END_REF] where the stress peaks were linked to the difference between (homogeneous) nucleation and propagation thresholds for internal instabilities.

The development of the computational tools accounting for the dislocational nature of plastic flows has become a priority because the inadequacy of the conventional methods based on CP were placing severe restrictions on the possibility to model plastically deforming ultra-small structural elements [START_REF] Zheng | Discrete plasticity in sub-10-nm-sized gold crystals[END_REF][START_REF] Zhang | Taming intermittent plasticity at small scales[END_REF][START_REF] Wang | Discrete shear band plasticity through dislocation activities in body-centered cubic tungsten nanowires[END_REF][START_REF] Parakh | Nucleation of Dislocations in 3.9 nm Nanocrystals at High Pressure[END_REF]. In particular, it was realized that to simulate discontinuous yielding in sub-micron crystals which involves massive nucleation of dislocations, the modeling approach cannot ignore lattice effects and must necessarily account for phenomena at the scale of dislocation cores. Capturing the attendant pattern formation is also challenging because the nucleated dislocations continue to interact strongly at many different scales [START_REF] Xia | Computational modelling of mesoscale dislocation patterning and plastic deformation of single crystals[END_REF][START_REF] Clouet | Dislocation locking versus easy glide in titanium and zirconium[END_REF][START_REF] Salvalaglio | Closing the gap between atomic-scale lattice deformations and continuum elasticity[END_REF][START_REF] Van Der Giessen | Roadmap on multiscale materials modeling[END_REF][START_REF] Bertin | Frontiers in the Simulation of Dislocations[END_REF].

Numerous computational alternatives to CP have been developed targeting different time and length-scales. They range from molecular dynamics (MD) [START_REF] Niiyama | Atomistic mechanisms of intermittent plasticity in metals: dislocation avalanches and defect cluster pinning[END_REF][START_REF] Zepeda-Ruiz | Probing the limits of metal plasticity with molecular dynamics simulations[END_REF] and similarly microscopic phase field crystal (PFC) method [START_REF] Salvalaglio | A coarse-grained phase-field crystal model of plastic motion[END_REF][START_REF] Chan | Plasticity and dislocation dynamics in a phase field crystal model[END_REF], to the more coarse grained phase field dislocation dynamics (PFDD) [START_REF] Hunter | Influence of the stacking fault energy surface on partial dislocations in fcc metals with a three-dimensional phase field dislocations dynamics model[END_REF][START_REF] Koslowski | A phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals[END_REF], discrete dislocation dynamics (DDD) [START_REF] Vattré | Modelling crystal plasticity by 3D dislocation dynamics and the finite element method: The Discrete-Continuous Model revisited[END_REF][START_REF] Ispánovity | Avalanches in 2D dislocation systems: plastic yielding is not depinning[END_REF][START_REF] El-Awady | Advances in Discrete Dislocation Dynamics Modeling of Size-Affected Plasticity[END_REF][START_REF] Bertin | GPU-accelerated dislocation dynamics using subcycling time-integration[END_REF], and continuum dislocation dynamics (CDD) [START_REF] Varadhan | Dislocation transport using an explicit Galerkin/leastsquares formulation[END_REF][START_REF] Hochrainer | Continuum dislocation dynamics: Towards a physical theory of crystal plasticity[END_REF][START_REF] El-Azab | Continuum Dislocation Dynamics: Classical Theory and Contemporary Models[END_REF]. The multi-scale quasi-continuum model (QC) attempts to bridge all the scales while accounting for each of them fully comprehensibly [START_REF] Shenoy | An adaptive finite element approach to atomic-scale mechanics-the quasicontinuum method[END_REF][START_REF] Miller | The Quasicontinuum Method: Overview, applications and current directions[END_REF].

MD approaches, including Molecular Statics (MS), and Density Functional Theory (DFT), accurately represent micro-mechanisms of plastic response while relying minimally on phenomenology [START_REF] Cia | Comprehensive Nuclear Materials[END_REF]. MD simulations have been particularly instrumental in the study of the homogeneous and heterogeneous dislocation nucleation [START_REF] Miller | On the nonlocal nature of dislocation nucleation during nanoindentation[END_REF][START_REF] Zepeda-Ruiz | Probing the limits of metal plasticity with molecular dynamics simulations[END_REF][START_REF] Parakh | Nucleation of Dislocations in 3.9 nm Nanocrystals at High Pressure[END_REF]. The only relative shortcoming of the atomistic simulations is that they are still computationally rather expensive in most applications, even at the small time and length scales of interest; also the problem of mapping to the macroscopic description in terms of the measurable quantities like stresses and strains is apparently not C. R. Physique, 0000, 1, n o 0, 000-000 yet fully resolved [START_REF] Zimmerman | Surface step effects on nanoindentation[END_REF][START_REF] Zimmerman | Deformation gradients for continuum mechanical analysis of atomistic simulations[END_REF][START_REF] Zepeda-Ruiz | Atomistic insights into metal hardening[END_REF][START_REF] Lim | Simulating dislocation plasticity in bcc metals by integrating fundamental concepts with macroscale models[END_REF]. Partial temporal averaging of atomistic molecular dynamics has emerged in the form a coarse-grained continuum PFC theory [START_REF] Elder | Modeling elasticity in crystal growth[END_REF]. The PFC approach was applied to the description of an ensemble of few interacting dislocation dipoles and was able to capture some cooperative features of crystal plasticity [START_REF] Chan | Plasticity and dislocation dynamics in a phase field crystal model[END_REF][START_REF] Skaugen | Separation of Elastic and Plastic Timescales in a Phase Field Crystal Model[END_REF][START_REF] Salvalaglio | A coarse-grained phase-field crystal model of plastic motion[END_REF]. However, such a detailed description of atomic lattices still remains rather demanding in terms of computer time, at least in the case of developed plastic flows with realistic number of interacting dislocations.

Discrete dislocation dynamics (DDD) approach was created to overcome the short time and length scales of atomistic methods. This model treats adequately long range interaction of dislocations without resolving the fine structure of the dislocation cores. To account for short range interactions, specific 'local rules' must be added governing, for instance, intersections and locks and some particular lattice scale effects have been indeed successfully modeled with the help of such additional phenomenological constructs, typically motivated by atomistic simulations [START_REF] Kubin | The modelling of dislocation patterns[END_REF][START_REF] Devincre | Three-Dimensional Simulations of Plastic Flow in Crystals[END_REF][START_REF] Cazacu | Multiscale Modeling of Heterogenous Materials: From Microstructure to Macro-Scale Properties[END_REF][START_REF] Po | Singularity-free dislocation dynamics with strain gradient elasticity[END_REF][START_REF] Wang | Reactions of lattice dislocations with grain boundaries in Mg: Implications on the micro scale from atomic-scale calculations[END_REF][START_REF] Geslin | Implementation of the nudged elastic band method in a dislocation dynamics formalism: Application to dislocation nucleation[END_REF][START_REF] Kohnert | Spectral discrete dislocation dynamics with anisotropic short range interactions[END_REF]. Being reinforced in this way, the DDD approach has emerged as an extremely useful approach for modeling evolution of many interacting dislocations [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF][START_REF] Dmitrieva | Investigation of the internal substructure of microbands in a deformed copper single crystal: experiments and dislocation dynamics simulation[END_REF]. However, it still remains challenging to account in the DDD framework for large deformations, crystal symmetry, lattice rotation, as well as the emergence of non-dislocational defects. Coarse graining of DDD, which opens the way towards modeling of dislocation patterning, has been attempted in the framework of the mesoscopic CDD methods where dislocational microstructures are modeled by continuum dislocational density fields [START_REF] Starkey | Theoretical development of continuum dislocation dynamics for finitedeformation crystal plasticity at the mesoscale[END_REF][START_REF] Acharya | Size effects and idealized dislocation microstructure at small scales: Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part I[END_REF][START_REF] Sandfeld | Pattern formation in a minimal model of continuum dislocation plasticity[END_REF]. While various phenomenological closure relations have been proposed to model the evolution of the dislocation density, the systematic development of this approach is hindered by the fact that rigorous statistical averaging in the ensemble of strongly interacting dynamic defects is still a big challenge.

In search of a micro-macro compromise a QC approach was proposed in [START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF] and then significantly developed in [START_REF] Shenoy | An adaptive finite element approach to atomic-scale mechanics-the quasicontinuum method[END_REF][START_REF] Miller | The Quasicontinuum Method: Overview, applications and current directions[END_REF][START_REF] Dobson | A multilattice quasicontinuum for phase transforming materials: Cascading Cauchy Born kinematics[END_REF][START_REF] Sorkin | A local quasicontinuum method for 3D multilattice crystalline materials: Application to shape-memory alloys[END_REF][START_REF] Kochmann | The Quasicontinuum Method: Theory and Applications[END_REF]. It is based on the observation that a fully atomistic resolution is necessary only in small spatial regions, while in most of the modeling domain the deformation fields can be represented by the classical continuum theory. A necessity of patching the continuum and discrete subdomains poses, however, a complex problem [START_REF] Miller | The Quasicontinuum Method: Overview, applications and current directions[END_REF][START_REF] Tadmor | Modeling Materials: Continuum, Atomistic and Multiscale Techniques[END_REF]. The QC method has been successfully used in many applications, including the study of nano-indentation, deformation of grain-boundaries and crack tip evolution [START_REF] Rodney | Structure and Strength of Dislocation Junctions: An Atomic Level Analysis[END_REF][START_REF] Knap | Effect of indenter-radius size on Au(001) nanoindentation[END_REF][START_REF] Yu | Edge dislocations interacting with a Σ11 symmetrical grain boundary in copper upon mixed loading: quasicontinuum method study[END_REF][START_REF] Jin | Quasicontinuum Simulation of the Effect of Lotus-Type Nanocavity on the Onset Plasticity of Single Crystal Al during Nanoindentation[END_REF].

Of particular interest to our study is the local version of the QC method representing a meso-scale compromise between continuum and atomistic approaches [START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF][START_REF] Sorkin | A local quasicontinuum method for 3D multilattice crystalline materials: Application to shape-memory alloys[END_REF]; a closely related approach is the interatomic potential FEM (IPFEM) [START_REF] Van Vliet | Quantifying the early stages of plasticity through nanoscale experiments and simulations[END_REF][START_REF] Zhu | Predictive modeling of nanoindentation-induced homogeneous dislocation nucleation in copper[END_REF]. In these effectively single scale methods the spatial domain is subdivided into discrete elements, as in a finite element method (FEM), and the deformation in each element is taken to be piece-wise affine. The energy density of an element, which then depends only on the displacement gradient at this point, is evaluated using interatomic potentials and the Cauchy-Born rule [START_REF]On the Cauchy-Born Rule[END_REF][START_REF] Weinan | Cauchy-born rule and the stability of crystalline solids: Static problems[END_REF][START_REF] Steinmann | Studies of validity of the Cauchy-Born rule by direct comparison of continuum and atomistic modelling[END_REF][START_REF] Podio-Guidugli | On (Andersen-)Parrinello-Rahman molecular dynamics, the related metadynamics, and the use of the Cauchy-born rule[END_REF]. These approaches were shown to capture adequately the basic structure of dislocation cores, even though they can misrepresent some truly atomistic features like deformations varying rapidly within the cut-off distance of the interatomic potential.

A mesoscale approach to crystal plasticity, which also allows one to treat dislocations in a fully continuum framework, and which is particularly relevant for our own development, is the phase field dislocation dynamics (PFDD) [START_REF] Koslowski | A phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals[END_REF][START_REF] Rodney | Phase field methods and dislocations[END_REF][START_REF] Beyerlein | Understanding dislocation mechanics at the mesoscale using phase field dislocation dynamics[END_REF]. This method has evolved from the original Landau approach to the modeling of phase transitions [START_REF] Chen | Phase-field models for microstructure evolution[END_REF][START_REF] Salman | Modeling of spatio-temporal dynamics and patterning mechanisms of martensites by phase-field and Lagrangian methods[END_REF][START_REF] Finel | Phase field methods: Microstructures, mechanical properties and complexity[END_REF][START_REF] Salman | The role of phase compatibility in martensite[END_REF][START_REF] Shchyglo | Martensitic phase transformations in Ni-Ti-based shape memory alloys: The Landau theory[END_REF][START_REF] Salman | Origin of stabilization of macrotwin boundaries in martensites[END_REF]. In PFDD lattice slips are described by scalar order parameters and the energy wells represent quantized lattice invariant shears [START_REF] Jin | Phase field microelasticity theory of dislocation dynamics in a polycrystal: Model and three-dimensional simulations[END_REF][START_REF] Rodney | Phase field methods and dislocations[END_REF][START_REF] Zheng | Improved phase field model of dislocation intersections[END_REF]; transition layers, separating regions with different amount of shear, represent the locations of dislocation lines. The Landau energy functional couples the tensorial linear elastic energy with the scalar lattice energy whose periodic structure is usually informed by atomic scale simulations based on the Cauchy-Born rule. The PFDD method enables simulations of much larger crystal sizes and much longer time scales than atomistic simulations [START_REF] Hu | Solute segregation and coherent nucleation and growth near a dislocation-a phase-field model integrating defect and phase microstructures[END_REF][START_REF] Chen | Phase-field models for microstructure evolution[END_REF][START_REF] Louchez | Microscopic Phase-Field modeling of hcp|fcc interfaces[END_REF][START_REF] Qiu | Predicting grain boundary structure and energy in BCC metals by integrated atomistic and phase-field modeling[END_REF]. While the extensions of PFDD to finite strains have recently appeared [START_REF] Biscari | Intermittency in crystal plasticity informed by lattice symmetry[END_REF][START_REF] Javanbakht | Phase field approach to dislocation evolution at large strains: Computational aspects[END_REF], the remaining challenge is that the structure of lattice invariant shears is resolved by scalar order parameters only approximately [START_REF] Xu | Phase-field-based calculations of the disregistry fields of static extended dislocations in FCC metals[END_REF][START_REF] Li | A review: applications of the phase field method in predicting microstructure and property evolution of irradiated nuclear materials[END_REF][START_REF] Salvalaglio | A coarse-grained phase-field crystal model of plastic motion[END_REF].

Mesoscopic tensorial model (MTM)

We have seen that despite many important advances in various specific problems, the sufficiently versatile computational method, allowing for natural coupling of different plastic 'mechanisms' while addressing realistic space and time scales2 , is still missing. The implied challenge is to build a synthetic mesoscopic tensorial approach dealing with large strains, while accounting correctly for both anisotropy and discreteness of the simplest Bravais lattices.

Different attempts along these lines can be found in the literature. The corresponding scalar approaches [START_REF] Landau | Application of a model of interacting atomic chains for the description of edge dislocations[END_REF][START_REF] Carpio | Edge dislocations in crystal structures considered as traveling waves in discrete models[END_REF][START_REF] Salman | Minimal integer automaton behind crystal plasticity[END_REF] can be viewed as generalizations of the minimalistic 1D Frenkel-Kontorova model [START_REF] Frenkel | On the theory of plastic deformation and twinning[END_REF][START_REF] Peierls | The size of a dislocation[END_REF][START_REF] Nabarro | Dislocations in a simple cubic lattice[END_REF] for the 2D case when only one slip system is activated. Despite their simplicity such models have been successful in describing dislocation cores [START_REF] Kovalev | Theoretical Description of the Crowdion in an Anisotropic Crystal Based on the Frenkel-Kontorova Model Including and Elastic Three-Dimensional Medium[END_REF], in simulate dislocation nucleation [START_REF] Plans | Homogeneous nucleation of dislocations as bifurcations in a periodized discrete elasticity model[END_REF][START_REF] Bonilla | Dislocations in cubic crystals described by discrete models[END_REF][START_REF] Lomdahl | Dislocation generation in the two-dimensional Frenkel-Kontorova model at high stresses[END_REF][START_REF] Srolovitz | Dislocation dynamics in the 2-d Frenkel-Kontorova model[END_REF] and even in modeling of plastic intermittency [START_REF] Salman | On the critical nature of plastic flow: One and two dimensional models[END_REF][START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF]. The tensorial models with linearized kinematics were proposed in [START_REF] Bulatov | A stochastic model for continuum elasto-plastic behavior. I. Numerical approach and strain localization[END_REF][START_REF] Minami | Nonlinear elasticity theory of dislocation formation and composition change in binary alloys in three dimensions[END_REF][START_REF] Onuki | Plastic flow in two-dimensional solids[END_REF][START_REF]Discrete models of dislocations and their motion in cubic crystals[END_REF]. They produced a more realistic picture of plastic flows but still could not account adequately for lattice-invariant shears associated with finite plastic slip [START_REF] Kaxiras | Energetics of large lattice strains: Application to silicon[END_REF].

The first attempts to incorporate the effects of geometrical nonlinearity were made in the context of a model of reconstructive phase transitions allowing for partial plastic accommodation [START_REF] Conti | A Variational Model for Reconstructive Phase Transformations in Crystals, and their Relation to Dislocations and Plasticity[END_REF]. A finite strain model focused directly on the modeling of plastic deformation was first developed for the case of highly anisotropic lattices of HCP type with a single slip system [START_REF] Salman | Minimal integer automaton behind crystal plasticity[END_REF]. Both of these papers can be viewed as different realizations of the original program of Ericksen who proposed that in nonlinear elasticity the energy periodicity should be made compatible with geometrically nonlinear kinematics of Bravais lattices [START_REF] Ericksen | Nonlinear elasticity of diatomic crystals[END_REF][START_REF]Loading Devices and Stability of Equilibrium[END_REF][START_REF]Special Topics in Elastostatics † †The research work herein reported was supported by a grant from the National Science Foundation[END_REF][START_REF]Some phase transitions in crystals[END_REF], see also the subsequent important developments of the mathematical formalism in [START_REF] Parry | On the elasticity of monatomic crystals[END_REF][START_REF] Folkins | Functions of two-dimensional Bravais lattices[END_REF][START_REF] Parry | Low-Dimensional Lattice Groups for the Continuum Mechanics of Phase Transitions in Crystals[END_REF][START_REF] Pitteri | Continuum models for phase transitions and twinning in crystals[END_REF].

Behind the coarse-grained approach of Ericksen was the general assumption that meso-scale material elements are exposed to an effective energy landscape which is globally periodic due to the presence of lattice invariant shears. From the perspective of such Landau-type continuum theory with an infinite number of equivalent energy wells, plastically deformed crystal emerges as a multiphase mixture of equivalent 'phases'. Plastic yield can be then interpreted as an escape from the reference energy well and plastic 'mechanisms' can be linked to low-barrier valleys in the energy landscape. Friction type dissipation controlling dynamics in continuum crystal plasticity emerges in such theory as a result of a homogenized description of an overdamped athermal dynamics in a rugged energy landscape [START_REF] Puglisi | Thermodynamics of rate-independent plasticity[END_REF][START_REF] Mielke | From Discrete Visco-Elasticity to Continuum Rate-Independent Plasticity: Rigorous Results[END_REF].

The 2D scalar model of crystal plasticity presented in [START_REF] Salman | Minimal integer automaton behind crystal plasticity[END_REF][START_REF] Salman | On the critical nature of plastic flow: One and two dimensional models[END_REF] was effectively a single slip system version of such theory. Building upon [START_REF] Conti | A Variational Model for Reconstructive Phase Transformations in Crystals, and their Relation to Dislocations and Plasticity[END_REF], its 2D tensorial version [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF] could already incorporate adequately the full GL(2, Z) symmetry of the Bravais lattice; the corresponding discrete symmetry group can be represented by 2 × 2 invertible matrices with integer entries and determinant ±1 [START_REF] Ericksen | Nonlinear elasticity of diatomic crystals[END_REF][START_REF]Some phase transitions in crystals[END_REF]. In such kinematically nonlinear theory the role of the order parameter is played by the metric tensor (characterizing local deformation), and the bottoms of the energy wells correspond to lattice invariant deformations.

Since the ground state in the continuum theory of this type is necessarily degenerate (hydrostatic) [START_REF]Loading Devices and Stability of Equilibrium[END_REF][START_REF] Fonseca | Variational methods for elastic crystals[END_REF], the regularization is necessary, as in the case of PFDD or any other Landau-type theory. In [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF] such regularizing internal length scale was associated with the size of the mesh generating discrete elastic elements which define the meso-scale. Given the large magnitude of the 'transformation strain', different 'phases' end up being localized at the scale of such mesoscopic elements and the domain boundaries appear macroscopically as linear defects mimicking dislocations.

In what follows we refer to the approach proposed in [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF] as the mesoscopic tensorial model (MTM). The main advantage of this approach is that it is formulated in terms of macroscopically measurable quantities (stress and strain) while being able to distinguishing between different symmetry induced configurations of dislocation cores. It can therefore account adequately for both long-and short-range interactions between dislocations. Most importantly, it allows for topological transitions associated with dislocation nucleation and annihilation even though the corresponding 'reactions ' appear as blurred on the scale of regularization. Last but not least, in the MTM approach the interaction of dislocations with various obstacles, including pinning by impurities and depinning from nucleation sources, can be handled without introducing ad-hoc relations.

In this paper the MTM is used to conceptualize the brittle-like response of nominally ductile submicron crystals. Through a series of numerical experiments we uncover the origin of the correlated dislocational response of such crystals with different symmetry at loading levels approaching the theoretical (ideal) strength. In particular, we reveal the mechanism of the observed orientation dependence of the degree of 'brittleness' and rationalize the system size dependence of the associated collective nucleation event.

GL(2, Z) energy density

For simplicity, we deal in this paper with the simplest two dimensional Bravais lattices: square (with lower, square symmetry) and triangular (with higher, hexagonal symmetry). In addition to its theoretical interest, the study of plasticity in 2D crystals is technologically relevant because such crystals have been recently found important in many applications [START_REF] Bhimanapati | Recent Advances in Two-Dimensional Materials beyond Graphene[END_REF][START_REF] Chen | Two Dimensional Ice from First Principles: Structures and Phase Transitions[END_REF][START_REF] Hoang | Formation of Two-Dimensional Crystals with Square Lattice Structure from the Liquid State[END_REF][START_REF] Zhang | Two dimensional hexagonal boron nitride (2D-hBN): synthesis, properties and applications[END_REF][START_REF] Akinwande | A review on mechanics and mechanical properties of 2D materials-Graphene and beyond[END_REF][START_REF] Chen | Two-Dimensional Metal Nanomaterials: Synthesis, Properties, and Applications[END_REF][START_REF] Kryuchkov | Complex crystalline structures in a two-dimensional core-softened system[END_REF][START_REF] Van Hoang | Compression-induced square-triangle solid-solid phase transition in 2D simple monatomic system[END_REF][START_REF] Ma | Atomic imaging of the edge structure and growth of a two-dimensional hexagonal ice[END_REF].

The main ingredient of the MTM approach is the objective elastic energy density which respects the GL(2, Z) symmetry. To construct such energy, we need to first introduce the deformation of a continuum body y = y(x), where y and x are the position vectors in the current (Eulerian) configuration and the reference (Lagrangian) configuration, respectively. Denote by F the deformation gradient ∇y = ∂y/∂x. The frame indifferent elastic energy density Φ(C) must depend on F through the metric tensor C = F T F. To account adequately for all deformations that map a Bravais lattice onto itself, we require that the strain energy density satisfies Φ(C) = Φ(M T CM) for any matrix M in GL (2, Z). This invariance follows from the fact that the same lattice can be generated by the two sets of lattice vectors ēA and e A if and only if ēA = M B A e B . In the presence of such symmetry, the space of metric tensors C K L = e K .e L partitions into periodicity domains. Therefore, when the energy is known in one of the domains, it can be automatically extended to all other equivalent domains by the use of global tensorial periodicity.

The space of metric tensors in the 2D case is three dimensional and in Fig. 1 we show a two dimensional section of this space defined by the condition det C = 1. For better visibility equivalent tensorial periodicity subdomains are shown in two alternating colors, gray and white. Each of these subdomains contains necessarily one and only one copy of a particular Bravais lattice as shown in Fig. ?? fig:detc1. In this figure, equivalent zero-parametric families of square lattice configurations are illustrated with solid squares S i , S ij , ...S ī, S ī j, ... while solid triangles T i , T ij , ...T ī, T ī j, ... corresponding to similar families of equivalent triangular lattices. Red lines correspond to one-parametric families of rectangular lattices and blue lines to similar families of rhombic lattices. Finally, the open 2D domains correspond to two-parametric families of equivalent oblique (monoclinic) lattices, see also [START_REF] Cayron | The transformation matrices (distortion, orientation, correspondence), their continuous forms and their variants[END_REF][START_REF] Gao | A Cayley graph description of the symmetry breaking associated with deformation and structural phase transitions in metallic materials[END_REF][START_REF] Gao | Determination of twinning path from broken symmetry: A revisit to deformation twinning in bcc metals[END_REF][START_REF] Gao | Phase Transformation Graph and Transformation Pathway Engineering for Shape Memory Alloys[END_REF]. Note that in this representation we do not see the 'volumetric' part of the metric tensor responsible, for instance, for void formation [START_REF] Marconi | Diffuse interface approach to brittle fracture[END_REF]. The darkened area in Fig. 1, denoted by D, is known as the 'fundamental' domain [START_REF] Engel | Geometric Crystallography: An Axiomatic Introduction to Crystallography[END_REF]. Among all the equivalent metrics describing symmetry-related lattices, this domain contains the "minimal" ones, known as 'reduced forms of Lagrange'. An arbitrary basis e i producing an arbitrary metric C has its unique 'Lagrange reduced' copy ẽi with the corresponding metric C in D.

The Lagrange reduction is obtained using the following rule: among all symmetry related copies of the original basis, the reduced one is formed by the shortest non-colinear vectors with an acute angle between them, see for instance [START_REF] Conti | A Variational Model for Reconstructive Phase Transformations in Crystals, and their Relation to Dislocations and Plasticity[END_REF]. Therefore, the domain D can be defined as:

D = 0 < C 11 ≤ C 22 , 0 < C 12 ≤ C 11 2 .
(

) 1
Inside this domain, the 'Lagrange reduced' square lattice, with basis vectors aligned with the close-packed directions e 1 = (1, 0), e 2 = (0, 1), corresponds to point S 0 in Fig. 1. Instead, the 'Lagrange reduced' triangular lattice with its basis vectors h 1 = γ(1, 0) and h 2 = γ(1/2, 3/2)), where γ = (4/3) 1/4 , corresponds to the point T 0 in Fig. 1.

In view of the GL(2, Z) periodicity, it is sufficient to know the elastic energy for C ∈ D. One can then use the symmetry to write Φ(C) = Φ 0 ( C), where C = M T CM is the Lagrange reduced metric, and the function Φ 0 is defined only inside D. The condition that the function Φ 0 is twice-differentiabile can be satisfied if the deviatoric part of Φ 0 is a sixth order polynomial [START_REF] Parry | On the elasticity of monatomic crystals[END_REF][START_REF] Parry | Low-Dimensional Lattice Groups for the Continuum Mechanics of Phase Transitions in Crystals[END_REF]. With hexagonal symmetry chosen as global, one can write the minimal expression for such polynomial [START_REF] Conti | A Variational Model for Reconstructive Phase Transformations in Crystals, and their Relation to Dislocations and Plasticity[END_REF]

Φ 0 d ( C/(det C) 1/2 ) = βψ 1 ( C/(det C) 1/2 ) + ψ 2 ( C/(det C) 1/2 ), (2)
where,

ψ 1 = I 1 4 I 2 -41 I 2 3 /99 + 7 I 1 I 2 I 3 /66 + I 3 2 /1056, and ψ 2 = 4 I 2 3 /11 + I 1 3 I 3 -8 I 1 I 2 I 3 /11 + 17 I 3
2 /528. and we used the following hexagonal invariants 2 , and 3 . The elastic energy density (2) depends on a single parameter β that may be used to enforce a particular symmetry on the reference state. For instance, the choice β = -1/4 guarantees that the global energy minimizers correspond to square lattices while the choice β = 4 shifts the bias towards the triangular lattice. In what follows, the total elastic energy will be then taken in the form plays only minor role in this study by affecting the structure of dislocation cores and controlling the formation of voids, will be chosen in the form Φ v (s) = µ(slog(s)). This choice, instead of a more realistic Lennard-Jones type potential, will allow is to avoid various volumetric instabilities while still excluding infinite compression. Moreover, to ensure that the strain field remains close to the surface det C = 1 in our numerical experiments we choose large value of the bulk modulus by setting

I 1 = 1 3 ( C11 + C22 -C12 ), I 2 = 1 4 ( C11 - C22 ) 2 + 1 12 ( C11 + C22 -4 C12 )
I 3 = ( C11 -C22 ) 2 ( C11 + C22 -4 C12 ) -1 9 ( C11 + C22 -4 C12 )
Φ 0 ( C) = Φ 0 d ( C/(det C) 1/2 ) + Φ v (det C).
µ = 25.
We remark that the particular potential (2) was chosen for illustrative purposes only and for most applications in crystal plasticity the piece-wise quadratic potentials with GL(2, Z) symmetry would be sufficient, see for instance [START_REF] Salman | On the critical nature of plastic flow: One and two dimensional models[END_REF][START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF]. Outside crystal plasticity, the general approach of MTM can be used as well with potentials satisfying less restrictive symmetry constraints, see for instance, [START_REF] Friesecke | Validity and Failure of the Cauchy-Born Hypothesis in a Two-Dimensional Mass-Spring Lattice[END_REF].

Energy landscape

To visualize the local energy landscape we use the plane with coordinates ((C 11 -C 22 )/ 2,C 12 ). In Fig. 2(a) we illustrate the potential (2) with β = -1/4 describing an energy landscape whose absolute minima are the equivalent square lattices. To model plastic flows in triangular lattices with higher, hexagonal symmetry, we choose in (2) the parameter value β = 4 and the resulting elastic potential is illustrated in Fig. 3(a).

In Fig. 2(a) and Fig. 3(a) one can clearly distinguish between the 'soft' directions located inside the deep valleys of the energy landscape (marked in blue), and the 'hard' directions pointing away from the valleys and forcing the system to climb high energy barriers separating the valleys (yellow and red regions). Inside each low energy valley the energy is periodic with the bottoms of the energy wells corresponding to equivalent lattices. Thus, in the case of square lattices (Fig. 2(a) the equivalent minima are denoted by S i , S ij , ...S ī, S ī j, ... and in the case of triangular lattices (Fig. 3(a)), by T i , T ij , ...T ī, T ī j, ... Note that while the energy is locally convex near the bottom of each energy well and can be approximated by the standard Hookean paraboloid, the global energy landscape is highly nonconvex and can be rather characterized as rugged. In particular this means that outside the immediate vicinity of the energy bottoms the system finds itself in spinodal regions where the constitutive behavior softens opening the possibility for various mechanical instabilities. Note that the local maxima of the one parametric energy landscape along the low energy valleys turn out to be located very close to the configurations corresponding to saddles describing perfect lattices with alternative symmetries: triangular lattices, in the case when we shear a square lattice (see triangles near the local maxima along the white dashed paths in Fig. 2(a)) and square lattices, in the case when we shear a triangular lattice (see squares near the local maxima along the white dashed paths in Fig. 3(a)).

CP theory obviously appreciates the implied complexity of the energy landscape but not in the same way as MTM. Thus, the low energy valleys are modeled as zero energy plastic 'mechanisms'. The periodicity of the energy inside the valleys is neglected but instead the phenomenological friction type lays are introduced to govern the 'working' of such mechanisms. Elasticity is usually assumed to be linear and is only accounted for around the reference energy well, which is then extended globally. This does not create problems in CP because at large strains plastic deformation along the valleys always wins against elasticity which becomes at such strains prohibitively expensive energetically. To describe in the CP framework the coupling between different plastic 'mechanisms, involved for instance in latent hardening [START_REF] Ortiz | Nonconvex energy minimization and dislocation structures in ductile single crystals[END_REF], additional phenomenological relations have to be formulated. In MTM those will be automatically accounted for due to the presence of various saddle points (mountain paths) that are apparent in Fig. 2(a) and Fig. 3(a).

To illustrate a more detailed structure of the energy landscape along the particular low energy valleys we now consider two parametric families of square lattices transformed to one another by the volume preserving shear deformations of the form

F(α, φ) = 1 + αa(φ) ⊗ a ⊥ (φ), (3)
where α is the amplitude of shear, a(φ) = R(φ)e 1 , R(φ) is a matrix describing a counterclockwise rotation by the angle φ, e 1 is a unit vector of the Cartesian coordinate system directed along the close-packed x direction and a ⊥ is a vector orthogonal to a.

For instance, if we fix one of the parameters at φ = 0 and consider a one parametric family involving the metric tensors C(α), we obtain the path 1-1, shown in Fig. 2(a) by a dashed white line facing left. It starts at the reference state S 0 , corresponding to the unstressed square lattice, and at integer values of parameter α generates an infinite sequence of equivalent unstressed replicas of the same square lattice. For instance, the 'closest' lattices configurations S 1 and S 1 can be reached from S 0 by the elementary shear deformations F(±1, 0).

If we choose instead φ = π/2, we obtain another one parametric family of shears described by the path 2-2, shown in Fig. 2(a) by a dashed white line facing right. It also starts at the reference state S 0 and produces at integer values of parameter α an infinite sequence of equivalent unstressed configurations of the square lattice. Once again, the 'closest' lattices S 2 and S 2 can be reached reached through the shear deformations F(±1, π/2).

The energy landscape along these symmetric paths is the same, see Fig. 2(b). Note that around the symmetric energy minima α = 0 and α = ±1, visible in such a graph, the energy is quadratic but then it loses convexity and the corresponding mechanical response starts to soften. Such softening is ultimately behind the break up of elastic response and the emergence of yield. If the system is quasi-statically driven through such an energy landscape, it undergoes a succession of snap-back instabilities which, under overdamped dynamics, merge into the dissipative plastic response which postulated phenomenologically in the classical CP approach, see [START_REF] Puglisi | Thermodynamics of rate-independent plasticity[END_REF] for details.

Note that a composition of shear deformations from the paths 1-1 and 2-2 can be interpreted as an activation of a double-slip. For instance, the mapping F(1, π/2) F(1, 0) brings the system from the reference lattice configuration S 0 to the equivalent lattice configuration S 12 after it slips in two perpendicular directions. In fact, all the bottoms of the energy wells S ij , ...S ī j, ..., corresponding to equivalent square lattices, can be reached similarly by combining quantized simple shears (3) with φ = 0 or φ = π/2. Moreover, as we have already mentioned, a simple shear path leading away from one of the energy wells passes near the saddles corresponding to ideal lattices with different symmetries. It is then natural to conclude that these saddles serve as a switching points activating double slip and inciting other composite shears. Based on these observations, one can argue that multiple alternative lattice symmetries, that are virtually invisible in the 'single symmetry' classical CP approach, may be contributing fundamentally to the complexity of plastic flows in crystals. Periodicity of the energy along the particular tensorial directions, illustrated in in Fig. 2(b), is usually directly postulated in PFDD, with each of these directions bringing in a scalar order parameter of its own. To describe a composition of the corresponding shear deformations, the PFDD order parameters need to be coupled phenomenologically. Instead, in the MTM, such coupling is automatic due to the global symmetry requirements. Moreover the MTM approach allows one to assess the height and the structure of the barriers separating the corresponding valleys and controlling the activation of double-slip. Note also, that a combination of scalars can represent a globally defined tensorial variable (metric tensor) only approximately.

(C 11 -C 22 )/ 2 C 12 (a) 1 1 3 3 2 2 T 0 T 1 T 2 T 1 T1 ¯1 ¯T¯2 1 Φ Figure 3.
Turning back to Fig. 2(a), we observe that the two other symmetric paths 3-3 and 4-4, shown by dashed blue lines, correspond to simple shears applied to the rotated crystal with either φ = π/4 or φ = 3π/4. While these symmetric shears are no longer aligned with the dense planes of our square lattice, the corresponding loading paths are still special since they bring the system from the reference energy well S 0 directly to the second nearest energy wells, represented in Fig. 2(a) by the small squares S 11 , S 22 , S 11 and S 22 . Along such paths the energy landscape is also periodic but the period doubles and the separating barriers become much higher as these paths are largely located outside the low energy valleys. While in CP such paths are not viewed as plastic 'mechanisms' that can be activated, they are still accounted for in the MTM which deals with the whole tensorial landscape and therefore encompasses all possible paths.

Finally, the red dashed line in Fig. 2(a) illustrates a non-symmetric loading path with φ = arctan(1/2), also originating in S 0 . This path exhibits the characteristic loading-unloading asym-metry of generic simple shear deformations. In particular, along this path the equivalent square energy well S 2 can be reached through 'unloading' at F(-1, arctan(1/2)) but none of the nearest square wells are reached upon 'loading' with α positive.

The analysis above suggests the strong dependence of the plastic response on the orientation of the crystalline sample inside a hard loading device, even if it always applies the same simple shear deformation. The fact that there are soft and hard orientations has been long known in crystal plasticity. MTM approach presents a more nuanced picture, showing that the soft paths may come not only with different size of the periodically placed barriers but also with different periodicity, while along the hard paths one may expect extremely high barriers and generic loading-unloading asymmetry. The crucial advantage of the MTM is that all this complexity does not have to be postulated and emerges instead from the global symmetry of the energy landscape.

The obtained general picture is confirmed if we turn from square to triangular lattices, however, as the analysis below shows, the important differences between the energy landscape in lower symmetry (square) and higher symmetry (triangular) lattices exists.

In the case of triangular lattices we can again illustrate the available plastic 'mechanisms' by considering again the one parametric simple shear paths. The single slip plasticity in triangular lattices can be modeled by the shears

F(α, φ) = 1 + αb 1 (φ) ⊗ b 2 (φ), (4)
where b A (φ) = R(φ)h A with A = 1, 2. The (non-unit) lattice vectors h A have been introduced earlier and the vectors b A are their duals defined by the conditions b A .b B = δ B A . Using (4), we can generate three symmetric shear paths originating in the unstressed lattice T 0 and describing plastic slips along the three close-packed directions. They correspond to crystal orientations φ = 0, π/3, 2π/3, see the white dashed lines in Fig. 3(a). Here again, for the integer values of α we obtain infinite families of equivalent replicas of the original triangular lattice. In particular, the lattice configurations T 1 and T 1, located on the path 1-1, can be reached by the shear deformations F(±1, 0) starting from the lattice configuration T 0 (see the dashed white line facing left) Upon further 'unloading' F(-2, 0), another triangular lattice T 11 can be reached, while the lattice configuration T 11 is reachable by the matching additional 'loading' F(2, 0) (not visible in Fig. 3(a)).

The symmetry related paths 2-2 and 3-3, also shown in Fig. 3(a) by dashed white lines, correspond to the shear deformations F(α, π/3) and F(α, 2π/3). Both of them also originate in the point T 0 and produce at integer values of α infinite families of equivalent replicas of the original triangular lattice. Thus, along the path 2-2, one can reach the equivalent lattices T 2 , T 1, and T 12 , all visible in Fig. 3(a). Similarly, the equivalent lattice configurations T 1 and T 2 are reachable along the path 3-3.

For each of these paths, which all correspond to the known plastic 'mechanisms', the one dimensional periodic energy landscape is exactly the same, see Fig. 3(b). The maxima of such energy landscapes are located close to the saddle points of the global energy landscape corresponding to square lattice configurations. Because they follow the low energy valleys of the global energy landscape (shown in blue in Fig. 3(a)), these simple shear paths correspond to 'soft' loading directions for a triangular lattice. Observe that such valleys are splitting into two in each of the energy wells which is clearly an easy mechanism for the activation of multi-slip in such higher symmetry lattices.

A non generic path with φ = π/2 is also shown in Fig. 3(a) (see the blue dashed line). In this case, the crystal is not driven along one of the low energy valleys, however it does pass through the bottoms of the energy wells (with period 3!). If we load the crystal in shear along this 'hard' direction, we can expected to reach the 'yellow' and even 'red' zones in Fig. 3(a) and therefore acquire considerable elastic energy before the ultimate breakdown of an elastic state. Curiously enough, as we show below, this does not happen.

To conclude this Section, we remark that, that since we used a particular expression for the energy (2), we did not deal with the most general structure of the energy landscape for square and tetragonal lattices. However, the discussed energy landscape was sufficiently general in the sense that the location of the energy wells and therefore the configuration of the energy valleys and the connecting saddle points were dictated exclusively by the GL(2, Z) symmetry. The main potentialspecific features of the energy landscape are the curvature of the wells near their bottoms and the height of the energy barriers. Such quantitative ambiguity is the element which MTM shares with all other Landau type theories.

Note also that the level of smoothness of the Landau potential can be improved if we replace the six order polynomials used in (2), by the higher order polynomials. Various exponential approximations were discussed in [START_REF] Folkins | Functions of two-dimensional Bravais lattices[END_REF] and an example of an analytic function with GL(2, Z) symmetry was presented in [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF]. However, in general, such potentials are analyzed numerically and the possibility of a relatively simple analytic representation is not an issue. Therefore, physically relevant potentials can be constructed directly inside the fundamental domain D by analyzing homogeneous deformations of atoms and applying the Cauchy-Born rule. To this end one can use interatomic potentials, embedded atom methods and even quantum mechanics. Such numerically constructed energy potentials can be then continued by symmetry using the Lagrange reduction mentioned above or other similar approaches [START_REF] Engel | Geometric Crystallography: An Axiomatic Introduction to Crystallography[END_REF][START_REF] Grosse-Kunstleve | Numerically stable algorithms for the computation of reduced unit cells[END_REF][START_REF] Andrews | A space for lattice representation and clustering[END_REF].

Computational approach

We now turn to the numerical implementation of the MTM approach. Its underlying idea is that rateindependent crystal plasticity can be modeled within the framework of the appropriately discretized nonlinear elasticity combined with athermal, overdamped dynamics. If the loading is sufficiently slow, the exact nature of such dynamics is not essential as long as it effectively performs incremental energy minimization [START_REF] Puglisi | Thermodynamics of rate-independent plasticity[END_REF][START_REF] Mielke | From Discrete Visco-Elasticity to Continuum Rate-Independent Plasticity: Rigorous Results[END_REF]. As the loading evolves the system remains on the same branch of local minima of the energy till the latter ceases to exist and then it switches to another equilibrium branch. Such switching is a dynamic process (avalanche) which can be considered instantaneous at the time scale of the loading but it still contributes to dissipation. In this representation, plastic flow emerges as a set of equilibrium solutions of nonlinear elasticity with dissipation taking place exclusively during the abrupt branch switching events [START_REF] Puglisi | Thermodynamics of rate-independent plasticity[END_REF].

We recall that solution of an elastic problem implies local minimization of the energy

W = Ω Φ(∇y)d x,
which is prescribed on a reference domain Ω with unit volume. We assume that the system is loaded by an affine displacement field prescribed on ∂Ω (hard device). The conditions of mechanical equilibrium are usually formulated in terms of the first Piola-Kirchhoff stress tensor P = ∂Φ/∂F. In the index free form they can be written as

∇ • P = 0.
Using the Eulerian i , j = 1, 2 and the Lagrangian K , L = 1, 2 indexes and assuming summation on repeated indexes, we can rewrite these equations in the form

A i K j L y j ,K L = 0,
where

A i K j L = ∂ 2 Φ 0 ( C) ∂F i K ∂F j L . (5)
Here C = M T CM, where the integer-valued matrix M with determinant equal to one is computed for each value of C using the Lagrange reduction algorithm. The matrix M depends on the current state of deformation F in a piece-wise constant manner and therefore can be considered as constant in [START_REF] Andrews | A space for lattice representation and clustering[END_REF] which can be then rewritten as

A i K j L = M M N M PQ ∂C M P ∂F j L ∂C AB ∂F i K M AC M B D H MQC D + ∂ 2 C M P ∂F i K ∂F j L Σ NQ . (6)
Here we introduced the tensors

Σ Σ Σ =   ∂Φ 0 ∂ C11 1 2 ∂Φ 0 ∂ C12 1 2 ∂Φ 0 ∂ C12 ∂Φ 0 ∂ C22   , (7)
and

H =           ∂ 2 Φ 0 ∂ C 2 11 1 2 ∂Φ 0 ∂ C12∂ C11 1 2 ∂Φ 0 ∂ C12∂ C11 ∂ 2 Φ 0 ∂ C22∂ C11   1 2   ∂Φ 0 ∂ C11∂ C12 1 2 ∂Φ 0 ∂ C 2 12 1 2 ∂Φ 0 ∂ C 2 12 ∂Φ 0 ∂ C22∂ C12   1 2   ∂Φ 0 ∂ C11∂ C12 1 2 ∂Φ 0 ∂ C 2 12 1 2 ∂Φ 0 ∂ C 2 12 ∂Φ 0 ∂ C22∂ C12     ∂ 2 Φ 0 ∂ C11∂ C22 1 2 ∂Φ 0 ∂ C12∂ C11 1 2 ∂Φ 0 ∂ C12∂ C11 ∂ 2 Φ 0 ∂ C 2 22          
, which are assumed as known inside the fundamental domain D. We also recall here for convenience that ∂C K L /∂F i j = δ K j F i L + δ L j F i K , and

∂ 2 C K L /∂F i M ∂F j N = (δ K M δ LN + δ K N δ LM )δ i j .
As we have already mentioned, the continuum elasticity problem, formulated above, is highly degenerate which is the property the MTM shares with other similar Landau type theories. Usually such theories are regularized through the introduction of an internal length scale. In contrast to the conventional Ginzburg-Landau approaches like PFDD, relying for regularization on higher gradients of the order parameters, in MTM the regularization is achieved by spatial discretization. More precisely, deformation is assumed to be piecewise linear and the elastic response is attributed to discrete material elements whose size is viewed as a physical parameter of the model [START_REF] Conti | A Variational Model for Reconstructive Phase Transformations in Crystals, and their Relation to Dislocations and Plasticity[END_REF][START_REF] Salman | Minimal integer automaton behind crystal plasticity[END_REF][START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF].

We therefore need to reduce the space of admissible deformations to compatible piece-wiseaffine mappings. To this end we build a network whose discrete nodes are labeled by integer valued coordinates a = 1, ..., N 2 . We assume that each element of the network is a deformable triangle and write the displacement field in the form u(x) = u a N a (x), where N a (x) are the compactly supported shape functions, u a are the amplitudes of nodal displacements and summation over repeated indexes effectively extends over elements containing or bounding point x. The mesoscopic deformation gradient is then F(x) = 1 + ∇u(x), and the equilibrium equations can be written in the form

∂W ∂u a = Ω P(F)∇N a (x)d x = 0, (8)
where P = 2FMΣ Σ ΣM T . This problem can be solved by quasi-Newton method followed by the so called NR 'refinement' when the initial guess is too far from the solution for Newton-Raphson method to converge initially [START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF].

More specifically, to solve [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF] for u a we first use the L-BFGS algorithm [START_REF] Bochkanov | Alglib", Available from: www alglib net[END_REF] which builds a positive definite linear approximation of (8) allowing one to make a quasi-Newton step lowering W . Such iterations continue till the increment in total energy W becomes sufficiently small. The obtained approximate solution is then used as an initial guess w a to solve, using LU factorization [START_REF] Sanderson | Armadillo: a template-based C++ library for linear algebra[END_REF][START_REF] Fishman | The ITensor Software Library for Tensor Network Calculations[END_REF], the equations for the correction dw a

K ab i j d w b j + R a i = 0, (9)
where

K ab i j = A i K j L (F) ∂N a ∂x K ∂N b ∂x L , R a i = P i K (F) ∂N a ∂x K . ( 10
)
The displacement field can be updated in this way till the value of the integral in ( 8) is sufficiently small and then the loading parameter can be advanced again. In the case when the applied deformation is, a homogeneous shear F(α, φ) with fixed orientation φ, the loading parameter is the shear amplitude α. By changing this parameter in increments of 10 -4 , we advance the displacement field u(α, φ) = ( F(α, φ) -1)x for all nodes a on the boundary of the body ∂Ω.

The outcome of such numerical experiments depends on the value of the internal length scale entering the discrete problem through the finite element size h. Each element contains n 2 interacting atoms, where n = h/a and a is a fixed atomic scale. For such element deformed in simple shear we can use one of the known ab initio methods to compute its energy, which becomes close to periodic as we increase the value of h. The value of this internal length scale used in MTM depends on the smallest range of energy periodicity required in the problem; in many crystals the periodicity at the level of the few first energy wells can be captured already for h ∼ 10a.

Since the linear size of the macroscopic sample is L = 1, the small dimensionless parameter of the MTM is h/L = 1/N where N 2 is the number of the (mesoscopic) elements. Thus, if h is in nm size range, the simulations with N ∼ 10 3 would describe a micrometer size crystal. When N is small, dislocation cores emerge as blurred because the scales smaller than 1/N are homogenized out. While such cut-offs may compromise the short-range interaction of dislocations, long-range interactions at distances larger than 1/N will be still captured correctly.

Numerical experiments

To understand the paradoxical crack-free brittle behavior of very small, initially dislocation-free crystals, we used MTM and conducted a series of numerical experiments with square and triangular crystals. The samples were subjected to simple shear, and we performed experiments with several different orientations of the same samples in the same hard loading device. We used square cut samples in all our experiments and could vary the system size up to N = 1024.

Macroscopic response

The obtained macroscopic responses for crystals described by the piece-wise smooth potential (3) are summarized in Fig. 4(a) and Fig. 4(b). We show separately the energy-strain and the stressstrain responses for square and triangular crystals loaded in different orientations. The energy was obtained by incremental local minimization starting from the reference state. To compute the stressstrain curve, we define the resolved-shear stress in the direction of loading as τ(α) = dW/dα = Ω P :

(d F/dα)d x.
In each of our numerical experiments, the initial phase of the deformation history is (nonlinear) elastic. During this stage, all elements have associated metrics that remain inside the extended fundamental neighborhood of the initial phase (Pitteri neighborhood [START_REF] Pitteri | Continuum models for phase transitions and twinning in crystals[END_REF]), comprised of four (in case of square symmetry) or six (in case of triangular symmetry) symmetric replicas of the fundamental domains. While the response of the crystal inside such domain is nominally elastic, our experiments show that the affine deformation becomes unstable before the boundary of this neighborhood is reached. The instability takes place at a critical value of the loading parameter α = α * c which depends not only on crystal symmetry and crystal orientation, see Fig. 4, but also on the value of N (see below). The breakdown of affine elastic regime takes the form of a catastrophic drop in both stress and energy. While these drops are still associated with plastic deformation, they are highly reminiscent of brittle fracture.

We now discuss how the macroscopic mechanical behavior during such discontinuous yield depends on crystal symmetry and how it varies when crystals of the same symmetry are differently orientated in the loading machine.

We start with a pristine square crystal that we load along the principal slip direction with φ = 0. In such an 'easy' glide, the instability is preceded by the purely elastic softening, and the yield takes place near the maximal load. However, instead of conventional continuous yielding, the crystal experiences discontinuous yield as it abruptly loses almost all accumulated energy with stress dropping almost to zero. This means that the crystal manages to expel almost all nucleated dislocations away from its bulk, either by annihilating or sending them towards the boundaries. If the boundaries were unconstrained, the transition would be pristine-to-pristine, but in the hard-device loading conditions, at least some of the nucleated dislocations end up forming low energy piles-up near the boundaries (see more about this in the next Section).

α = 0 = arctan( 1 2 ) = π/4 ϕ ϕ ϕ = 0 = arctan( 1 2 ) = π/4 ϕ ϕ ϕ W (c) α (d) = 0 = π/3 = π/2 ϕ ϕ ϕ 0 0 .1 0.2 0 .3 0.4 0
A similar almost pristine-to-pristine discontinuous-yield occurs in the case of a less symmetric but still non-generic shear with φ = π/4 but only after both the elastic energy and the stress reach a much higher value. As we have already seen, along this loading path, the periodicity strain doubles, and the barriers become much higher than in the case of shear loading directed along the low energy valley. Instead, the crystal is first driven away from such valley and, before yielding, manages to accumulate a considerable amount of energy without an apparent softening. It is then rather remarkable that the eventual breakdown of the elastic state leads to almost complete relaxation of considerable elastic stress. The explanation probably lies in the fact that the crystal's orientation is still special vis a vis the applied load.

Loading of the crystal in a non-generic tensorial direction φ = arctan(1/2) leads to a moderate softening in the elastic range, which again ends with a discontinuous yield. The resulting dislocation avalanche relaxes the stored elastic energy only minimally. This suggests that the crystal is left with considerable amount of dislocations which may not be strongly correlated because the stress relaxes only minimally.

In triangular lattices with hexagonal crystallographic symmetry, the macroscopic mechanical response is different. For two loading paths with φ = 0 and φ = π/3, passing inside the low energy valleys, the overall softening nonlinear elastic response terminates with a yielding avalanche which takes place at exactly the same level of energy and the same value of the loading parameter. The pristine state is not recovered after the discontinuous yield because higher symmetry gives rise to geometric frustration, and the stored elastic energy cannot be fully resolved through dislocation self-organization. However, the resulting stress drop is considerable, which suggests that long-range correlations have been created. The slight difference in the terminal post-avalanche state for the paths with φ = 0 and φ = π/3 is a finite size effect suggesting that the corresponding instability modes interact differently with the incommensurate square shape of the sample.

Along the non-generic 'hard' loading path with φ = π/2, the crystal shows hardening nonlinear elastic response, however, the yielding stress is only slightly higher than in the case of the 'soft' paths considered above. The discontinuous yield leads again to relatively small energy relaxation, suggesting that a high symmetry environment favors the considerable accumulation of dislocations. However, the stress drop is smaller than along more symmetric 'soft' paths, which hints towards a higher level of frustration and weaker long-range correlations in the post-buckling equilibrium state.

To summarize, our numerical experiments suggest that nominally ductile crystals, which are expected to yield continuously in the bulk form (large N , pre-existing defects), can still undergo a brittle-to-ductile (BD) transition when the value of N and the initial dislocation density both drop below certain thresholds. Such brittleness has been indeed realized to be a characteristic feature of sub-micron samples with high initial purity (absence of solutes, precipitates, and dislocations). The implied discontinuous yield results from massive homogeneous nucleation of dislocations in the form of a highly cooperative avalanche. As we have shown, the MTM can successfully simulate such avalanches and captures not only the catastrophic stress drop but also, in some cases, the recovery of almost pristine post-avalanche states. Even more importantly, MTM allows one to study the subtle quantitative dependence of the parameters of such discontinuous yield on crystal symmetry and sample orientation. Some features of the observed behavior may depend on our particular choice of the Landau potential, like for instance, the crossover from elastic softening to elastic hardening depending on the orientation of the sample; the question whether such behavior is instead a generic property of all GL(2, Z) symmetric potentials will be addressed elsewhere.

Ideal shear strength

To rationalize the observed behavior, we now study analytically the linear stability of an affine elastic response in homogeneously deformed crystals. From the classical continuum elasticity standpoint, the instability would mean that a homogeneous configuration is no longer a weak local minimum of the elastic energy. In a hard device, such instability results from the local loss of rank-1 convexity of the elastic energy density which also means the loss of strong ellipticity of the equilibrium equations [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Grabovsky | Normality Condition in Elasticity[END_REF]. The instability in the regularized model is necessarily delayed which is one of the sources of the 'smaller is stronger' size effect in sub-micron crystals.

Consider a homogeneous configuration of an elastic crystal with deformation gradient F rigidly imposed on its boundary. Suppose that this state becomes linearly unstable at the critical value of the imposed strain Fc . To find this threshold, we need to linearize the equilibrium equations around the homogeneous state. If we write them in terms of small incremental displacements u superimposed on the homogeneous state F (used as the reference state), we obtain the system A pi q j (F)u j ,pq = 0, where A pi q j (F) = F pK F qL A K i L j is the (Eulerian) tensor of incremental moduli. The unstable mode can be represented as a combination of Fourier components u = η η η exp (i kn n nx), where n n n is the unit normal selecting the modulation direction, η η η is the amplitude, and k is the wave number. In the Fourier space, this equation can be re-written as

[Q Q Q(F,n n n)η η η]η η η = 0,
where Q j l (F,n n n) = A i j ml n i n m (F) is the Eulerian acoustic tensor [START_REF] Merodio | Material instabilities in fiber-reinforced nonlinearly elasti solids under plane deformation[END_REF][START_REF] Kumar | On the hyperelastic softening and elastic instabilities in graphene[END_REF]. Stability of the homogeneous state F is lost when there exists a non-trivial n n n such that

det Q Q Q(F,n n n) = 0.
This condition can be used to identify the unstable values of C [START_REF] Hill | Acceleration waves in solids[END_REF][START_REF] Rice | Localization of plastic deformation[END_REF][START_REF] Ogden | Non-linear elastic deformations[END_REF]. Generalizing the classical definition of Frenkel [START_REF] Anderson | Theory of Dislocations[END_REF], we can associate the ideal shear strength with such C and compute them by performing a sequence of monotone loading tests. The ensuing critical orientation n n n and the associated polarization η η η can be interpreted as the characteristics of the incipient defects, see, for instance, [START_REF] Van Vliet | Quantifying the early stages of plasticity through nanoscale experiments and simulations[END_REF][START_REF] Zhu | Predictive modeling of nanoindentation-induced homogeneous dislocation nucleation in copper[END_REF][START_REF] Zhong | Simulating nanoindentation and predicting dislocation nucleation using interatomic potential finite element method[END_REF][START_REF] Miller | On the nonlocal nature of dislocation nucleation during nanoindentation[END_REF][START_REF] Bigoni | Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability[END_REF][START_REF] Garg | Universal scaling laws for homogeneous dislocation nucleation during nano-indentation[END_REF]. To identify the sub-domain in the space of reduced metric tensors (inside the extended fundamental domain) where the homogeneous deformations are stable, we conducted a series of numerical experiments with differently oriented lattices loaded by shear deformations of the form (3) and (4). Interpolating over a family of one parametric paths we obtained the surfaces of ideal strength for both square and triangular lattices, see the light grey contours in Fig. 5(a) and Fig. 5(b), respectively. These contours represent upper bounds for the sets which can be interpreted as continuum yield surfaces.

Observe that such surfaces may have rather different shapes: more elongated in the case of less symmetric square lattices and more round in the case of more symmetric triangular lattices. To maintain the symmetry of the configurational space, we showed in the case of triangular lattice two symmetric yield surfaces. Similar stability boundaries can be, of course, constructed around each of the equivalent replicas of the unstressed lattices [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF].

The strong asymmetry of such yield surface in the case of square lattices, see Fig. 5(a), shows that if an anisotropic crystal is loaded in 'soft' and 'hard' tensorial directions it can exhibit very different strength. Thus, the instability can take place either inside the low energy valley or get delayed till after the crystal has first stored considerable amount of elastic energy. Instead, in high symmetry crystals the apparent yield surface is almost isotropic, see Fig. 5(b), and since in this case all tensorial loading directions are almost equally 'soft', the expected theoretical strength is low independently of the orientation of the applied shear.

In our numerical experiments, the parameters of discontinuous yield were practically indistinguishable from the computations based on the idea of ideal shear strength, see small stars on the stress-strain and energy-strain graphs in Fig. 4. This means that inside the elastic range the system with N = 1024 is adequately described by the continuum limit. This also justifies the experimental observations that pristine sub-micron crystals yield close to the limits of theoretical strength which therefore cannot be improved without changing the nature of the crystal. It is instructive to look in more detail at the instability modes for crystals reaching different points α = α c (φ) along the surface of ideal strength. The knowledge of the corresponding vectors n n n c and η η η c allows one to identify the primary instability modes which, at least for sufficiently large N , guide the eventual nucleation of dislocational dipoles. Since the value of the bulk modulus µ in the strainenergy density was chosen to be sufficiently large, our crystals are effectively incompressible and therefore we can always expect that approximately n n n c ⊥ η η η c and the instability mode is close to a simple shear.

α α θ θ α θ (α c , θ c ) (α c , θ c ) (α c , θ c ) det(Q) det(Q) det(Q) (a) (b) (c) (d) (e) (f)
Consider first square lattices, and assume that the macroscopic loading paths are again given by (3) with φ = 0, arctan(1/2), π/4, see Fig. 2. In Fig. 6(a-c) we show the evolution of the determinant of the Eulerian acoustic tensor along each of these path as a function of the parameter θ defining the orientation of the unit vector n n n = (cos θ, sin θ). The external boundaries of the dark blue regions define the stability limits and the smallest value of the loading parameter α on one of such boundaries, defines the maximal homogeneous strain achievable along the corresponding loading path. Note that the whole pattern is periodic in θ and only one period is shown in Fig. 6(a-c).

We recall that the 'soft' path with φ = 0 corresponds to a simple shear along one of the dense planes of the square lattice. According to Fig. 6(a), at the critical value of the loading parameter α c ≈ 0.134, the associated θ c ≈ 3.035rad. Since the unit vectors n n n c and η η η c ≈ n n n ⊥ are almost parallel the current (deformed) lattice vectors, see Fig. 6(d), the unstable mode essentially activates a single slip in the 'vertical' direction.

The fact that this first unstable mode is not the one associated with θ c ≈ π/2, which would mean a single slip in the 'horizontal' direction is interesting because that is what one would expect based on the naive application of the plastic 'mechanisms' approach of CP. Note, however, that the mode with θ c ≈ 1.5746rad, which corresponds to the anticipated 'horizontal' direction of slip is destabilized almost simultaneously. We remark that the splitting of these two modes is the effect of geometric nonlinearity. Indeed, if kinematics is linearized, which essentially means that the parameter α is small, simple shears along the tensorial directions φ = 0 and φ = π/2 become indistinguishable when written in terms of the strain tensor E(α, φ)

= (1/2)(F T (α, φ)F(α, φ) -I) since E(α,0) = 1 2 0 α α α 2 E(α,π/2) = 1 2 α 2 α α 0 . ( 11
)
The geometrically linearized theories, neglecting the terms quadratic in α, are essentially confined to the plane tr (C) = 2, which is tangential to the hyperboloid det C = 1 at S 0 . In such theories, see, for instance, [START_REF] Onuki | Plastic flow in two-dimensional solids[END_REF][START_REF]Discrete models of dislocations and their motion in cubic crystals[END_REF][START_REF] Minami | Nonlinear elasticity theory of dislocation formation and composition change in binary alloys in three dimensions[END_REF][START_REF] Geslin | Investigation of coherency loss by prismatic punching with a nonlinear elastic model[END_REF], the paths φ = 0 and φ = π/2 would merge into a single line

C 11 = C 22 = 1.
As a result, the two minima S 1 and S 1 would collapse on each other. In our geometrically nonlinear approach the applied simple shear biases one of these modes and the two plastic 'mechanisms' end up activated consecutively rather than simultaneously. Similar degeneracy occurs also in the case of triangular lattices if the configurational space is reduced to the plane tangent to the point T 0 as in [START_REF] Onuki | Sheared solid materials[END_REF]. For our two other loading paths with φ = arctan(1/2) and φ = π/4, see Fig. 6(b,c), a single unstable mode can be clearly isolated, however, the vectors n n n c and η η η c ≈ n n n ⊥ c are again almost parallel to the deformed lattice vectors, see Fig. 6(e,f). Therefore, in both cases, a single slip system is activated ('vertical') and initially nucleated dislocation dipoles will be of one type only. Here we imply that the post-bifurcational localization will take the form of the nucleation of dislocation pairs along the plane selected by the condition of continuum instability as in Peierls-Nabarro model [START_REF] Nabarro | Dislocations in a simple cubic lattice[END_REF]. We also observe that in the hard device loading conditions, the preference of the 'horizontal' slip system is only apparent even for our 'quasi-horizontal' loading directions, which in fact become progressively 'harder' as we move away from the 'soft' path φ = 0. Along such paths the instability takes place at higher values of energy density and since the unstable mode is not perfectly aligned with deformed lattice vectors the ensuing geometrical frustration may lead to higher complexity of the resulting dislocation patterns.

α α θ θ α θ (α c , θ c ) (α c , θ c ) (α c , θ c ) det(Q) det(Q) det(Q) (a) (b) (c) (α c , θ c ) (d) (e) (f)
In the case of triangular lattice, we again consider non-generic simple shear loading paths along the dense crystallographic planes with φ = 0 and π/3 and a generic path with φ = π/2, see [START_REF] Anderson | Theory of Dislocations[END_REF]. The first two represent 'soft' shearing directions, while the last one is the 'hard' one facing potentially higher energy barrier. Figure 7(a,b,c) show the level sets of the determinant of the Eulerian acoustic tensor in the plane of parameters α and θ for each of the three paths.

Note first that, in contrast to the case of square symmetry, here the instability thresholds along different paths do not differ considerably. Moreover, the apparently 'hard' path with φ = π/2, which does not correspond to any straightforward plastic 'mechanism', gets destabilized before the apparently 'easy' ones with φ = 0 and φ = π/3. We reiterate that such quasi-uniformity of the instability conditions is a natural property of higher symmetry lattices with all thresholds collapsing onto one in the case of ideal isotropic solids (with superimposed fluctuations in the case of amorphous glasses or polycrystals).

We also observe that in all three cases the instability modes are characterized by the vectors n n n c and η η η c ≈ n c n c n c ⊥ that are no longer co-linear with the orientation of the deformed lattice basis. This suggests that the initial patterning, controlled by the continuum instability, is incommensurate with lattice slips and therefore the mechanism of dislocation nucleation may be rather different from the one predicted by the Peierls-Nabarro model [START_REF] Nabarro | Dislocations in a simple cubic lattice[END_REF]. The fact that still η η η c ≈ n n n ⊥ c suggests that the macroscopic will appear in the form of simple shear, however, the misalignment of the instability mode with crystallographic slip planes makes the prediction of the actual activated slip planes difficult. To show the effect of finite N , we also studied the stability limits in the original discrete setup, see also [START_REF] Sanderson | Armadillo: a template-based C++ library for linear algebra[END_REF][START_REF]A User-Friendly Hybrid Sparse Matrix Class in C++[END_REF]. To this end we computed the smallest eigenvalue of the 2N 2 × 2N 2 Hessian matrix K ab i j . Our Fig. 8 shows that although the critical value of the loading parameter α * c (N ) in the discretized problem, chosen by the condition that such eigenvalue is equal to zero, is slightly larger than the ideal shear strength α c in the continuum problem, the two approach quickly as the element size h = 1/N decreases. In the case of square lattices the gap closes up already at N ∼ 100 as the instability directions are almost perfectly aligned with the deformed crystal, see Fig. 8(a). In the case of triangular lattices, the misalignement is much stronger and the finite size effect in the form of a gap between the predictions of discrete and continuum theories remains apparent for much smaller element sizes, see Fig. 8(b). Note that the asymptotic behavior of the smallest eigenvalues near zero, shown in Fig. 8, is different from the prediction of the theory of amorphous plasticity where the corresponding eigenvectors are quasi-localized [START_REF] Dasgupta | Universality of the plastic instability in strained amorphous solids[END_REF][START_REF] Bonfanti | Elementary plastic events in amorphous silica[END_REF][START_REF] Richard | Predicting plasticity in disordered solids from structural indicators[END_REF]. The linear response of pristine crystalline solids with zero disorder is global and that is why its elastic instabilities can be largely captured by the classical continuum theory.

The dependence of the instability threshold α * c on the system size N is illustrated in Fig. 9. The theoretical limit, effectively reached by the samples with N = 1024, is shown by the red dashed line. At small N we observe the emergence of a nontrivial asymptotics α * c ∼ N -1/3 . It suggest that the spatial scale of the dislocation microstructure is of the order ∼ h 1/3 L 2/3 , which implies hierarchical 'domain splitting' near the external boundaries. For square and triangular lattices with N = 40 we illustrate in Fig. 10(a,b), the 2N 2 dimensional eigenvectors u a i of the corresponding Hessian matrices K ab i j when their lowest eigenvalues crosses zero at α * c . These eigenvectors are presented as vector fields with the homogeneous component subtracted. In both cases, the modulations have orientations n n n * c that are close to the ones predicted by the continuum theory. However, if the unstable wavenumber remains arbitrary in the scale-free continuum theory, the discrete theory selects a particular length scale. The strong misalignment, in the case of triangular lattices between the orientation of the macro modulations and the lattice vectors, see Fig. 10(b), can delay the transition to the ultimate lattice scale instability culminating in the nucleation of dislocation dipoles. In the case of square lattices, shown in Fig. 10(a), we observe lattice scale modulations present already in the original unstable mode. They correspond to the wave vectors at the boundary of the Brillouin zone, see some elementary examples of such multiscale instabilities in [START_REF] Truskinovsky | Quasicontinuum modelling of short-wave instabilities in crystal lattices[END_REF][START_REF] Bertoldi | Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformations[END_REF].

To summarize, our numerical experiments show that pristine sub-micron crystals exposed to affine deformation yield discontinuously near the thresholds of theoretical strength. In the case of simple shear loading, the obtained thresholds show dependence on the orientation of the sample in the loading device, which is in qualitative agreement with experimental results, e.g. [START_REF] Ziegenhain | Influence of crystal anisotropy on elastic deformation and onset of plasticity in nanoindentation: A simulational study[END_REF][START_REF]Effect of crystal orientation on the size effects of nano-scale fcc metals[END_REF]. Our study also reveals a strong link between the structure of the apparent yield surface and the lattice symmetry.

Discontinuous yielding

In this Section, we turn to the study of the fine structure of the microscopic response following the loss of stability of the homogeneous state. The goal is to reveal the detailed unfolding of the catastrophic dislocation nucleation avalanche and trace how the final configuration, containing a large number of defects, emerges in the process of energy minimization. In particular, we show that complementary pictures of the evolving defect configurations emerge in the real and the configurational spaces.

All numerical experiments start with a dislocation-free crystal u = 0, and we drive the system using an athermal quasistatic protocol. We use hard device boundary condition controlling the positions of surface nodes and impose in this way an affine deformation F on the boundary of a square domain and choose the loading step (strain increment) sufficiently small to ensure continuity of the minima outside the avalanches. The macroscopic stress-strain and energy-strain curves obtained in such numerical experiments are summarized in Fig. 4(a) and Fig. 4(b) for the square and triangular crystals, respectively.

For α < α * c the picture is rather simple. In real space, the deformation is affine F(x) ≡ F and in the configurational space of the metric tensors C we observe the perfect overlap for all elements. As we have already mentioned, there is an excellent agreement between the numerical value of α * c and the predicted value of ideal shear strength, see Fig. 4. As soon as stability is lost, the homogeneous configuration breaks down. The pristine crystal transforms into a highly defected one, with most of the dislocations eventually escaping from the bulk and forming various pile-ups near the fixed boundaries. The configurational points, corresponding to different elastic elements, spread in the space of metric tensors C, reaching the neighborhoods of distant energy wells. The emerging dislocation patterns, which we discuss in detail below, were obtained from the simulations involving more than a million elements.

Square lattices

Consider first a square crystal loaded in a crystallographically exact simple shear with φ = 0. In Fig. 11, we show the snapshots of the distribution of the shear component of the Cauchy stress σ = (det F) -1 FP T in the physical space as the crystal evolves (in fast computational time) through the avalanche after the loading parameter has reached the value α * c . In the insets, we illustrate the concurrent evolution of the cloud of configurational points in the space of metric tensors which tracks different stages of the implied energy minimization process.

At positive values of the loading parameter α the crystal is driven along the energy valley away from the energy minimum S 0 towards the minima S 1 , S 11 , etc. This could suggest the development of plastic slip along the horizontal plane, however, as we have seen, the first instability is 'vertical' which means activation of the slip system represented by the minima S 2 , S 22 , etc., see Fig. 6. While the initial states of the avalanche are indeed dominated by the activation of the vertical dislocation dipoles, see Fig. 11 (a), horizontal dipoles appear later as well because of local stress concentration and consequent involvement of the second slip system predicted by the instability analysis to be activated at slightly larger loading levels, see Fig. 11(b).

Note that the two low energy valleys merge around the point T 0 , which describes a triangular lattice. The structure of the energy landscape around this point is close to a monkey saddle and passing of configurational points through such saddle creates non-trivial coupling between horizontal and vertical slip systems [START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF]. As a result, even though most of the units move towards the energy well S 2 some of the configurational points eventually also populate the energy well S 1 . As the crystal is driven further, it is pulled towards another monkey saddle describing the triangular lattices T 1 . Around this intersection of the low energy valleys, the configurational points take a turn towards the energy well S 12 instead of continuing towards apparently more natural well S 11 . Similarly, around the monkey saddle T 2 the configurational flow is directed towards the energy well S 21 , see Fig. 11(bd). All these choices depend sensitively on the boundary conditions and will be rationalized in a separate publication. To see how the spreading of the configuration points in the space of metric tensors transforms when we move into the physical space, we focused in Fig. 12 on three randomly selected fragments of the post-avalanche equilibrium configuration shown in Fig. 11(d). The emerging equilibrium pattern contains dislocational pile ups near the boundaries and few locked up dislocational dipoles in the bulk of the specimen. Note, in particular, the pairs of dislocations on two parallel slip planes blocking each other, see Fig. 12(b), and on two perpendicular slip planes forming characteristic locks, see Fig. 12(d). The fact that after the system spanning avalanche most of the nucleated defects manage to either annihilate or escape to the boundaries of the crystal, explains why the crystal becomes almost pristine after such a devastating discontinuous yield event.

To check the generality of these conclusions, we also performed numerical simulations of the discontinuous yield at α ∼ α * c with initially perturbed displacement field which was 'dirtied' by random Gaussian perturbations with zero mean and small variance (∼ 0.0001). Such disorder is too small to suppress brittleness [START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF] and, we show in Fig. 13 the equilibrium stress field after the catastrophic avalanche. Since the pre-avalanche state was not 'pristine', the energy minimizing self organization of dislocations is compromised by the geometrical frustration imposed by the disorder. The resulting dislocation pattern appears to be structurally similar at several scales, see the successive insets shown in Fig. 13(b,c). The two main perpendicular slip planes are now activated at almost the same level. Note the formation of self screening dislocation-rich walls, reminiscent of low-angle grain boundaries, which separate dislocation-poor, low-stress domains where the lattice undergoes simple rotations, see small rotated squares in Fig. 13(d); all these structural features have been previously recorded in physical experiments. Note also, that since samples containing imperfections show diminished brittleness [START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF], the discontinuous yield can disappear in cyclic loading already after few cycles.

Consider next a generic 'hard' shear loading path with φ = arctan(1/2) which does not follow the low energy valley. As we have already mentioned, in this case a homogeneously deformed crystal will store a bit more elastic energy before instability than in the case of the symmetric 'soft' paths. The breakdown scenario is shown in Fig. 14: after reaching the (regularized) ideal strength limit

F(α * c , arctan 1
2 ), dislocations nucleate along the 'vertical' slip plane which agrees with the prediction C. R. Physique, 0000, 1, n o 0, 000-000 of the stability analysis. In the configurational space of metric tensors the stream of points describing such dislocations is directed from S 0 to S 2 , see Fig. 14(a). Eventually few 'horizontal' dislocations form as well, see Fig. 14(b,c), and because of geometrical frustration they are not all expelled to the boundaries of the crystal, see Fig. 14(d). Therefore the system is left with considerable residual energy and the stress drop is relatively small. The most interesting yield scenario is obtained for the non-generic φ = π/4 loading path which leads from the reference state S 0 to the second closest energy minimum S after passing a much higher energy barrier than in the case φ = 0. Here even more elastic energy is stored before the instability than in the case φ = arctan(1/2). The breakdown takes place at F(α * c , π/4) with a single instability mode activated, see Fig. 15 (a). However, the unfolding collective dislocation nucleation process pattern quickly becomes very complex. Almost instantaneously, in addition to the original energy well S 0 , all four neighboring square wells (S 1 , S 1) and (S 2 , S 2) become engaged and few elements even reach more distant square wells, see Fig. 15(b,c). The flow of configurational points in the space of metric tensors apparently uses the monkey saddles 3 corresponding to the triangular lattice configurations T 0 and T 1 to access the four nearest square energy wells. The compositions of the shear strains from the two main slip systems, leading to the square energy wells S 11 and S 11 , are getting activated as well. The post-avalanche dislocation pattern is then multi-slip with minimal pile up on the boundaries and most of dislocations self-organizing in the bulk of the crystal, see Fig. 15 (d). A close-up of the highlighted region shown in Fig. 16(a-b) down to the scale of finite element nodes puts in evidence the formation of dislocation walls separating almost unstressed cells with correlated lattice misorientations. In other words, we observe the formation of dislocation-rich bands which self screen long-range elastic fields allowing the strains inside the dislocation-poor regions to relax by reducing the deformation gradient to pure rotation. In Figure 16(c), we show such differently re-oriented crystal (almost) stress free regions in more detail. Interestingly, molecular dynamics simulations also suggest that at almost theoretical strength, dislocations alone can no longer release the elastic energy, and the reorientation of crystal lattices emerges as an important mechanism of plastic response [START_REF] Zepeda-Ruiz | Probing the limits of metal plasticity with molecular dynamics simulations[END_REF]. In other words, a perfect sub-micron crystal appears to be yielding by catastrophic strain-induced reorientation of the crystal lattice.

To summarize, our numerical experiments conducted on square lattices suggests that while in specially oriented crystals, pristine-to-pristine brittle yield at ideal shear strength threshold is a possibility, the generically oriented, dislocation free sub-micron crystals, can be expected to yield discontinuously with massive dislocation nucleation culminating in the formation of complex dislocational patterns. During such instabilities a large amount of stored elastic energy is released and the connection between the initial instability mode and the developing deformation pattern is quickly lost. To minimize elastic energy, dislocations do not only nucleate cooperatively but also self organize hierarchically forming multi-scale cell structures with presumably complex statistical properties, to be studied separately using, for instance, the approach presented in [START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF].

Triangular lattices

Consider next the case of higher symmetry triangular lattices. As we have already seen, in this case the macroscopic mechanical response is in overall agreement with the analytical yield surface, modulo the fact that, differently from what we have seen in the case of square lattices, discontinuous yield is observed for values α * c which were about a 10% higher than the corresponding analytical predictions α c , see Fig. 8(b). This gap, which is a function of the finite element size h, results from the disregistry between the orientation of the unstable modes and the 'soft' directions of the energy landscape.

Consider first the loading path corresponding to a simple shear (4) with φ = 0. An application of the idea of plastic 'mechanisms' suggests that crystal should evolve from the initial energy well T 0 to the next closest triangular well T 1 ; the corresponding low energy valley is shown in Fig. 3 as the path 1-1. However, as we have already seen, the actual instability of the homogeneously deformed lattice state T 0 occurs strictly after the determinant of the continuum acoustic tensor has become negative for some orientation n n n c . To deal with encountered soft directions, the implied dislocation nucleation is postponed. The anticipated low energy valleys in the energy landscape cannot be used because the apparently low barriers leading towards the state T 1 along the path 1 -1 are still too high. Instead the softness of the energy landscape drives the system in the direction of the saddle point S 1 which corresponds to the square lattice while still remaining within the fundamental domain. Therefore in our numerical experiments we observed that during the gradient descent type energy minimization at α c that periodically-spaced modulations develop prior to the dislocation nucleation, see Fig. 17(a). The orientation of the modulation bands agrees with the directions n n n c obtained using the continuum instability condition.

As the instability develops, the non-equilibrium configurations found by the minimization algorithm show the sharpening of the band boundaries and then the secondary symmetry breaking through the homogeneous nucleation of dislocation dipoles, see Fig. 17(b). Note that the orientation of these dipoles agrees with the lattice (T 0 to T 1 path) but not with the orientation of the bands (T 0 to S 1 path). Note also that the implied pairs of dislocations with opposite signs nucleate at the centers of the modulation bands where the displacement gradients are large which apparently helps the system to overcome the energy barriers leading to T 1 lattice state. The emerging dislocation nucleation scenario is slightly different from the one suggested by the Peierls-Nabarro model where, effectively, the bands are assumed to be atomically sharp and lattice oriented. In Fig. 17(b-d)) we see that the modulation eventually completely breaks down with only one slip system ending up being activated. Due to the misalignment of the bands with the slip planes and the associated geometrical frustration, the nucleated dislocations interact with each other strongly causing the activation of the double slip represented by the triangular lattice state T 11 . With dislocations either annihilating or escaping to the boundaries, the resulting picture is in basic agreement with the prediction of a single slip CP theory.

Along the other symmetric loading paths with φ = π/3 and φ = 2π/3, which also describe shears on dense planes of the triangular lattice, the general scenario of discontinuous yield remains the same. While the degree of misalignment of the initial modulation bands with the lattice does not change, the orientation of these bands relative to the square computational box is slightly different in each of these cases which (mildly) affects the outcome of the instability. The succession of events, constituting the system spanning avalanche in the case of loading orientation φ = π/3, is illustrated in Figs. 18(a-d). Here again the macroscopic modulations create the nucleation sites responsible for the secondary instability which leads to the formation of almost ideal dislocational dipolar mats [START_REF] Hansen | Low energy dislocation structures due to unidirectional deformation at low temperatures[END_REF], see Figs. 18(b). Dislocations eventually self organize into bands, see Figs. 18(c), but eventually all end up on the surface of the crystal, see Figs. 18(d).

Here again the primary, purely elastic instability develops in the direction from triangular well T 0 to the square well S 0 . The secondary instability associated with barrier crossing leads to the flow of configurational points from T 0 to T 1 which therefore creates plastic slip. As dislocations pile up near the boundaries, they activate the secondary slip associated with the well T 11 and even the second slip system represented by the well T 12 . Ultimately, all four triangular lattices T 1 , T 2 , T 1 1 and T 12 get involved, however, the ensuing complexity remains localized near the sharp corners of the sample. At the end of the avalanche the bulk of the sample appear to be free of dislocations which explains the dramatic stress drop experienced by the crystal. In view of the high symmetry of the triangular lattice, the structure of the yielding avalanche for the generic path with φ = π/2, which also starts in the original energy well wells T 0 , is not very different from what we have already seen in the case of the non-generic path with φ = 0. The influence of the neighboring square well S 0 is again felt in the structure of the initial elastic modulations, see Fig. 19(a). The deformation ultimately localizes inside the low energy valley connecting the triangular wells T 0 and T 1. Transient dislocation entanglements eventually get dissociated with dislocations either annihilating or escaping to the boundaries, however, the stored elastic energy is not fully relaxed by such limited plasticity, see Fig. 19(b-d). In general, this example suggests that in sub-micron crystals with high symmetry and accordingly 'round' surfaces of theoretical strength, the outcome of the discontinuous yield is decided not so much by the orientation of the sample in the loading machine, but rather by the orientation of the sample boundaries vis-a-vis lattice directions. To summarize, using MTM, we were able to show that, similar to what is observed in experiments, pristine sub-micron crystals undergo discontinuous yielding close to theoretical strength. During the catastrophic avalanche, massive nucleation of dislocations first transforms such crystals from defect-free to defect-saturated, but then most of the defects either annihilate or escape to the fixed boundaries producing again almost pristine crystals. The access through MTM to transient nonequilibrium states allowed us to study for the first time the sensitive dependence of the phenomenon of discontinuous yield on crystal symmetry, crystal orientation in the loading machine, and even crystal shape. The general feature of the simulated discontinuous yield is that the instability starts as a long-wave purely elastic modulation which breaks the affinity of the original homogeneous state. The actual dislocation nucleation emerges as a secondary instability, triggered by non-affinity and proceeding collectively. A large number of defects appear almost simultaneously and the 'auto-catalytic' nature of such avalanche is due to the fact that already nucleated dislocations trigger new nucleation events. The process terminates only when the elastic repulsion from the walls finally blocks the nucleation in the bulk.

The pristine-to-pristine yield takes place if a non-generic loading allows the system to completely avoid frustration and channel almost all dislocations away from the bulk. Otherwise the avalanche jams with dislocations forming self-locking patterns. Such patterns can be interpreted as cell structures because relatively narrow dislocation walls bound extended domains where lattice undergoes simple rotation. In the case of open boundaries, the dislocation pile-up, playing an important role in this scenario, will be replaced by the global strain localization. Here N = 500 and φ = 0.

Beyond the principal avalanche

The analytic and numerical results presented in the previous Sections are relevant for conceptualizing the puzzling discontinuous response of nominally ductile crystals. Such brittle-like response is generic for ultra-small and, therefore, dislocation-starved crystals in the sense that it is observed routinely in compression and nano-indentation experiments [START_REF] Corcoran | Anomalous plastic deformation at surfaces: Nanoindentation of gold single crystals[END_REF][START_REF] Bei | Effects of focused ion beam milling on the nanomechanical behavior of a molybdenum-alloy single crystal[END_REF][START_REF] Chisholm | Dislocation starvation and exhaustion hardening in Mo alloy nanofibers[END_REF][START_REF] Wang | Pristine-to-pristine regime of plastic deformation in submicron-sized single crystal gold particles[END_REF].

As we have shown, to ensure such pseudo-brittle behavior, the crystals should be structurally perfect (dislocation and defect-free). In this case, the main effect is an explosive system spanning avalanche at the end of the elastic regime with massive collective nucleation of dislocations, which results in a catastrophic load drop known as discontinuous yield. In this Section, we briefly discuss what happens if the loading resumes and the crystal, whose purity is now compromised, continues to yield. We limit the analysis to just a few observations, while the detailed study will be presented elsewhere.

As we have seen, after the termination of the system size avalanche and the re-stabilization of the inhomogeneous state, the generic crystalline sample can be expected to contain a highly correlated configuration of lattice defects. If the quasi-static loading continues, such an imperfect crystal will first undergo another quasi-elastic deformation which will be succeeded by another avalanche, necessarily less dramatic in view of the ubiquity of nucleation sources and the availability of the locking sites. Subsequent monotone loading will be punctuated by an intermittent succession of such avalanches, with most of them small, but some reaching again the size of the whole system [START_REF] Friedman | Statistics of dislocation slip avalanches in nanosized single crystals show tuned critical behavior predicted by a simple mean field model[END_REF][START_REF] Zaiser | Statistical aspects of microplasticity: experiments, discrete dislocation simulations and stochastic continuum models[END_REF][START_REF] Ispánovity | Avalanches in 2D dislocation systems: plastic yielding is not depinning[END_REF][START_REF] Derlet | The stress statistics of the first pop-in or discrete plastic event in crystal plasticity[END_REF][START_REF] Sethna | Deformation of Crystals: Connections with Statistical Physics[END_REF][START_REF] Papanikolaou | Avalanches and plastic flow in crystal plasticity: an overview[END_REF][START_REF] Cui | Spatio-temporal plastic instabilities at the nano/micro scale[END_REF][START_REF] Sparks | Avalanche statistics and the intermittent-to-smooth transition in microplasticity[END_REF][START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF]. To illustrate these general observations and show that MTM approach is well suited to study the effects of intermittent discontinuous yield, we now present few illustrations of the mechanical response of square crystals subjected to simple shear with φ = 0 following directly the initial plastic avalanche. Since, naturally, the complexity of the dislocational patterns increases with subsequent loading, we reduce the precision of the description by adopting a smaller value of N = 512 but preserve the small annealed disorder as in the experiment shown in Fig. 13.

The simulated macroscopic response is summarized in Fig. 20, where we show the computed We first observe that after the major stress/energy drop, the smooth elastic response turns into jerky plastic yield. More specifically, the post-catastrophic response can be decomposed into conservative elastic steps that are interdigitated by dissipative stress/ energy drops representing dislocation avalanches broadly distributed in size, see Fig. 20 (b,d). While the averaged stress response in this range shows monotone hardening, the averaged energy first decreases and then increases. This is due to the fact that the energy losses cannot even be seen on such stress-strain curves, see [START_REF] Puglisi | Thermodynamics of rate-independent plasticity[END_REF] for more details. The jerkiness of the response originates in transitory elastic deformation of self-locked microstructures which always ends with partial or complete unlocking and dynamical restructuring, accompanied by energy dissipation.

To trace such restructuring, we show in Fig. 21(a-d) a succession of four snapshots of the stress field in the arbitrary chosen moments of 'time' as indicated in by dots marked by numbers 1, 2, 3, 4 in Fig. 20(b,d). Afterwards, in Fig. 22 we zoom into two spatial regions A and B and follow the 'time' evolution of the particular groups of dislocations over several avalanches. The numerical experiment, illustrated by these figures, is similar to the one shown in Fig. 12 but with even smaller annealed initial disorder. We observe that continuing loading is accommodated by additional nucleation of dislocations which remain in the bulk due to apparent repulsion from the already saturated boundaries. The configurational points in the space of metric tensors spread more uniformly among the three equivalent square energy wells S 0 , S 1 and S 2 , being distributed almost evenly by the monkey saddle T 0 . The ensuing dislocation pattern includes individual dislocation locks as well as dislocation walls separating domains where the lattice is almost unstressed. This pattern results from dynamic self-organization of dislocations and can be viewed as generic.

Upon further loading, the crystal continues to harden as seen in Fig. 20(d). The plastic yield takes the form of an irregular sequence of stress-drops associated with partial restructuring of the dislocation pattern, see Fig. 21(b-d). The configurational points continue to spread over the metric space reaching distant wells corresponding to various compositions of simple shears. This means that despite the simple shear strain applied on the boundary, the deformation in the bulk of the sample takes the manifestly multi-slip form. Interestingly, the patterns identified in Fig. 21(b-d) are reminiscent of the ones obtained by x-ray diffraction method in [START_REF] Jakobsen | Formation and subdivision of deformation structures during plastic deformation[END_REF] where the authors observed the formation of dislocation boundaries separating (nearly) dislocation-free regions with almost perfect lattices. Such subgrains were shown to exhibit intermittent dynamics, appearing and disappearing as well as displaying transient splitting behavior. Note that most small-sized avalanches are associated with the transitional motion of depinning dislocations between the dislocation-rich wall-like patterns. During such transitions, dislocations either get locked again or annihilate, however, at least some of the existing large-scale dislocation structures persist during the loading. To illustrate this conclusion, we show in Fig. 22 the enlarged view of the time evolution for the two rather arbitrarily selected square sub-domains, which we identify as A and B . In A, we see an almost rectangular dislocation structure based on dipolar walls, which ensure an ideal screening of the long-range elastic fields [START_REF] Zaiser | Statistical aspects of microplasticity: experiments, discrete dislocation simulations and stochastic continuum models[END_REF]. With loading, this structure deforms elastically but otherwise remains relatively stable, experiencing only minor rearrangements (addition and subtraction of dislocations from local entanglements). Instead, the incomplete dipolar wall shown in B does not survive the loading, and finally, the participating dislocations disappear entirely due to an abrupt collective annihilation. Usually, such restructuring events give rise to bigger avalanches, and we checked that the largest avalanche between the deformation states (c) and (d) shown in Fig. 20(b,d), can be indeed linked to the disappearance of the multi-dislocation structure B , see Fig. 22.

We observe, however, that the avalanches associated with partial restructuring of the dislocation pattern, are still much smaller in size than the system spanning avalanche responsible for the major restructuring behind the original discontinuous yield of pristine crystal. This is associated with a 'special preparation' of the pristine (homogeneous) crystal whose high degree of degeneracy causes instability to happen simultaneously all over the sample. The dynamic randomization produced by the discontinuous yield generates annealed disorder and reduces the degeneracy, producing more conventional yield represented by mostly small avalanches. However, due to long-range elastic interactions, strong correlations in such apparently random system remain, and they lead to occasional recurrence of the system size avalanches, known as 'dragon-king event' [START_REF] Sornette | Dragon-kings: Mechanisms, statistical methods and empirical evidence[END_REF]. This leads to the supercritical avalanche distribution in sub-micron crystals which was observed in experiment and has been recently successfully simulated in the framework of a scalar version of MTM [START_REF] Zhang | Variety of scaling behaviors in nanocrystalline plasticity[END_REF].

Conclusions

In this paper, we presented the first systematic demonstration of the effectiveness of the mesoscopic tensorial model (MTM) by showing that it can effectively simulate complex plastic phenomena in crystals involving a large number of interacting dislocations. A characteristic feature of the MTM is that it resolves lattice scale dislocation cores while operating with the engineering concept of stress and strain. It can be viewed as a tensorial version of the scalar phase-field model, accounting for large rotations and finite lattice invariant shears in geometrically precise way. Rather remarkably, the MTM can deal adequately with short-range interactions of dislocations while accounting for full crystallographic symmetry, without any direct use of interatomic potentials or other ab initio methods.

The key feature of the MTM is that it is formulated as a nonlinear field theory of Landau-type respecting the discrete GL(3, Z) symmetry. It can also be viewed as a geometrically and physically nonlinear anisotropic elasticity theory. It may first sound paradoxical that elasticity framework was chosen to describe plasticity. However, this viewpoint is fully consistent with the fact that dislocations are elastic defects whose interaction is largely elastic. Under quasi-static loading, inelastic dissipation takes place during fast switching of elastic branches. Such branch switching events are ubiquitous due to the rugged nature of the generic MTM energy landscape and the presence of elastic driving.

To generate such energy landscape, the MTM relies on local validity of the Cauchy-Born rule even though some features of the nonlocal atomistic description are then inevitably lost. Behind the emergence of the continuum stresses and strains is then the spatial averaging over a microstructural length scale: the latter characterizes the size of a cluster of atoms (an element) assumed to deform homogeneously. The elastic branches are switched when, as a result of elastic instabilities, such elements switch between the neighboring wells of the globally periodic Landau energy.

An important feature of the MTM approach to crystal plasticity is that plastic 'mechanisms' are not assigned a priori, but take the form of low energy valleys connecting the wells of Landau energy. Plastically deformed state emerges in this representation as a mixture of equivalent Landau phases, with dislocations appearing as incompatible components of phase boundaries. In such a model slip directions are not rigidly pre-defined and loading along one slip system can, by itself, initiate another slip system. Topological transitions, like the ones involved in the activity of Frank-Read and SA (single arm) sources, are also incorporated automatically.

In this paper we used the simplest 2D version of the MTM to simulate numerically a benchmark phenomenon: the catastrophic discontinuous yielding in pristine sub-micron crystals due to cooperative nucleation of dislocations. We showed that in simulations, the scenario of such 'brittle-like' dislocation nucleation event is similar to the one observed in experiments involving imperfectionfree sub-micron crystals. Moreover, the use of the MTM allowed us to quantify, for the first time, the dependence of the ensuing dislocation patterning on the crystal symmetry and on the orientation of the crystal in the loading device. The critical load was shown to be in agreement with the analytical predictions of the ideal shear strength which we associated with the violation of the Legendre-Hadamard (strong ellipticity) condition. We also showed that this criterion can successfully predict the observed orientation of the incipient non-affine patterns.

Vastly different dislocation nucleation scenarios were shown to characterize the discontinuous yield in square and triangular lattices. In the case of square symmetry, shearing in the 'soft' directions causes pristine-to-pristine transitions with stress dropping almost to zero as a result of the catastrophic avalanche. Instead, loading of the same crystals in the 'hard' directions leads to the formation of a complex, highly correlated dislocation pattern and results in residual stress. In higher symmetry hexagonal crystals, the theoretical strength threshold depends only weakly on the orientation of the crystal. Due to ensuing geometrical frustration, the nucleation of dislocations is preceded by the development of long-wave elastic modulations. While most of the dislocations manage to escape to the boundaries of the crystal, the post avalanche stress remains considerable independently of the sample orientation. In both square and hexagonal cases, lattice configurations with alternative symmetries (hexagonal and square, accordingly) appear as saddle points of the Landau energy and play an important role in fomenting multi-slip by re-directing the flow of configurational points between converging energy valleys.

Our study shows the potential of the MTM to deal with plastic flows, involving strongly interacting dislocations, while relying only on minimal phenomenological assumptions. It opens new paths in the study of spatial and temporal complexity associated with developed plastic flows. Already the first concrete results, obtained in this paper, represent promising steps towards harnessing brittle events in sub-micron crystals compromising forming processes, endangering the load-carrying capacity of micro-machine parts, and jeopardizing reliability in various other micromechanical applications.

Figure 1 .

 1 Figure 1. The structure of the GL(2, Z) periodicity domains in the space of metric tensors C. The presented portion of the infinite surface det C = 1 represents simple lattices with det C > 0. The projected section of the original hyperbolic surface is defined by the inequalities 0 < C 11 < 6, 0 < C 22 < 6, and -3 < C 12 < 3.

22 Figure 2 .

 222 Figure 2. (a) The level sets of the strain-energy density (2) with β = -1/4. Colors indicate the energy level: blue, low; red, high. Dashed lines correspond to the simple shear loading paths F(α, φ) defined in (3). (b) The periodic strain-energy profile along the two symmetric paths 1 -1 and 2 -2 and and the current state of lattice vectors.

  Figure 3. (a) The level sets of the strain-energy (2) with β = 4. Colors indicate the energy level: blue, low; red, high. Dashed lines correspond to the simple shear loading paths F(α, φ) defined in (4). (b) The periodic strain-energy profile along the three symmetric paths 1 -1, 2 -2, and 3 -3 and the current state of lattice vectors.

Figure 4 .

 4 Figure 4. Macroscopic mechanical response of square (a-c) and triangular crystals (b-c) subjected to simple shear loading: (a,b) the equilibrium elastic energy W(α) (c,d); the equilibrium stress defined as τ(α) = dW/dα = Ω P : (d F/dα)d x.

Figure 5 .

 5 Figure 5. Surfaces of ideal shear strength for square (a) and triangular (b) crystals. Silver contours delimit the domains of linear stability for homogeneous states. The strain-energy density is taken from (2).

Figure 6 .

 6 Figure 6. Instability limits for square lattices: (a-c) level sets for the determinant of the Eulerian acoustic tensor, (d-f) parameters of the unstable modes vis a vis the deformed lattice vectors. (a) φ = 0, (b) φ = arctan 1/2 and (c) φ = π/4. The parameters (α c ,θ c ) are calculated from the condition det(Q) = 0.

Figure 7 .

 7 Figure 7. Instability limits for triangular lattices: (a-c) level sets for the determinant of the Eulerian acoustic tensor, (d-f) parameters of the unstable modes vis a vis the deformed lattice vectors. (a) φ = 0, (b) φ = π/3 and (c) φ = π/2. The parameters (α c ,θ c ) are calculated from the condition det(Q) = 0.

Figure 8 .

 8 Figure 8. Smallest eigenvalue of the discrete stiffness matrix K as a function of the loading parameter α for different values of element size h = 1/N : (a) square, (b) triangular lattice. The crystal orientation is always φ = 0. In the continuum case, we show the smallest eigenvalue of the acoustic tensor Q Q Q(F,n n n c ).

Figure 9 .

 9 Figure 9. System size dependence of the critical loading parameter α * c (N ) in numerical experiments with triangular crystals oriented at φ = 0.

Figure 10 .

 10 Figure 10. Unstable modes at the instability point for square (a) and triangular (b) lattices. Blue arrows show the modes associated with the lowest eigenvalue of the stiffness matrix K. Here N = 40 and φ = 0.

Figure 11 .

 11 Figure 11. Snapshots of the stress field during discontinuous yield in pristine square crystals. Colors indicate the level of the shear component of the Cauchy stress. Insets show the same images in the configurational space. Black dots indicate the value of the metric tensor in individual finite elements. Here N = 1024 and φ = 0.

Figure 12 .

 12 Figure 12. Final stress configuration after the discontinuous yield in pristine square crystals: (b-d) show the enlarged versions of the square indicated by the same letters in (a). Here N = 1024 and φ = 0.

Figure 13 .

 13 Figure 13. Final stress configuration after the discontinuous yield in a square crystal which was randomly perturbed at the point of instability. Colors indicate the level of the shear component of the Cauchy stress. Squares show the rotated versions of the unstressed lattice. Here N = 1024 and φ = 0.

Figure 14 .

 14 Figure 14. Snapshots of the stress field during discontinuous yield in pristine square crystal. Colors indicate the level of the shear component of the Cauchy stress. Insets show the same images in the configurational space. Black dots indicate the value of the metric tensor in individual finite elements. Here N = 1024 and φ = arctan(1/2).

Figure 15 .

 15 Figure 15. Snapshots of the stress field during discontinuous yield in pristine square crystal. Colors indicate the level of the shear component of the Cauchy stress. Insets show the same images in the configurational space. Black dots indicate the value of the metric tensor in individual finite elements. Here N = 1024 and φ = π/4.

Figure 16 .

 16 Figure 16. Stress field after the avalanche for a pristine square crystal with orientation φ = π/4 (a); zoomed view of the marked area (b); further zoom to the scale of the discrete finite element nodes (c).

-0. 2 -Figure 17 .

 217 Figure 17. Snapshots of the stress field during discontinuous yield in pristine triangular crystal. Colors indicate the level of the Cauchy shear stress σ x y . Insets show the same images in the configurational space. Black dots indicate the value of the metric tensor in individual finite elements. Here N = 1024 and φ = 0.

Figure 18 .

 18 Figure 18. Snapshots of the stress field during discontinuous yield in pristine triangular crystal. Colors indicate the level of the shear component of the Cauchy stress. Insets show the same images in the configurational space. Black dots are placed according to the value of the metric tensor's components in individual finite elements. Here N = 1024 and φ = π/3.

Figure 19 .

 19 Figure 19. Snapshots of the stress field during discontinuous yield in pristine triangular crystal. Colors indicate the level of the shear component of the Cauchy stress. Insets show the same images in the configurational space. Black dots indicate the value of the metric tensor in individual finite elements. Here N = 1024 and φ = π/2.

Figure 20 .

 20 Figure 20. Macroscopic mechanical response of square crystal upon loading beyond the catastrophic system-size avalanche.

Figure 21 .

 21 Figure 21. Evolution of the shear component of Cauchy stress in physical space after the catastrophic system-size avalanche. The snapshots correspond to the selected points in Fig. 20 marked by numbers 1, 2, 3, 4. The corresponding distributions in the configurational space are shown in the insets. Here N = 500 and φ = 0.

Figure 22 .

 22 Figure 22. Enlarged view of the dislocation configurations inside the square regions (A) and (B) shown in Fig. 21 as they evolve with loading. Colors indicate the shear component of the Cauchy stress tensor.

Recently continuum models were proposed with internal fluctuations accounted for through stochastic constitutive equations[218,

223].C. R. Physique, 0000, 1, n o 0, 000-000

C. R. Physique, 0000, 1, n o 0, 000-000

Allowing one, for instance, to resolve the system size dislocation nucleation avalanche for micron size crystals.C. R. Physique, 0000, 1, n o 0, 000-000

Rather shallow energy maxima surrounded by the non-degenerate saddles describing rhombic lattices as shown in[START_REF] Baggio | Landau-Type Theory of Planar Crystal Plasticity[END_REF].C. R. Physique, 0000, 1, n o 0, 000-000

Acknowledgments

The first two authors contributed to this paper equally. All authors acknowledge helpful discussions with S. Conti, M. Jędrychowski, N. Greneche, J. Weiss, P. Zhang and S. Queyreau. The work was partially supported by the grants ANR-17-CE08-0047 (O. U. S. and L. T.) and ANR-19-CE08-0010-01 (O. U. S). The possibility to access computational resources (MAGI) provided by Sorbonne Université Paris Nord is highly appreciated.