Crystal elasto-plasticity on the Poincar\'e half-plane

TL;DR

This paper models 2D crystal plasticity using a nonlinear variational approach on the Poincaré half-plane, revealing how lattice symmetry influences defect dynamics and plastic flow without auxiliary hypotheses.

Contribution

It introduces a symmetry-invariant strain energy framework on the Poincaré half-plane for 2D crystals, linking lattice symmetry to defect interactions and plasticity mechanisms.

Findings

Symmetry constrains the energy landscape and defect behavior.

Plastic flow occurs via basin-hopping strain activity.

Numerical simulations show symmetry's impact on defect patterns.

Abstract

We explore the nonlinear variational modelling of two-dimensional (2D) crystal plasticity based on strain energies which are invariant under the full symmetry group of 2D lattices. We use a natural parameterization of strain space via the upper complex Poincar\'e half-plane. This transparently displays the constraints imposed by lattice symmetry on the energy landscape. Quasi-static energy minimization naturally induces bursty plastic flow and shape change in the crystal due to the underlying coordinated basin-hopping local strain activity. This is mediated by the nucleation, interaction, and annihilation of lattice defects occurring with no need for auxiliary hypotheses. Numerical simulations highlight the marked effect of symmetry on all these processes. The kinematical atlas induced by symmetry on strain space elucidates how the arrangement of the energy extremals and the possible bifurcations of the strain-jump paths affect the plastification mechanisms and defect-pattern complexity in the lattice.

Figures

Equations

$$\mathbf{C}={\mathbf{F}}^T\mathbf{F}={\mathbf{C}}^T>0,\mathbf{C}=2\mathbf{E}+\mathbf{I},$$

$$\bar{\mathbf{C}}=(\det \mathbf{C}{)}^{-1/2}\mathbf{C},$$

$$\sigma =\sigma (\mathbf{C}),$$

$$\sigma (\mathbf{C})=\sigma ({\mathbf{G}}^T\mathbf{CG}),$$

$$\sigma (\mathbf{C})={\sigma }_{\text{v}}(\det \mathbf{C})+{\sigma }_{\text{d}}(\bar{\mathbf{C}}).$$

$$\hat{z}(\mathbf{C})=\frac{C_{12}}{C_{11}}+i\frac{\sqrt{\det \mathbf{C}}}{C_{11}}\in \mathbb{H},$$

$$\hat{z}(\bar{\mathbf{C}})={\bar{C}}_{11}^{-1}({\bar{C}}_{12}+i).$$

$$\tilde{z}=\{\begin{array}{ll}\frac{m_{22}z+m_{12}}{m_{21}z+m_{11}}& \text{if }\det \mathbf{m}=1\\ \text{vskip}12.0ptplus4.0ptminus4.0pt& \\ \frac{-m_{22}\bar{z}+m_{12}}{-m_{21}\bar{z}+m_{11}}& \text{if }\det \mathbf{m}=-1.\end{array}$$

$$\mathcal{D}=\{z\in \mathbb{H}:|z|\ge 1,0\le \text{Re}(z)\le \frac{1}{2}\}.$$

$${\sigma}_{\text{d}}(\mathbf{\bar{C}})={\sigma}_{\text{d}}\big{(}J(\hat{z}(\mathbf{\bar{C}}))\big{)}.$$

$${\sigma }_{\text{d}}={\sigma }_{\text{d}}(J_0(z)),J_0(z)=|J(z)-J(z_0)|,$$

$${\sigma}_{\text{d}}\big{(}\mathbf{\bar{C}}\big{)}=\mu J_{0}(z)^{2}=\mu|J(z)-J(z_{0})|^{2}\quad\text{if}\quad J^{\prime}(z_{0})\neq 0,$$

$$\sigma_{\text{d,plast}}\big{(}\mathbf{\bar{C}}\big{)}=\mu J_{0}(z)^{2/\kappa(z_{0})}=\mu\big{|}J(z)-J(z_{0})\big{|}^{2/\kappa(z_{0})},\quad z=\hat{z}(\mathbf{\bar{C}}),$$

$$\sigma_{\text{sq,aniso}}\big{(}\mathbf{\bar{C}}\big{)}=\mu_{1}\;\text{Re}^{2}\sqrt{J(z)-1}+\mu_{2}\;\text{Im}^{2}\sqrt{J(z)-1}\,,$$

$$j(\tau )=\frac{1}{q}+744+196884q+21493760q^2+864299970q^3+20245856256q^4+\dots ,$$

$${\sigma }_{\text{v}}(\det \mathbf{C})=\lambda (\det \mathbf{C}-\log \det \mathbf{C}),$$

$$\sigma_{\text{d,sq}}=\mu\big{|}J(z)-1\big{|}.$$

$$i\to \zeta \to i+1\text{vs.}i\to i+1,$$

$$\sigma_{\text{d,hex}}=\mu\big{|}J(z)\big{|}^{2/3}.$$

$$\mathbf{H}\nabla \sigma -\xi \nabla \sigma =0,\xi \in \mathbb{R},\text{and}\mathbf{He}\cdot \mathbf{e}\ge 0,$$

Full text

1 Introduction

Crystal elasto-plasticity, especially at the microscales, represents a peculiar meeting point between coexisting effects such as solids elasticity and plastic flow with dislocation-driven intermittency, patterning, and hardening, see Miguel et al., (2001), Kocks and Mecking, (2003), Dimiduk et al., (2006), Csikor et al., (2007), Fressengas et al., (2009), Uchic, (2009), Irastorza-Landa et al., (2016), Papanikolaou et al., (2017), Sills et al., (2018). Within the framework of continuum mechanics, linear elasto-plasticity satisfactorily explains most of crystals’ behavior at the macroscopic scales. As the elastic and plastic distortions are largely independent at these scales, a number of ad-hoc modelling assumptions need to be introduced, including yield conditions and plastic flow rules (Gurtin et al.,, 2010; Chung and Lee,, 2018).

At the microscopic scales, however, elasto-plastic effects are inherently coupled, and there has been extensive research on unified modeling approaches. Phase-field models have been successful in reproducing the evolution of the ground-state configurations which is at the basis of crystal plasticity, especially in the more recent nonlinear implementations by Levitas, (2013, 2018) or Biscari et al., (2015), Vattré and Denoual, (2016), Xie et al., (2019). Atomistic and discrete dislocation dynamics (DDD) models, as in Chan et al., (2010), Alava et al., (2014), Healy and Ackland, (2014), Xiang et al., (2017), Záležák et al., (2017), Zepeda-Ruiz et al., (2017), Sills et al., (2018), have also greatly expanded our knowledge of crystal microplasticity phenomena.

However, on the one hand both phase-field and DDD methods need auxiliary hypotheses on plastic flow and/or defect nucleation and interaction. Atomistic-type modelling, on the other hand, while free from such drawbacks, concerns only short time scales, and is also not efficient in elucidating how crystal symmetry and kinematics influence defect and microstructure formation and evolution in the distorted solids. For instance, the fundamental role of kinematic compatibility in the mechanical behavior of crystalline materials may be appreciated only through an in-depth analysis of the deformation-gradient map representing the lattice distortion (Song et al.,, 2013; Biscari et al.,, 2015; James,, 2018; Feng et al.,, 2019).

To sidestep the above-mentioned issues in the modelling of crystal mechanics, in this paper we pursue an approach to crystal elasto-plasticity based on a nonlinear continuum-type variational framework explicitly informed by global crystal symmetry. The behavior of the deforming material is understood in terms of energy-minimizing strain fields on a suitable energy landscape, whose topography is determined by the full invariance of the underlying lattice. This makes non-convex the strain energy density, with a symmetry-prescribed countable infinity of minimizers. The associated valley floors (Hoffman et al.,, 1986) give the preferential avenues on which occurs the coordinated basin-hopping of the local strains, strongly directing the progress of the body’s deformation field on the bumpy landscape of the total energy functional. In our lattice-based, and thus discrete, systems such evolution is regularized (Baggio et al.,, 2019), and bursty material response emerges as a discrete sequence of local minimizers is visited under the external driving. In this context crystal plasticity results from the formation and propagation of crystal defects, such as dislocations, with no need for extra assumptions on their nucleation, annihilation and movement, as these derive from concentrated slip processes which allow large strain relaxation in the crystal while seeking to locally preserve the lattice structure.

Inertia and the possible presence of oscillatory external forces may influence the bursty dynamics of crystals, for instance affecting the scaling behaviour and shape of plastic avalanches (Truskinovsky and Vainchtein,, 2008; Salerno and Robbins,, 2013; Weiss,, 2019; McFaul et al.,, 2019). Such effects may become relevant or even dominate in particular situations, especially when the characteristic time of the driving forces approaches the avalanche durations (McFaul et al.,, 2019), or when the dislocations are induced to move at high velocities (Tang,, 2018). However, when the driving and the relaxation time scales are sufficiently separated, these effects may be neglected, and it is possible to determine the equilibria through local energy-minimization as a limit case of over-damped dynamics driving the system towards some nearest optimal configurations on a very slowly varying energy landscape. Also the viscous contributions that may in general affect the way and the characteristic times in which local equilibrium is approached, may be expected to have a less significant role in these quasi-static processes where the system visits a sequence of equilibria under the very slow driving. The role of time-scale separation and the possibility of using energy minimization for the slower relaxation phenomena have been discussed for instance in relation to slip concentration by Freddi and Royer-Carfagni, (2016), or in shear band propagation in the presence of fracture by Arriaga and Waisman, (2017).

In more detail, we model crystal elasto-plastic behavior based on the original proposal by Ericksen, (1977, 1980) concerning the material symmetry of a crystalline substance. Accordingly, we assume the stain energy density $\sigma$ to be invariant under all the deformations which map the underlying lattice onto itself, see also Folkins, (1991), Parry, (1998). These lattice-invariant distortions thus dictate, in the space of strain tensors, the location of infinitely-many relaxed states for the crystal. Explicitly, these are given by the elements of (a conjugate of) the infinite and discrete group $\text{GL}(n,\mathbb{Z})$ collecting all the invertible $n\times n$ matrices with integral entries, where $n=2,3,$ is the dimension of the lattice under study. For brevity, we refer to this property as GL-invariance.

Extensive research has successfully used the ensuing variational models. Special attention has been given to reversible martensitic transformations, with the aim of better understanding and enhancing the performance of shape-memory alloys (Pitteri and Zanzotto,, 2002; Bhattacharya,, 2004; Bhattacharya and James,, 2005; Song et al.,, 2013; James,, 2015, 2018). In these cases the associated finite deformations are largely confined to suitable ‘Ericksen-Pitteri neighborhoods’ (EPNs) in strain space, whereon GL-symmetry reduces to point-group symmetry (Ericksen,, 1980; Pitteri,, 1984; Pitteri and Zanzotto,, 2002; Conti and Zanzotto,, 2004). This much reduces the difficulties related to the full GL-invariance of the theory. GL-symmetry has also been used to study reconstructive structural transformations, in which strains may reach or overcome the EPN bundaries, producing lattice-defects and plasticity phenomena (Pitteri and Zanzotto,, 2002; Conti and Zanzotto,, 2004; Bhattacharya et al.,, 2004; Caspersen et al.,, 2004; Lew at al,, 2006; Pérez-Reche et al.,, 2007; Denoual et al.,, 2010; Vattré and Denoual,, 2016; Pérez-Reche et al.,, 2016; Pérez-Reche,, 2017). Full lattice invariance has furthermore been used for instance in the study of slip and twinning in helical structures (Feng et al.,, 2019). In general, GL-invariance plays a key role in the determination of lattice twinning mechanisms (Pitteri and Zanzotto,, 2002; Bhattacharya,, 2004), in relation to both structural transformations and plasticity in crystals. Depending on circumstance, dislocation-based irreversible plastic phenomena may be assisted also by multiple-cell twinning deformations, related to the fact that crystals can de-couple internal degrees of freedom behaving as ‘multi-lattices’ (unions of congruent Bravais lattices), rather than always as a single Bravais lattice. The details of the more complex kinematics and energetics needed for the extension of the present GL-approach to encompass also these deformation modes can be found in Pitteri and Zanzotto, (2002), see also James, (2018). Sun et al., (2017) give the corresponding treatment of twinning proliferation in magnesium, while Lew at al, (2006) use this framework to establish a bridge between GL-invariant elasto-plasticity and quantum mechanical approaches to strain-energy computations. Thus, the proposed line of research based on GL-energetics provides an attractive avenue also for the investigation of actual crystal plasticity. Besides the mentioned benefits afforded by the continuum approach, a further main point of interest here is the possibility of capturing in an intrinsic way the differences reported in experiments on the plastic behavior of crystals with different symmetry, orientation, or loading (Brinckmann et al.,, 2008; Uchic,, 2009; Weiss et al.,, 2015; Zhang et al.,, 2017; Papanikolaou et al.,, 2017; Cui et al.,, 2017; Alkan et al.,, 2018; Cho et al.,, 2018; Sparks and Mass,, 2018, 2019).

In the present work we advance this theoretical approach, evidencing further basic effects GL-elasto-plasticity in 2D. While interesting per se, this also provides a necessary road-map for the much more complex study of the 3D case. For recent results on 3D crystal plasticity via phase-field modelling under full lattice symmetry, see Biscari et al., (2015).

As in Folkins, (1991), Parry, (1998), Baggio et al., (2019), here we use a natural parameterization of 2D strain space by means of the upper complex Poincaré half-plane $\mathbb{H}$, one of the best known models of the hyperbolic plane (Kilford,, 2008; Terras,, 2013). On $\mathbb{H}$ the well-known Dedekind tessellation (Mumford et al,, 2002; Ye et al.,, 2005) transparently displays the action of GL-symmetry on strain space, and thus on a crystal’s energy landscape. By elaborating on the ideas by Parry, (1998), we write explicitly 2D-lattice strain energies which automatically comply with GL-invariance by using suitable modular forms, a well known class of complex functions considered in various branches of mathematics and physics (Schoeneberg,, 1975; Apostol,, 1976; Stillwell,, 2001; Mumford et al,, 2002; Kirsten and Williams,, 2010). By design, in the present study of crystal plasticity we restrict ourselves to GL-energies with a minimum complexity of the associated topography. We thus consider a family of potentials with a single minimizer (periodically replicated through GL-periodicity) and a very simple explicit expression, see also Baggio et al., (2019). This is motivated by the aim of obtaining plastic-flow results without possibly spurious effects due to the presence of other metastable lattice configurations. Similar methods can be used to write also GL-invariant potentials with two or more minimizers coexisting in the fundamental domain of GL-periodicity (see Patriarca, (2019) for preliminary results), thus allowing for the modelling of crystal plasticity and structural phase transformations, as in Conti and Zanzotto, (2004).

The performed numerical simulations highlight the influence of GL-constrained energy topography on plastification and defect-generation mechanisms. In particular it is evidenced the key role played by the symmetry-constrained arrangement of the ground states and of the possible bifurcations of the strain-jump paths.

This paper is organized as follows. In Section 2 we give preliminaries on the GL-invariant energies for crystalline materials. In Section 3 we present the Poincaré half-plane $\mathbb{H}$ as a parameterization of strain space in 2D, and we introduce the Klein invariant $J$, by means of which we construct a class of simple GL-invariant energies in Section 4. We use the latter for the numerical simulations of crystal plasticity presented and discussed in Section 5. In Section 6 we stress the role on the plastification mechanisms of the GL-constrained networks built by means of the energy extremals and valley floors.

2 Strain energies for crystalline materials

In a two-dimensional (2D) setting, we consider a body whose points $\mathbf{X}$ have coordinates $X_{i}$ ($i=1,2$) with respect to an orthonormal basis. The deformations are one-to-one maps $\mathbf{x}=\mathbf{x}(\mathbf{X})$, where the $x_{i}$ ($i=1,2$) identify current point positions. The deformation gradient $\mathbf{F}=\bm{\nabla}\mathbf{x}$ is a $2\times 2$ matrix with elements $F_{ij}=\partial x_{i}/\partial X_{j}$, and the associated symmetric positive-definite Cauchy-Green tensor is

[TABLE]

where $\mathbf{E}$ is the nonlinear strain (Pitteri and Zanzotto,, 2002; Gurtin et al.,, 2010). The unimodular tensors

[TABLE]

with $\det\mathbf{\bar{C}}=1$, describe a 2D hyperboloid within the 3-dimensional cone of $2\times 2$ symmetric positive-definite tensors $\mathbf{C}$ in (1).

The total strain-energy functional is the integral of a strain-energy density over the reference configuration of the body. The energy density $\sigma$ is a smooth Galilean-invariant scalar function, which depends on $\mathbf{F}$ only through $\mathbf{C}$, i.e.

[TABLE]

with material symmetry requiring that

[TABLE]

for any $\mathbf{C}=\mathbf{C}^{T}>0$ and any deformation $\mathbf{G}\in\cal G$. Here $\cal G$ is a suitable group contributing to the characterization of the material response. As mentioned in the Introduction, for 2D crystalline substances we follow Ericksen’s earlier proposal and assume $\cal G$ to be conjugate, via the choice of reference basis, to the group $\text{GL}(2,\mathbb{Z})$ describing the global symmetry of 2D Bravais lattices (Ericksen,, 1980; Michel,, 2001; Pitteri and Zanzotto,, 2002).

The strain energy is split into a convex volumetric part $\sigma_{\text{v}}$ penalizing the departure of $\det\mathbf{C}$ from 1, plus a distortive term ${\sigma}_{\text{d}}$ which only depends on $\mathbf{\bar{C}}$:

[TABLE]

By GL-invariance, ${\sigma}_{\text{d}}$ is non-convex with countably-many local minimizers, whose arrangement in strain space we discuss below (see Fig. 1).

3 Parameterization of 2D distortions on the Poincaré half-plane

The Poincaré (upper complex) half-plane is the set $\mathbb{H}=\{x+iy\in\mathbb{C}$, $y>0\}$, endowed with the metric (ds)^{2}=\big{(}(dx)^{2}+(dy)^{2}\big{)}/y^{2}, see Kilford, (2008); Terras, (2013). A complex representation of the cone of tensors $\mathbf{C}=\mathbf{C}^{T}>0$ is provided by associating to any $\mathbf{C}$ the complex number as in Folkins, (1991), Parry, (1998):

[TABLE]

where $C_{ij}$ are the matrix elements of $\mathbf{C}$ in the given orthonormal basis. From Eq. (6) we obtain $\hat{z}(\mathbf{C})=\hat{z}(\mathbf{C}^{\prime})$ if and only if $\mathbf{C}$ and $\mathbf{C}^{\prime}$ are proportional. This means that the 2D $\mathbf{\bar{C}}$-hyperboloid of unimodular strains in Eq. (2) is smoothly mapped one-to-one onto $\mathbb{H}$ by:

[TABLE]

The local body deformation, described by the four entries $F_{ij}$ of the deformation gradient $\mathbf{F}$, is thus equivalently accounted for by using $\det\mathbf{F}$, the angle $\theta$ in the polar decomposition of $\mathbf{F}$, plus two independent variables in $\mathbf{\bar{C}}$ (or the corresponding $\hat{z}(\mathbf{\bar{C}})$ from Eq. (7)) which enter ${\sigma}_{\text{d}}$ in Eq. (5).

The identification in Eq. (7) is very useful in the present approach to crystal mechanics because the material-symmetry maps $\mathbf{C}\mapsto\mathbf{G}^{T}\mathbf{C}\mathbf{G}$, for $\mathbf{G}\in\cal G$, correspond via (7) to the natural action on $\mathbb{H}$ by the classical fractional linear (Moebius) isometries of $\mathbb{H}$ with integral entries, augmented by the map $z\mapsto-\bar{z}$, where $\bar{z}$ is the conjugate to $z\in\mathbb{C}$ (Folkins,, 1991; Parry,, 1998). Explicitly, if $\mathbf{m}\in\text{GL}(2,\mathbb{Z})$ and $\tilde{\mathbf{C}}=\mathbf{m}^{T}\mathbf{C}\mathbf{m}$, the corresponding map $z\mapsto\tilde{z}$ on $\mathbb{H}$ via (7) is given by

[TABLE]

A very effective geometrical representation of this action is given by the $\text{GL}(2,\mathbb{Z})$-generated Dedekind tessellation of $\mathbb{H}$ shown in Fig. 1 (Mumford et al,, 2002; Ye et al.,, 2005). It contains all the $\text{GL}(2,\mathbb{Z})$-related mutually congruent (in the sense of the hyperbolic metric) copies of the fundamental domain

[TABLE]

The points in the interior of $\mathcal{D}$ correspond through Eq. (7) to strain tensors (and thus lattices) with trivial symmetry and thus to oblique lattices, while the points on the boundary $\partial\mathcal{D}$ of $\mathcal{D}$ correspond to metrics possessing nontrivial symmetries. The latter include (see Fig. 1) rectangular lattices, represented by boundary points on the imaginary axis; fat-rhombic lattices (whose inner angles are greater than $\frac{\pi}{3}$), represented by boundary points with $\lvert z\rvert=1$; and skinny-rhombic lattices, represented by boundary points with $\text{Re}(z)=1/2$. We refer to Parry, (1998), Michel, (2001), Pitteri and Zanzotto, (2002), Conti and Zanzotto, (2004) for a detailed description of the global and point symmetry of the Bravais lattice types in 2D, and of the $\text{GL}(2,\mathbb{Z})$ fundamental domain in $\mathbf{\bar{C}}$-space.

The identification (7) of the $\mathbf{\bar{C}}$-hyperboloid with $\mathbb{H}$ leads to the remarkable observation by Folkins, (1991), Parry, (1998), that the GL-invariant strain energies ${\sigma}_{\text{d}}$ for 2D crystalline materials in (4)-(5) are closely related to the modular functions on $\mathbb{H}$. Particularly useful for us is the existence (Schoeneberg,, 1975; Apostol,, 1976) of a unique complex function $J$, also known as the Klein invariant, holomorphic on $\mathbb{H}$ with the following properties: $J$ is periodic under the action (8)1 of the ’modular group’ $\text{SL}(2,\mathbb{Z})$ (which denotes the subgroup collecting the positive-determinant elements of $\text{GL}(2,\mathbb{Z})$, see Ye et al., (2005)); $J$ is one-to-one between the fundamental domain $\cal D$ and the complex half-plane $\mathbb{C}^{+}$ with non-negative imaginary part; furthermore, $J$ satisfies $J(z)=\overline{J(-\bar{z})}$ for any $z$, it assumes values $J(i)=1$ and $J(\rho)=0$, and it diverges when $\text{Im\,}z\to+\infty$, with a simple pole at infinity. Its only stationary points are $z=i$ and $z=\rho$ (i.e. $J^{\prime}(i)=0$ and $J^{\prime}(\rho)=0$), with also $J^{\prime\prime}(\rho)=0$. The symmetry requirements imply that $J$ assumes real values at the boundary $\partial\cal D$ of $\cal D$, and that $J$ diverges also when $z$ approaches any rational point on the real axis, see also Schoeneberg, (1975), Apostol, (1976), Parry, (1998), and Eq. (15) below.

4 A class of simplest GL-invariant strain-energies for 2D crystal elasto-plasticity

Based on the strain parameterization (7) and the above properties of $J(z)$, Parry, (1998) proposed that potentials satisfying (4)-(5) be written as suitable functions of $J$:

[TABLE]

As mentioned in the Introduction, depending on the specific properties and number of minimizers and other extremals in the domain $\cal D$ of Eq. (9), energies of this form can describe plasticity and phase transformations in 2D crystals, as well as their coupling. Our present aim is deriving a family of simplest potentials in (4)-(5)-(10) by using functions ${\sigma}_{\text{d}}(J)$ in (10) which have a single minimizer $z_{0}$ in $\cal D$, see also Baggio et al., (2019). A single minimizer, GL-related to $z_{0}$, is therefore also on each tile in the Dedekind tessellation of Fig. 1.

Given any energy-minimizing $\mathbf{\bar{C}}_{0}$, we set $z_{0}=\hat{z}(\mathbf{\bar{C}}_{0})$. Then $J_{0}(z)=|J(z)-J(z_{0})|$ provides a basic measure of the distance between any strain and the ground state, via their counterparts $z\in\mathbb{H}$. By considering in (10) functions

[TABLE]

it is possible to obtain a simple class of smooth potentials whose sole minimizers are $\mathbf{\bar{C}}_{0}$ and all its GL-related copies when $z_{0}\in\partial\mathcal{D}$111As remarked in Schoeneberg, (1975), Apostol, (1976), $J(z)$ is SL-invariant under (8)1. However, a function ${\sigma}_{\text{d}}(J_{0}(z))$ in (11) is in fact GL-invariant whenever the minimizer $z_{0}=\hat{z}(\mathbf{\bar{C}}_{0})$ is on the boundary of $\cal D$, as in this case $J(z_{0})$ is real. Ground states lying on $\partial\cal D$ correspond to lattices with non-trivial symmetries as in square, hexagonal, rhombic, or rectangular 2D Bravais lattices. Indeed, in order to have GL-invariance, strain functions as in (11) must depend on (Im $J(z))^{2}$, while they can depend arbitrarily on Re $J(z)$. If $z_{0}$ has non-trivial symmetry, any function of $J_{0}(z)$ automatically depends on (Im $J(z))^{2}$, and does exhibit full GL-symmetry. If $z_{0}$ belongs to the interior of $\cal D$, thus corresponding to an oblique 2D-lattice ground state, potentials (11) depending on $J_{0}$ would be SL- but not GL-invariant. In this case we can instead consider in (11) functions for instance of J_{0}^{*}(z)=\big{(}[\text{Re}\,J(z)-\text{Re}\,J(z_{0})]^{2}+[(\text{Im}\,J(z))^{2}-(\text{Im}\,J(z_{0}))^{2}]^{2}\big{)}^{1/2}, which gives an alternative GL-invariant definition of the distance from the ground state. To agree with standard linear elasticity when the strained lattices are sufficiently close to the well bottoms, the Taylor expansion of the functions in (11) must have positive-definite quadratic behavior close to $z_{0}$, i.e. $\sigma_{\text{d}}(J_{0}(z))\approx(z-z_{0})^{2}$ for $z\approx z_{0}\in\mathbb{H}$. Since $J_{0}(z)\approx|J^{\prime}(z_{0})||z-z_{0}|$ as $z\approx z_{0}$, when $J^{\prime}(z_{0})\neq 0$ the simplest potential satisfying our requirements is

[TABLE]

where $\mu>0$ is an elastic modulus, and $z=\hat{z}(\mathbf{\bar{C}})$. As $J^{\prime}(z_{0})\neq 0$ for all $z_{0}$ different from the corner points $i$ or $\rho$ of $\cal D$ in Fig. 1, the potential in (12) is thus the desired one for all $z_{0}\in\partial\cal D$, $z_{0}\neq i,\rho$.

The particularly relevant cases $z_{0}=i$ (square) and $z_{0}=\rho$ (hexagonal), which represent the two maximally-symmetric lattices in 2D, need care because they are the only stationary points of $J$, so (12)2 does not hold. The Klein invariant indeed satisfies $J(z)=1+O(z-i)^{2}$ as $z\to i$, and $J(z)=O(z-\rho)^{3}$ as $z\to\rho$, because $J^{\prime}(i)=0=J^{\prime}(\rho)=J^{\prime\prime}(\rho)=0$. To ensure a correct linear-elastic behavior we must then consider suitably modified potentials in (11)-(12), i.e. $(|J(z)-J(i)|^{1/2})^{2}=|J(z)-1|$ for square lattices ($z_{0}=i$), and $(|J(z)-J(\rho)|^{1/3})^{2}=|J(z)|^{2/3}$ for hexagonal lattices ($z_{0}=\rho$) .

Summarizing, a unified expression for the simplest elasto-plastic GL-potentials for 2D crystals with square, hexagonal, rhombic or rectangular ground state $z_{0}=\hat{z}(\mathbf{\bar{C}}_{0})\in\partial\mathcal{D}$ is

[TABLE]

where $\mu>0$ and $\kappa(z_{0})$ is the order of zero of $J(z)-J(z_{0})$ near $z_{0}$ (in particular, $\kappa(z_{0})=1$ if $J^{\prime}(z_{0})\neq 0$, i.e. for all $z_{0}\in\partial\mathcal{D}\setminus\{i,\rho\}$).

We notice that writing ${\sigma}_{\text{d}}$ in (13) in terms of the above $J$-based modular order parameters further develops the notion of transcendental order parameters considered in the extended Landau theory for crystalline materials examined in Dmitriev et al., (1988), Horovitz et al., (1989), as $J_{0}(z)$ endows the strain energy density $\sigma$ in (5) with the full invariance of the underlying (2D) lattice. The strain-energies $J_{0}(z)$ in (13) are a simplest class of fully GL-invariant regular functions possessing a unique minimizer $z_{0}$ in the fundamental domain $\cal D$. Given this, $J_{0}(z)$ then has in $\cal D$ at most three stationary points, i.e. the imposed minimizer $z_{0}$, and the two further points $i$ and $\rho$ where any regular GL-invariant function, in particular as in (13), is necessarily stationary due to GL-symmetry. The two cases where either $z_{0}=i$ (energy with square ground states) or $z_{0}=\rho$ (energy with hexagonal ground states) are explicitly considered below. Besides providing simple strain-energies exhibiting the full lattice symmetry and a tunable ground state, the functions (13) may be used as building blocks within the wider family of energies (10). They can for instance provide explicit expressions, both analytic and GL-invariant, for potentials in structural phase transformations and their associated elasto-plastic effects, as in Patriarca, (2019). They can also be used to describe the generalized stacking fault energies ($\gamma$-surfaces) derived through ab-initio (Bakst et al.,, 2018; Kamimura et al.,, 2018) or atomistic phase field (Qiu et al.,, 2019) computations. Polynomial combinations of $\sigma_{\text{d,plast}}$ in (13) give the simplest class of functions for these purposes, with coefficients which can be determined from the values of the elastic constants, the Peierls stresses, the properties of possible supplementary stable states, or other relevant experimentally determined mechanical quantities.

We mention explicitly that the potentials in (13) describe lattices with isotropic elastic moduli. Anisotropic elasticities may be introduced (except for the hexagonal case, where they are forbidden by symmetry) by assuming a suitable separate dependence of the potentials on the real and imaginary parts of $J(z)-J(z_{0})$ rather than on its modulus $J_{0}(z)=|J(z)-J(z_{0})|$ as in (11)-(13). For example, an energy exhibiting anisotropic square elasticities is obtained by writing:

[TABLE]

with $\mu_{1},\mu_{2}$ positive moduli.

5 Plastification under shearing

We illustrate some relevant features of the model through the results of two numerical simulations of plastic-flow initiation in quasi-static conditions, for a square and a hexagonal lattice. To emphasize the basic effects of symmetry in these plastifying crystals, we consider the simplest case of loading each body along a primary-shear direction (i.e. a densest crystallographic line) in the lattice.

The body is square-shaped, and contains an initially homogeneous material in which the underlying lattice is defect-free. The bottom side is aligned with a dense (primary-shear) lattice direction, in turn coinciding with one of the given orthonormal basis vectors. The strain-controlled shear boundary conditions are imposed to the top side with the bottom side kept fixed, while the remaining two sides are left stress-free. The boundary shear parameter is $\gamma\in\mathbb{R}$, with $\gamma=0$ in the ground state $z_{0}$ for both the square ($z_{0}=i$) and the hexagonal ($z_{0}=\rho$) case. For each $\gamma$ we search for the local minimizers of the body’s total strain-energy functional computed through the density (5)-(13), and complying with the imposed boundary conditions.

The minimization algorithm follows a preconditioned energy-gradient descent, implemented through a continuous FreeFEM code (Hecht,, 2012) using an unstructured mesh for a total of about $3\times 10^{4}$ degrees of freedom. At each step, the initial guess coincides with the metastable state associated with the previous value of the imposed loading. The algorithm then looks for metastable states which are close (in the sense of the preconditioned gradient-descent) to the previous equilibrium. The FEM discretization introduces in our non-convex problem a regularizing length scale related to lattice discreteness. Therefore, the mesh is fixed, and the solutions are the piecewise-differentiable functions with constant gradient over FEM cells and possible gradient jumps over mesh sides.

A well-known property of $j(z)\equiv 1728J(z)$ is that it has all integers in its Fourier-expansion coefficients (Schoeneberg,, 1975; Apostol,, 1976):

[TABLE]

where $q=\exp(2\pi i\tau)$. In particular, we used the first 25 terms (Van Wijngaarden,, 1953) in (15) for our minimization algorithm. The volumetric energy in (5) has been chosen of the form

[TABLE]

with $\lambda/\mu\sim 30$ to enforce quasi-incompressibility.

5.1 Shearing of a square lattice

In the first simulation (see Fig. 2) the energy ${\sigma}_{\text{d}}$ in (5) is given by (13) with $z_{0}=i$, $J(i)=1$, and $\kappa(z_{0})=2$, i.e.:

[TABLE]

The imposed primary-shear path $i\rightarrow 1+i$ is the straight horizontal dashed-blue line in Figs. 2(a)-(b)-(c), with parameter $\gamma$ such that $\gamma=0$ in the ground state $z_{0}=i$, and $\gamma=1$ in the neighboring equivalent sheared-square configuration $i+1$ (see Figs. 1 and 2(a)).

Figs. 2(d)-(e)-(f) show three snapshots of the resulting $\gamma$-dependent plastification field. Figs. 2(a)-(b)-(c) show the clustering of $\mathbf{\bar{C}}$-values, which, via (7) is represented by a $\gamma$-evolving cloud of points on the Dedekind tessellation of $\mathbb{H}$. This distribution of local-strain density gives a very informative picture of the main properties of the plastification field, see also Balandraud et al., (2015), Baggio et al., (2019). The Supplementary Video V1 shows the entire simulation, reporting also the $\gamma$-dependent histograms for the local rotation angle $\theta$ and of $\det\mathbf{F}$. Figs. 2(a)-(b)-(c) show the defects mediating plastic flow. They are heuristically identified through the centrosymmetry parameter used by the visualization software Ovito (Stukowski,, 2010), which counts the number of nearest neighbors of lattice points. In this setting dislocations emerge when nearby cells undergo concentrated slips originating from shears in GL$(2,\mathbb{Z})$ which cannot be accommodated by elastic deformations. The stress distribution of a trapped lattice defect in Figs. 2(d)-(e)-(f) shows the correct features of dislocations’ elastic field both far from and in the proximity of the core region, as in Fig. 4 in Baggio et al., (2019).

Fig. 2(g) shows that the initially defect-free lattice goes through a significant elastic load-up with a quasi-homogeneous initial deformation. As $\gamma$ increases, the strain cloud in $\mathbb{H}$ widens due to the growing strain heterogeneity originating from the unloaded body-sides.

The end of the elastic regime is marked by a first large plastification event at about $\gamma\simeq 0.17$, when the strain cloud suddenly jumps in $\mathbb{H}$ away from the initial point $i$, see Fig. 2(b). This first large plasticity burst occurs as local strains cross over to both the two symmetry-equivalent square energy-basins closest to $i$ in the direction of the imposed boundary condition, i.e. $i+1$ and $\zeta$.

This complex plasticity mechanism originates from the way global GL-symmetry locally moulds the energy landscape, and how this interacts with the driving boundary condition. Both the square points $i+1$ and $\zeta$ are involved because the optimal barrier-crossing path from $i$ to $i+1$ goes through the hexagonal point $\rho=e^{i\pi/3}$, which is a degenerate monkey saddle.222Since as noted earlier $J^{\prime}(\rho)=J^{\prime\prime}(\rho)=0$, for the square energy function $\sigma_{\text{d,sq}}$ in (17) there are degenerate monkey saddles in $\rho$ and in all other hexagonal points of $\mathbb{H}$, with null Hessian and third-order terms, up to a suitable rotation, of the form $x^{3}-3xy^{2}$. There is thus in $\rho$ a bifurcation of the optimal square-to-square barrier-crossing path $i\rightarrow i+1$ also to the optimal path and $i\rightarrow\zeta$ going to the other nearby equivalent square point $\zeta$ (see the dashed-red fat-rhombic curves in Figs. 2(a)-(b)-(c)).

Due to this bifurcation, plastic-flow initiation does not occur through the activation of the primary lattice shear $i\rightarrow i+1$ (on the [0,1] lattice line) pertaining to the external driving $\gamma$ (dashed-blue line), but rather through a strain avalanche involving mainly the other primary shear path $i\rightarrow\zeta$ (dotted-blue line in Fig. 2(a)), which is on the [1,0] line and is symmetry-related to $i\rightarrow i+1$ (note that these two shear paths coincide in the linear approximation at $i$). As a consequence, although the body is being loaded in the horizontal shear direction $i\rightarrow i+1$, we observe the formation of the vertical slip bands in Fig. 2(e) as a first complex event marking the onset of plastification in the square lattice (see also the Supplementary Video V1). The width and shape of the associated strain cloud in Fig. 2(b) actually indicates that not only are both the primary shears $i\rightarrow i+1$ and $i\rightarrow\zeta$ becoming simultaneously active as the imposed boundary condition follows the path $i\rightarrow i+1$, but that plastic flow indeed takes place with the concurrent activation of many other local deformation paths, with the elastic stabilization of strain values also belonging to the non-convexity region near the degenerate saddle in $\rho$ (see Fig. 2(b)).333The locally-destabilizing shear directions originating from the loss of positive-definiteness of the acoustic tensor examined by Baggio et al., (2019) further contribute to specifying the bifurcation behavior of the energy-relaxing strain field. We notice there the creation of higher strain inhomogeneities and a finer spatial distribution of defects in the body already at the first nucleation (see Fig. 5 in Baggio et al., (2019)), due to the imposition of hard boundary conditions. By contrast, in the present study we observe the creation of a broad-band defect pattern upon the first strain burst (Fig. 2(a)-(b)). This is due the less stringent boundary constraints presently considered, involving free body sides.

A large stress drop corresponds to this first relaxation event, see Fig. 2(g), giving rise to the nucleation peak in the stress-strain relation, as expected in this originally clean system (Fedelich and Zanzotto,, 1992; Truskinovsky and Vainchtein,, 2004; Ding et al.,, 2012). The orange spikes in Fig. 2(g) show that after the first event a discontinuous strain activity is ever present in the body due to the coordinated basin hopping of local strain values. Correspondingly, intermittent defect nucleation and dynamics spontaneously progress in the lattice by energy minimization, with no need for extra assumptions, see also the Supplementary Video V1. Bursty plastic flow ensues, as is typically observed in the microscale plasticity of crystalline materials (Miguel et al.,, 2001; Dimiduk et al.,, 2006; Csikor et al.,, 2007; Uchic,, 2009; Papanikolaou et al.,, 2017).

Under a further shear increase the body eventually separates (for $\gamma\simeq 0.28$) with a large plastic slip into two clearly defined regions in Fig. 2(f), one almost undistorted and the other deformed in the direction of the boundary shear, at the higher shear value pertaining to the minimizer $i+1$. The width of the slip band is dictated by the shear value $\gamma$ imposed at the boundary. This large slip event is again associated to a marked stress relaxation, see the second large drop in the stress-strain curve in Fig. 2(g).

Of particular interest in Fig. 2(f) is the defect pattern produced within the slip band. This band originates as points in the upper part of the body shear to $i+1$ from both near $i$ and near $\zeta$, rather than directly from $i$, see Figs. 2(b)-(c) and the Supplementary Video V1. When neighboring lattice cells reach the same strain-energy well $i+1$ in $\mathbb{H}$ through the paths

[TABLE]

different values in their local rotation angle $\theta$ may occur (and are energetically cost-free due to Galilean invariance (3)). Such rotation differences in nearby cells experiencing the same strain $i+1$ create kinematic incompatibilities which give rise to the slip-band defects shown in Fig. 2(f), see also Biscari et al., (2015), Gao et al., (2019).

We have clear evidence here of the importance of the bifurcations on the paths followed by nearby local strains while changing basin driven on the GL-topography. These bifurcations behave as ‘disorder engines’, enhancing defect creation in the body, in the first nucleation events and all along the flow. The basic energy $\sigma_{\text{d,sq}}$ in (17) presents such effects in an archetypical way for square lattices, due to the bifurcations occurring through the degenerate monkey saddles at the hexagonal configurations.

5.2 Shearing of a hexagonal lattice

The second simulation is performed, under the same conditions as the previous one, for an initially defect-free hexagonal crystal with ground state $\rho$ (Fig. 1). The energy ${\sigma}_{\text{d}}$ in (5) is now given by (13) for $z_{0}=\rho$, $J(\rho)=0$, and $\kappa(z_{0})=3$:

[TABLE]

As in the square case, we see that after the initial elastic load-up, plastic flow initiates through a large plastification avalanche. Figs. 3(a)-(b) present a snapshot of the corresponding body shape and strain cloud in $\mathbb{H}$. The inset in Fig. 3(a) shows the strong stress relaxation associated to this first strain event, and the bursty character of the ensuing plastic flow, as indicated by the orange strain-activity spikes. Full results are reported in the Supplementary Video V2.

Unlike with the square case discussed above, plasticity in this hexagonal lattice largely follows the primary-shear path imposed by the boundary conditions on the body. This is again a direct consequence of the GL-symmetry constraints. Indeed, for the hexagonal energy $\sigma_{\text{d,hex}}$ in (19), with minimizers on the hexagonal points $\rho$, $\rho+1$, $...$, the saddles at the square points $i$, $i+1$, $\zeta$, $...$, are necessarily standard, with an indefinite non-degenerate Hessian. For instance, in Fig. 3(b) the square saddle in $i+1$ is the mountain pass for the red-dashed hexagonal-to-hexagonal optimal path $\rho\rightarrow\rho+1$ (the dashed blu line indicates the primary shear path). In general, each square saddle is here the mountain pass on a unique hexagonal-to-hexagonal optimal path, with no barrier-crossing bifurcations. We remark that in this case the slipped portion of the crystal is largely defect free, see Fig. 3(a), unlike with the square case in Fig. 2(f).

6 Networks of energy extremals and valley floors: Bethe trees and Husimi cactuses

The simulations in the previous Section extend the observations in Baggio et al., (2019) concerning the complexity of collective defect nucleation in the present plasticity framework. They especially highlight the role played in crystal plasticity by the symmetry-imposed constraints on the energy extremals and on the basin-hopping strain paths on $\mathbb{H}$, with their possible bifurcations. In particular, we see such paths are largely directed along the valley floors of the energy landscape. To better elucidate these effects, we explicitly consider the gradient-extremal curves on $\mathbb{H}$ pertaining to our energy functions. These are loci defined by the condition (20)1 that $\bm{\nabla}\sigma$ be an eigenvector of $\mathbf{H}$ (Hoffman et al.,, 1986; Quapp,, 1989; Sun and Ruedenberg,, 1993), where $\bm{\nabla}\sigma$ and $\mathbf{H}$ are respectively the gradient and the Hessian on $\mathbb{H}$ of the energy $\sigma$. Inequality (20)2 then identifies the valley floors of the energy topography, as the gradient-extremal curves occupying the energy landscape bottom:

[TABLE]

where $\mathbf{e}$ is a vector orthogonal to the energy gradient in the hyperbolic metric. The arrangement of these energy-surface features is summarized by an infinite GL-invariant network whose nodes are the energy extremals and edges are along the energy valley floors of $\sigma$ (see also related considerations by Vattré and Denoual, (2016), Gao et al., (2019)). The node-coordination values $c$ on this network, with the corresponding edge multiplicity, are determined by GL-symmetry. In particular, nodes with $c>2$ give valley-floor bifurcations. The topological and metric features of this network underpin a crystal’s plastification mechanisms. They make explicit the geometric scaffolding on which the GL-energetics largely leads, via the short and long range elastic interactions, the coordinated local strain-jump activity on $\mathbb{H}$ producing the plastification bursts, with the boundary conditions contributing to select the activated deformation paths.

In particular, for both the square and hexagonal energies (17) and (19) considered above, the valley-floor network coincides with the Bethe-like tree (de Miranda-Neto and Moraes,, 1992, 1993) shown in red in Fig. 4(a).444The fact that the valley-floor networks for both our energies(17) and (19) coincide with the Bethe tree in Fig. 4(a) whose edges are on the boundary of $\cal D$, can be derived directly from the GL-symmetry relations across each boundary portion of $\cal D$ and the monotonic growth properties of $J$ on $\partial\cal D$. A more familiar representation of this Bethe tree is obtained when considering a disc model of the hyperbolic plane, as in Conti and Zanzotto, (2004) or Patriarca, (2019). See also Harlow et al., (2011). This is characterized by 3-connected hexagonal nodes, with further 2-connected square nodes mid-way along the edges, given by the fat-rhombic valley floors of both (17) and (19). The plastification and defect-generation scenarios for square vs. hexagonal lattices in Section 5 originate from the different positions of the saddles on this Bethe tree. In the hexagonal case (19) the (standard) saddles are located at non-bifurcating nodes with $c=2$. Valley-floor bifurcations only take place here at the 3-connected minimizer (hexagonal) nodes, with local plastic strains largely confined near the primary-slip paths (see Fig. 4(b)). By contrast, for the square energy (17) the (degenerate) saddles coincide with the 3-connected bifurcating nodes of the Bethe tree. This promotes the complex deformation mechanisms involving the conjunct activation of several plastification paths at once and sequentially in the body.

It is interesting to briefly examine also the structure of the valley-floor networks associated to more generic GL-invariant strain potentials (4)-(10) than the simplest ones (17)-(19) considered above. A family of twice-differentiable GL-energies produced by Conti and Zanzotto, (2004) allows for a one-parameter unfolding of the hexagonal monkey saddles of $\sigma_{\text{d,sq}}$ in (17), see the square-hexagonal minimizers’ bifurcation diagram in Fig. 4 of that work. GL-invariant $J$-energies with an analogous square-hexagonal bifurcation are also considered in Patriarca, (2019), as a linear combinations of $\sigma_{\text{d,sq}}$ and $\sigma_{\text{d,hex}}$ in (17) and (19). As discussed in Pitteri and Zanzotto, (2002), the bifurcaton features, with the related unfolding of the possibly degenerate extremals, depend on the symmetry properties of the bifurcation points: in the present case for $i$ (square) we have a subcritical pitchfork to fat-rhombic configurations, and for $\rho$ (hexagonal) the bifurcation is transverse to three rhombic branches, see Conti and Zanzotto, (2004).

As a consequence, the hexagonal monkey saddles first unfold generically into ’monkey regions’ involving three standard saddles at fat [skinny] rhombic points in the vicinity of a hexagonal minimum [maximum] of the energy,555For a two-parameter unfolding of the monkey saddle and the corresponding phase diagram, see for instance Shtyk et al., (2017). see also Fig. 6 in Baggio et al., (2019). Plastification mechanisms and defect generation are largely unaffected by these small deviations from the degenerate monkey-saddle configuration.

On the other hand, when the unfolded extremals are skinny-rhombic saddles located away from the hexagonal maximizer, these become 2-connected (standard) mountain passes each between a suitable pair of minimizing square nodes. This is analogous to what happens in the hexagonal case with $\sigma_{\text{d,hex}}$ in (17), where only 2-connected standard-saddle nodes are present, and no valley-floor bifurcations other than at the minimizers. However, the valley-floor network topology associated to this new configuration of unfolded extremals is very different from the previous cases, as it now involves triangle-loops with vertices on the 4-connected square minimizers. The resulting graph is thus a Husimi cactus (de Miranda-Neto and Moraes,, 1992, 1993), rather than a loop-free Bethe tree, see the blue network in Fig. 4(a). This points to two distinct scenarios for plastification and defect generation in square lattices. By contrast, the bifurcation diagram shows there is a single scenario possible for GL-energies with hexagonal ground states, as already described above.

We conclude that interesting insight on the behavior of crystalline materials is gleaned from the energy-surface features summarized by the associated network of extremals and valley floors. Analogous studies of energy features are done in other fields, as in organic chemistry and mechanochemistry (Valtazanos and Ruedenberg,, 1986; E et al.,, 2005, 2007; Ess et al.,, 2008; Rehbein and Carpenter,, 2011; Quapp and Bofill,, 2016), protein folding and other biomolecules’ dynamics (Wales,, 2003; Mallamace et al.,, 2016), optics (Bociort and van Turnhout,, 2005), or superconductivity (Shtyk et al.,, 2017). The techniques for energy-topography analysis there developed will thus help shed light also on crystal mechanics phenomena (see also the graph analysis methods in Edelsbrunner, (2001)). In particular, taken together with the results by Baggio et al., (2019), the discussion above concerning different plasticity scenarios for different lattice symmetries should aid in the understanding of the complex behavior and marked dissimilarities experimentally observed in the microplasticity of materials with distinct crystal structures, including sensitivity to orientation and loading type, as well as defect-pattern formation and non-universal scaling exponents for dislocation avalanches, as referenced in the Introduction.

APPENDIX

Caption to Supplementary Video V1

Numerical simulation of the plastic flow initiation in a homogeneous square body containing an initially defect-free square lattice. The imposed loading is along a primary shear lattice direction, aligned with the parallel square-cell side and bottom body side. The shearing is imposed through the constrained horizontal sides of the body, with the two remaining sides free. The strain-energy is given by Eq. (17). See also Figs. 1-2 in the main text.

(a) Bursty deformation field in the body for increasing shear parameter $\gamma$, indicated by the moving green dot along the $\gamma$-axis in (e). Lattice points are color-coded according to (c) depending on the energy basin in the Poincaré half-plane $\mathbb{H}$ visited by the strain of the associated lattice cell during loading. Notice the defect dynamics intrinsically produced in this model by energy minimization.

(b)-(d)-(f) Intermittent evolution of the $\gamma$-dependent histograms of the deformation-gradient parameters recorded during the simulation.

(b) Evolution of the histogram (strain cloud) of the density of strain parameters on the Dedekind tessellation of $\mathbb{H}$ during loading. The histogram color-coding provides the percentage of body FEM cells with strain at each point of $\mathbb{H}$. The straight horizontal dashed-blue line between the two neighboring square configurations $i$ and $i+1$ is the image in $\mathbb{H}$ of the primary shear path imposed as boundary condition. The initial configuration is in $i$ for $\gamma=0$, while $\gamma=1$ corresponds to $i+1$. The energetically optimal barrier-crossing path from $i$ to $i+1$ is the red-dashed path passing through the hexagonal monkey saddle in $\rho=e^{i\pi/3}$. Such path bifurcates in $\rho$ also to the optimal barrier-crossing path from $i$ to $\zeta=\frac{1}{2}(i+1)$ (the blue-dotted line indicates the corresponding primary shear path from $i$ to $\zeta$). Shading shows the convexity domains around the hexagonal minimizers in $\mathbb{H}$ of the square energy function.

(c) Color-code map used in (a) for the (symmetry-related copies of the fundamental domain $\cal D$ which make up the various) energy-basins on $\mathbb{H}$ pertaining to the square strain-energy in Eq. (17) of the main text.

(d) Evolution of the histogram for the values of $\det\mathbf{F}$, indicating volumetric effects in the lattice.

(e) Bursty plastic flow. Jagged stress-strain behavior of the body (blue), and the associated spikes (orange) showing the percentage of basin-hopping strain values during loading, as $\gamma$ grows. The body response is elastic to about $\gamma\simeq 0.17$, after which bursty plastic flow begins, with a large relaxation event, continuing then for growing $\gamma$. A second substantial plastification takes place at $\gamma\simeq 0.28$, where the body separates, in a large plastic slip, into two clearly defined regions.

(f) Evolution of the histogram for the values of the angle $\theta$ in the polar decomposition of the deformation gradient $\mathbf{F}$, indicating local rotation in the lattice.

Notice in (d) and (f) the effects associated with the slip band formation. Both the local dilations and local rotations relax to almost homogeneous values after the second plastification event. In particular, the rotation-angle $\theta$ distribution is bimodal, with peaks corresponding to the (zero) rotation in the undistorted portion of the body and to the value of $\theta$ pertaining to the complete shear induced by the boundary condition, concentrated in the slip band.

Caption to Supplementary Video V2

Numerical simulation of plastic flow initiation in a homogeneous square body containing an initially defect-free hexagonal lattice. The imposed loading is along a primary shear lattice direction, aligned with a hexagonal-cell side parallel to the bottom body side. The shearing is imposed through the constrained horizontal sides of the body, with the two remaining sides free. The strain-energy is given by Eq. (19). See also Figs. 1 and 3 in the main text. All panels for this Supplementary Video V2 refer to the same variables and distributions as in the caption to Supplementary Video V1.

In this case the slip-band formation results from a single large plastification avalanche for $\gamma\simeq 0.29$, close to the value where we observe the slip-band formation also in the square lattice in Supplementary Video V1. The cells involved in the present slip event, driven by the boundary conditions, jump directly from the initial well at $z=\rho$ to the neighboring well at $z=\rho+1$, without the bifurcation activating other paths as observed in the square case of Supplementary Video V1. This limits the number of defects created in the slip process, and the resulting slip-band is largely defect-free, see also Fig. 3 of the main text.