Ericksen-Landau Modular Strain Energies for Reconstructive Phase Transformations in 2D Crystals

Abstract

By using modular functions on the upper complex half-plane, we study a class of strain energies for crystalline materials whose global invariance originates from the full symmetry group of the underlying lattice. This follows Ericksen’s suggestion which aimed at extending the Landau-type theories to encompass the behavior of crystals undergoing structural phase transformation, with twinning, microstructure formation, and possibly associated plasticity effects. Here we investigate such Ericksen-Landau strain energies for the modelling of reconstructive transformations, focusing on the prototypical case of the square-hexagonal phase change in 2D crystals. We study the bifurcation and valley-floor network of these potentials, and use one in the simulation of a quasi-static shearing test. We observe typical effects associated with the micro-mechanics of phase transformation in crystals, in particular, the bursty progress of the structural phase change, characterized by intermittent stress-relaxation through microstructure formation, mediated, in this reconstructive case, by defect nucleation and movement in the lattice.

1 Introduction

Ericksen’s early proposal [1, 2] of an infinite and discrete invariance group for a crystalline material’s strain energy aimed at expanding Landau-type variational approaches to encompass structural phase transformations and twinning in crystals, with the associated phenomena of fine microstructure, and possibly defects, forming in the lattice. Accordingly, the material invariance of the crystalline substance should reflect the global symmetry of the underlying lattice, with the strain energy density $\sigma$ invariant under all the deformations mapping the lattice onto itself. See also [3] for a similar viewpoint. This invariance dictates the location of countably-many ground states for the crystal in strain space, including those produced by the lattice-invariant shears and rotations which play a key role in twinning mechanisms when in the presence of structural phase changes, as well as in lattice-defect creation and ensuing plastification phenomena [4–11].

A wide-ranging extension of non-linear elasticity theory originated in this way, with special attention initially given to suitable ranges of finite but not too-large deformations, i.e., to ‘Ericksen-Pitteri neighborhoods’ (EPNs) in strain space [2, 4–6, 12], whereon the global lattice invariance reduces to point-group symmetry.¹ A large body of literature originated from such EPN-based approach, especially aiming at modelling reversible martensitic transformations [4, 5, 11, 14–17], also with the goal of improving the mechanical properties of shape-memory alloys, for instance to enhance their reversibility performance through the control of twinned-microstructure formation [18, 19].

Another line of research used the above Ericksen-Landau framework to model a wider class of phenomena in crystal mechanics, including reconstructive structural transformations where strains may attain or go beyond the EPN boundaries, producing defect nucleation and evolution in the lattice, and, in general, also to model phenomena directly related to the plastic behavior of crystalline materials, where the large deformations are not confined to any EPNs in strain space. Although discussions of the behavior of 3D crystals based on global lattice symmetry can be found in [7, 20–23], more systematic research has been done been done on crystal elasto-plasticity only in the 2D case. A family of Ericksen-Landau energies 2D crystals, obtained by patching suitable polynomials to obtain $\mathcal{C}^{2}$-smoothness and global invariance, was proposed in [6]. This allowed an improved understanding of the behavior of crystalline materials also in regimes of large deformations possibly outside the EPNs [8–10, 24–26].

A parallel line of work on the 2D case was based on the observations in [3, 27], where the natural tools proposed for a theory encompassing full lattice symmetry in 2D are modular functions, the well-known class of complex maps arising in diverse branches of Mathematics and Physics [28–32]. This led to the formulation of potentials suitable for 2D crystal plasticity, by following, in particular, the suggestion in [27] to construct Ericksen-Landau energies by means of the Klein modular invariant $J$ [28, 29, 33, 34]. This is akin to using a modular order parameter for crystal mechanics, extending earlier related notions such as the transcendental order parameter in [35, 36]. In this spirit, [25, 37] proposed an explicit class of smooth $J$-based strain energies with a unique ground state, up to full lattice symmetry, exploring the ensuing variational modelling of 2D crystal elasto-plasticity. We notice that smoothness, as obtained here and in [25, 37], is the standard assumption for potentials in Landau theory and elasticity theory, while the reduced $\mathcal{C}^{2}$-smoothness of [6] would induce jump discontinuities in the third- and higher-order moduli [38], affecting the connected material properties.

Here we continue these investigations by examining an explicit, simplest class of Ericksen-Landau $J$-based strain potentials for reconstructive transformations in 2D lattices, focussing on the most relevant case of energy functions exhibiting ground states with the two maximal symmetries, square and hexagonal ($s$-$h$), of 2D Bravais lattices. We explore some basic properties of the $s$-$h$ strain-energy landscapes moulded by global symmetry, in particular their bifurcation and valley-floor network, which are important in the selection of the activated deformation paths under total-energy minimization. We thus use a global $s$-$h$ potential in the simulation of a quasi-static shearing test, obtaining typical effects associated with the micro-mechanics of phase transformations in crystals. In particular, we observe strain avalanching underpinned by bursty coordinated basin-hopping activity of the local strain values under the slowly changing boundary conditions. This produces the inhomogeneous progress of the structural phase change, characterized by jagged stress relaxation via bursty microstructure development in the body, also mediated, in the reconstructive case, by defect nucleation and movement in the lattice. The present simulations also confirm the role, highlighted yet in [37], of the energy’s valley floors as largely establishing the deformation pathways for a crystal’s intermittent evolution under an external driving.

2 Strain Energies for 2D Crystalline Materials

2.1 The Strain Energy of 2D Crystals and Ericksen’s Proposal for Its Invariance

We consider a two-dimensional (2D) hyperelastic material, whose deformations are one-to-one maps $\mathbf{x}= \mathbf{x}(\mathbf{X})$, where the Cartesian coordinates $(x_{1},x_{2})$ identifying the current positions of material points $\mathbf{X}=(X_{1},X_{2})$ in a given reference state are considered with respect to a given ortho-normal basis $\{\mathbf{u}_{1},\,\mathbf{u}_{2}\}$. The deformation gradient $\mathbf{F}=\boldsymbol{\nabla }\mathbf{x}$ has matrix elements $F_{ij}=\partial x_{i}/\partial X_{j}$, and $\mathbf{C}=\mathbf{F}^{T}\mathbf{F}= \mathbf{C}^{T}>0$ is the symmetric, positive-definite the Cauchy-Green strain tensor. The strain-energy density $\sigma$ is a smooth real function of $\mathbf{C}$, and satisfying the material-symmetry requirements (1)₁-(1)₃:

$$ \sigma = \sigma (\mathbf{C})=\sigma (\mathbf{G}^{T}\mathbf{C} \mathbf{G}), \quad \mathbf{G}\in {\mathcal{G}}, \quad {\mathcal{G}} = E^{-1}GL(2,Z)E, $$

(1)

to hold for any $\mathbf{C}$ and for any tensor $\mathbf{G}$ in a suitable group $\mathcal{G}$ characterizing the response of the material. For crystalline substances we assume with Ericksen that the invariance group $\mathcal{G}$ be dictated by the material’s underlying lattice structure [1, 2], see also [3, 5, 11, 39]. In the 2D case under consideration here, this means that $\mathcal{G}$ should be a suitable conjugate to the group describing the global symmetry of 2D Bravais lattices, as is made explicit in Eq. (1)₄ above, with $E = (e^{h}_{j})$, for $\mathbf{e}_{j} = e^{h}_{j} \mathbf{u}_{h}$ (summation understood, with $j,h = 1, 2$), where $\{\mathbf{e}_{1},\,\mathbf{e}_{2}\}$ are the lattice basis in the reference state, and $\text{GL}(2,\mathbb{Z})$ denotes the group of unimodular (thus invertible) 2 by 2 matrices with integral entries [5]. For brevity we refer to assumption (1)₄ as to the GL-invariance of the density $\sigma$ in (1), which we split into the sum of a convex volumetric part $\sigma _{\text{v}}$, penalizing the departure of $\det \mathbf{C}$ from 1, and a distortive term ${\sigma }_{\text{d}}$ depending on the unimodular tensor $\bar{\mathbf{C}}=(\det \mathbf{C})^{-1/2}\mathbf{C}$:

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

(2)

Due to the GL-periodicity (1)₄, ${\sigma }_{\text{d}}$ in (2) is non-convex and only needs to be defined on a GL-fundamental domain in the space of unimodular strains, such as $\mathcal{D}$ made explicit in (5) below. As mentioned in [25], depending on the material and the effects of interest, in the numerical implementations of (1) the elements of the computational grid may represent suitable multiples of the actual lattice cell, in relation for instance to the conditions and scale of applicability of the Cauchy-Born Rule [40, 41] in the determination of the GL-periodic strain potential.

2.2 Modular Forms and GL-Invariant Strain Energies on the Poincaré Half Plane

Due to their GL-invariance, smooth potentials as in (1)-(2) are closely related to the modular functions on the Poincaré upper complex half-plane ℍ [3, 27]. This is best seen by smoothly mapping the space of 2D unimodular (positive-definite, symmetric) strain tensors $\bar{\mathbf{C}}$ bijectively to ℍ:

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

(3)

$$ \mathbb{H} = \{ x+iy\in \mathbb{C}, y>0 \}, \qquad (ds)^{2}=\big[(dx)^{2}+(dy)^{2} \big]/y^{2} \;, $$

(4)

where (4)₂ is the 2D hyperbolic metric on ℍ [42–44], and $\bar{C}_{ij}$ are the components of $\bar{\mathbf{C}}$ in the basis $\{\mathbf{u}_{1},\, \mathbf{u}_{2}\}$. Through the complex parameterization (3) of the unimodular-strain space, the material-symmetry maps $\bar{\mathbf{C}}\mapsto \mathbf{G}^{T}\bar{\mathbf{C}}\mathbf{G}$, for $\mathbf{G}\in {\mathcal{G}}$, become maps $\hat{z}(\bar{\mathbf{C}}) \mapsto \hat{z}(\mathbf{G}^{T} \bar{\mathbf{C}}\mathbf{G})$ which are all the isometries of ℍ given by the group of linear fractional transformations with integral entries, supplemented by the map $z \mapsto -\bar{z}$ (see [37] for details). Figure 1 represents geometrically this action, where the fundamental domain²

$$ \mathcal{D}=\{ z\in \mathbb{H} : \mathopen{\lvert }z \mathclose{\rvert }\ge 1,\ 0\le \text{Re}(z)\le \tfrac{1}{2} \}, $$

(5)

together with the collection of all its non-overlapping hyperbolically-congruent GL-related copies $m^{t}\mathcal{D}m$, $m \in \text{GL}(2,\mathbb{Z})$, produce the Dedekind tessellation [31, 45] of ℍ. In light of this, [27] suggested that, as a main building block to obtain GL-invariant smooth potentials ${\sigma }_{\text{d}}$ in (1), one can use the Klein invariant $J$. This complex map, given explicitly in the Appendix, is a $\text{SL}(2,\mathbb{Z})$-periodic holomorphic function on ℍ [28, 29, 33, 34], the modular group $\text{SL}(2,\mathbb{Z})$ being the positive-determinant subgroup of $\text{GL}(2,\mathbb{Z})$. Indeed, owing to the properties of $J$, it is possible to consider $J$-based Ericksen-Landau strain-energy functions ${\sigma }_{\text{d}}$ as in (1)-(2) by setting

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

(6)

where the smooth function ${\sigma }_{\text{d}}(J)$ should be such that it guarantees [37]: (a) the full GL-periodicity (1) for (6), rather than the sole invariance under $\text{SL}(2,\mathbb{Z})$ exhibited by $J$; and, (b) ensure the existence of a positive-definite elastic tensor for any stable lattice configuration, thus avoiding the existence of equilibria with degenerate Hessian, as required by elasticity theory. Examples of such potentials, suitable for the elasto-plasticity of 2D crystals, and for their reconstructive transformations, are discussed explicitly hereafter.

Fig. 1

(a) Dedekind tessellation of the Poincaré half-plane ℍ [30–32]. The positions are indicated of nine GL-equivalent square points ($i$, $i + 1$, $\zeta = \frac{1}{2}(i+1)$, $\zeta +1$, …, in blue), and four GL-equivalent hexagonal points ($\rho = e^{i \pi /3}$, $\rho -1$, …, in purple). The GL-equivalent points appear to be getting closer to each other the closer they get to the real axis, but are equidistant in the hyperbolic metric (4)₂ on ℍ. This rectangle and the indicated red dashed arc of the unit circle about the origin of ℍ refer also to Fig. 2. (b) Plot of the GL-periodic strain energy with square minimizers, see the function $\sigma _{i}$ in (8) with $\mu =1$. The plot is on the rectangle in panel (a), so that we observe nine energy wells, with bottoms at the nine blue square GL-equivalent points of panel (a)

2.3 Strain Energies for 2D Crystal Plasticity

In [25, 37] were analyzed some simplest forms of GL-invariant $J$-based strain-energy functions ${\sigma }_{\text{d}}$ as in (6), exhibiting a single ground-state configuration in each GL-copy of the fundamental domain $\mathcal{D}$. In this case, let the equilibrium lattice be given by the strain $\bar{\mathbf{C}}_{0}$, and set $z_{0} = \hat{z}(\bar{\mathbf{C}}_{0}) \in \mathcal{D}$, so that this is the only minimizer of ${\sigma }_{\text{d}}$ in $\mathcal{D}$. Then for any $z_{0}$ such that $J'(z_{0})\neq 0$ (that is, for any ground state except for square and hexagonal ones: $z_{0} \neq i, \rho$) the simplest $J$-based GL-invariant strain-energy function in (6) has the form

$$ \sigma _{z_{0}}\big(\bar{\mathbf{C}}\big) = \mu |J(z) - J(z_{0})|^{2}, $$

(7)

where $\mu >0$ is an elastic modulus, and $z = \hat{z}(\bar{\mathbf{C}})$ as in (3). In the case of the maximally symmetric ground states $z_{0}=i$ (square) or $z_{0}=\rho$ (hexagonal), taking into account that $J(i) = 1$, $J(\rho )= 0$, $J'(i) = 0$, $J'(\rho ) =J''(\rho ) = 0$, we have that the simplest $J$-based GL-potentials with non-degenerate elastic moduli at their minimizers are respectively given by [25, 37]:

$${\sigma }_i\left(\overline{\mathbf{C}}\right)=\mu |J(z)-1|\text{for }{z}_0=i\text{(square)},$$

(8)

$${\sigma }_{\rho }\left(\overline{\mathbf{C}}\right)=\mu {|J(z)|}^{2/3}\text{for }{z}_0=\rho \text{(hexagonal)}.$$

(9)

As an example, we show in Fig. 1 a portion of the GL-periodic energy landscape on ℍ given by the square energy (8).

3 Ericksen-Landau Theory for Reconstructive Transformations in Crystalline Materials

3.1 Simplest $J$-Based Strain Energies for the Square-Hexagonal Transformation

The functions in (7)-(9) can be used to construct GL-potentials of the type (6) with more than one minimizer in the fundamental domain $\mathcal{D}$, so that they are suitable for crystals which may undergo structural phase changes between different stable equilibrium configurations of their lattice (see Footnote (3)). Here we consider explicitly the case of reconstructive transformations, and their most relevant case, in which the two energy minimizers in $\mathcal{D}$ are the maximally symmetric points $i$ (square) and $\rho$ (hexagonal). This will produce a GL-invariant potential suitable for the $s$-$h$ transformation in 2D Bravais lattices.³

A simplest class of such $s$-$h$ densities is given by the following normalized linear combination of the two functions in (8)-(9):

$$ \begin{aligned} \sigma _{\text{R}}\big(\bar{\mathbf{C}}\big)&=\sigma _{i} \big(\bar{\mathbf{C}}\big)+\beta \sigma _{\rho}\big(\bar{\mathbf{C}} \big) \\ &=\mu |J(z)-1|+\beta \mu |J(z)|^{2/3} , \end{aligned} $$

(10)

where $z=\hat{z}(\bar{\mathbf{C}})$ as in (3) and (6), and where we consider $\beta >-\frac{3}{2}$. The modulus $\mu$ is here a scale factor which will henceforth be set to 1, so that (10) defines a one-parameter family of potentials with critical points $i$ and $\rho$, whose relative height is controlled by $\beta$ as shown in Fig. 2.

Fig. 2

GL-invariant strain energy for the square-hexagonal phase transformation. (a) Plot of the $s$-$h$ potential $\sigma _{\text{R}}$ in Eq. (10) for $\beta = 1$ and $\mu =1$, on the domain in Fig. 1(a), whereon the energy $\sigma _{\text{R}}$ exhibits thirteen energy wells: nine well bottoms are located at the nine equivalent square points, and four at the equivalent hexagonal points, shown in Fig. 1(a). (b) Section of the plot of the $s$-$h$ energy $\sigma _{\text{R}}$ in (10) with $\mu =1$ and varying $\beta$, taken along the unit-circle shown in red in Fig. 1(a). The value $\phi =\frac{\pi}{2}$ corresponds to the square point $i$, while $\phi =\frac{\pi}{3}$ and $\phi =\frac{2\pi}{3}$ correspond to the two neighbouring hexagonal points $\rho -1$ and $\rho$, with rhombic points given by generic values of $\phi$, see also Fig. 1 and Footnote (2)

3.2 Bifurcation and Valley Floors

We show in Fig. 3(a) the bifurcation on the plane ($\beta$, $y$), with $x=\frac{1}{2}$, for the critical points of the $s$-$h$ energy $\sigma _{\text{R}}$ in (10). This is obtained from the way the global GL-symmetry of the potential constrains, via the implied local (point-group) symmetry, the second-derivatives of its critical and bifurcation points [2, 5]. This diagram expectedly has the same main features as the one pertaining to the polynomial-based $s$-$h$ energy in [6]. The actual GL-periodicity of the bifurcation pattern of the energy in (10) is sketched in Video V1 of the Supplementary Material (SM).

Fig. 3

Bifurcation and valley floors for the $s$-$h$ energies. (a) Section on the plane ($\beta$, $y$), for $x=\frac{1}{2}$, of the GL-invariant bifurcation diagram for the critical points of the $s$-$h$ energy $\sigma _{\text{R}}$ in Eq. (10) for $\beta > -\frac{3}{2}$. Dotted and solid lines indicate unstable and stable critical-point branches, respectively, with the square-$i$ (solid red line) and hexagonal-$\rho$ (solid green lines) which coexist as local minimizers for $0<\beta <\frac{3}{2}$, with rhombic saddles in between (dotted blue lines). See also Video V1 in the SM. (b) GL-invariant hyperbolic network (a Bethe-like tree) of the valley floors on ℍ for the strain energy $\sigma _{\text{R}}$ in Eq. (10), for $\beta =1$. Nodes (blue and purple symbols) are at the $s$-$h$ minimizers, and edges are along those geodesics on ℍ which contain a pair of $s$-$h$ minimizers (for clarity only the arcs joining such $s$-$h$ points are marked, by dotted red lines). The rhombic saddles mentioned in panel (a) are marked by black dots. See also Figs. 4-5, and Videos V2, V3 in the SM

In Fig. 3(a) we see that the square critical point $i$ of $\sigma _{\text{R}}$ is stable for $\beta <\frac{3}{2}$, losing stability at $\beta = \frac{3}{2}$ through a subcritical pitchfork to two rhombic critical points (saddles). On the other hand, $\rho$ is a minimum for $\beta >0$, becoming unstable at $\beta =0$, where symmetry dictates the presence of a monkey saddle, unfolding [47] via a transverse bifurcation to three rhombic critical-point branches (standard saddles, only one of which belongs to the plane ($\beta$, $y$) of the figure; see also Video V1 in the SM). The points $i$ and $\rho$ are the only local minimizers of $\sigma _{\text{R}}$ in $\mathcal{D}$ for $0<\beta <\frac{3}{2}$, and in this range the energy $\sigma _{\text{R}}$ in (10) is thus suitable to model the $s$-$h$ transformation. The coexisting minima $i$ and $\rho$ have the same energy at the Maxwell value $\beta = \beta _{M} = 1$, so that the global minimum is $i$ for $0\leq \beta \leq 1$, while it is $\rho$ for $1\leq \beta \leq \frac{3}{2}$.

For the $s$-$h$ energy $\sigma _{\text{R}}$ in (10) with $0<\beta <\frac{3}{2}$, it is also interesting to highlight the structure of the infinite GL-invariant network of the valley floors connecting the energy extremals on ℍ. As mentioned in the Introduction, this considerably inform us regarding the evolution of the strain field in ℍ. Precisely, the valley floors of an energy $\sigma$ are the gradient-extremal loci on ℍ locally satisfying $\mathbf{H}\,\boldsymbol{\nabla }\sigma - \xi \, \boldsymbol{\nabla } \sigma =0$, with $\xi \in \mathbb{R}$ and $\mathbf{H}\mathbf{v}\cdot \mathbf{v}\geq 0$ [48–50], where $\boldsymbol{\nabla }\sigma$ and $\mathbf{H}$ are respectively the gradient and Hessian of $\sigma$, and $\mathbf{v}$ is any vector orthogonal to the energy gradient, with derivatives and orthogonality considered in relation to the hyperbolic metric (4)₂. This produces the following explicit gradient-extremal equation:

$$ \frac{2}{y^{3}}\,|\boldsymbol{\nabla }\sigma |^{2}\, \mathbf{v}_{2} - \frac{2}{y^{2}} \mathbf{H}\,\boldsymbol{\nabla }\sigma + \xi \, \boldsymbol{\nabla }\sigma =0, $$

(11)

where $\mathbf{v}_{2}$ is the unit vector in the $y$-direction on ℍ, and derivatives are intended in its standard atlas $(x, y)$. For $\beta =1$ the valley floors of the $s$-$h$ energy $\sigma _{\text{R}}$ obtained from (11) compose the hyperbolic network highlighted in Fig. 3(b), whose edges lie on those geodesics of ℍ [42, 43] which contain both the $s$-$h$ minimizers.

4 Shear-Driven $s$-$h$ Transformation

We investigate numerically the behavior of an $s$-$h$ phase-transforming crystal in quasi-static simple shear, imposed to the top side of a square body containing a square lattice aligned with the body sides, with fixed bottom and the remaining two body sides free. In this incremental test, for each value of the shear parameter $\gamma$ a local minimizer of the body’s total strain-energy functional is computed through the density $\sigma _{\text{R}}$ in (10) with $\beta =1$, and complying with the imposed boundary conditions (see the Appendix and [37] for details).

Figures 4(d)-(e)-(f) show three snapshots of the resulting $\gamma$-dependent strain field in the sheared crystal. Figures 4(a)-(b)-(c) highlight the corresponding strain clustering as a cloud of points which evolves with $\gamma$ on the GL-domains of the Dedekind tessellation of ℍ. See Video V2 in the SM for the numerical simulation of shearing up to $\gamma = 0.24$.

Fig. 4

Shearing of an $s$-$h$ phase-transforming crystal, with strain energy $\sigma _{\text{R}}$ in (10) for $\beta =1$. The imposed loading is along a primary-shear direction in the square lattice, parallel to the driven horizontal body-sides, see the Appendix. The associated path in ℍ is the straight dashed blue line $i \rightarrow 1+i$ in panels (a), (b), (c), with increasing shear parameter $\gamma$ (green dot), from the defect-free reference configuration in the ground state $z_{0}=i$ ($\gamma =0$) to the GL-equivalent fully-sheared square configuration $i+1$ ($\gamma =1$). Convexity domains around each $s$-$h$ energy minimizer are shaded gray, and the valley floors of $\sigma _{\text{R}}$ are in dashed-red as in Fig. 3(b). The snapshots (a), (b), (c) show the evolution of the strain clustering during shear, given by the heatmap 2D-histogram for the cell-strain density evolving on the Dedekind tessellation of ℍ in Fig. 1(a). Panels (d), (e), (f) show the associated body deformation (with same color coding for the cell strains in ℍ as in panel (a)), with panel (g) displaying the stress-strain relation (blue jagged line, with black dots marking the snapshots) for increasing $\gamma$. The response is elastic to about $\gamma = 0.13$, after which a bursty phase-transformation regime begins. Panels (d), (e), (f) show this is marked by the formation of $s$-$h$ phase mixtures (twin bands and lath-type microstructures) mediated by evolving lattice defects, as seen in the detail inset to panel (e). The deformation’s intermittency is tracked in panel (g), through the sequence of relaxation events given by the observed jumps in the jagged stress-strain diagram, as well as the orange spikes indicating the percentage of cell-strain values that are changing energy basin at each increment of $\gamma$. Panels (a), (b), (c) show that the strain-cloud path on ℍ under the imposed boundary condition here follows the directions of the valley floors in the energy landscape, as is the case in crystal plasticity [37]. More information on the bursty deformation triggered in this shear test is also in Figs. 5-6 and Videos V2, V3, V4 in the SM

The stress-strain relation in Fig. 4(g), where $F_{x}$ is the integrated shear stress on the top boundary, shows that the initially defect-free lattice begins shearing with a significant elastic charge, the associated strain cloud widening away from $i$ in ℍ, as $\gamma$ moves away from 0, due to the growing strain heterogeneity caused by the unloaded body-sides, see Figs. 4(a)-(d). A large transformation event at about $\gamma = 0.13$ ends the elastic regime, with a large stress drop taking place as part of the strain cloud in ℍ splits away from its initial location near the reference state $i$ towards the neighbouring well in $\rho$. A portion of the cells’ strains remain far from the well bottoms, elastically stabilized on the intermediate non-convex regions, see Figs. 4(b)-(c).

From there on, the imposed shearing induces a bursty deformation process in the body, characterized by an intermittent sequence of stress-relaxation events due to avalanching $s$-$h$ microstructure formation assisted by lattice-defect evolution, as can be seen in Figs. 4(e)-(f)-(g), and Video V2 in the SM. These phenomena occur as the strain field in the lattice locally takes advantage, for the relative minimization of the total energy, of the available $s$-$h$ GL-wells of the density $\sigma _{\text{R}}$. Twin-type bands and dislocations emerge as neighboring lattice cells suitably stretch or shear and rotate while satisfying Hadamard’s kinematic compatibility, together with the imposed boundary conditions.⁴ Long-range elastic interactions correspondingly produce coordinated basin-hopping on ℍ which results in strain avalanches within the shearing body under the slow driving, see also [25, 26, 37]. These complex deformation mechanisms involving both phase transition and defect evolution, occur here with no need for auxiliary hypotheses: they originate directly from energy minimization not only due to the GL-arrangement of the density minimizers in strain space, but also to the way in which the global GL-symmetry shapes the whole energy landscape. Indeed, we observe here in Figs. 4(d)-(e)-(f), as is the case also with crystal plasticity [37], that the GL-energy valley floors act as strain-cloud deformation pathways on ℍ. For instance, the creation and evolution of the vertical shear bands in Figs. 4(e)-(f), is due to strain avalanches occurring as suitable lattice domains leave the square reference state $i$, with the strain cloud following the valley-floor path $i \rightarrow \rho \rightarrow \zeta$ in ℍ, see Figs. 4(b)-(c), although the body is being externally loaded in the horizontal principal shear direction $i \rightarrow i+1$ (green dot in the same figures) for the square crystal. The activation of the deformation pathway $i \rightarrow \rho \rightarrow \zeta$ implies $s$-$h$ phase transformation events happening together with, and assisted by, dislocational effects in the lattice, as some cells’ strains respectively reach the $\rho$-well (hexagonal) or the $\zeta$-well (fully sheared square through a principal lattice-invariant shear) when the driving forces them away from $i$. We see here that in the present variational GL-modelling plastification may arise in the lattice via defect nucleation through lattice-invariant shears [6–10, 24, 25], because in the reconstructive case the barriers to the these shears are only as high as the barriers relative to the phase transformation itself. The stress distribution of the dislocations matches well the linear elastic predictions for the far field, see Fig. 1 of [23] and Fig. 4 of [25].

Figure 5 and Video V3 in the SM display explicitly the shearing body’s strain cloud as it flows in quasi-static intermittent fashion along the energy-surface valley floors. Most of the cell’s strains are located near the involved well bottoms, with a fraction elastically stabilized on the non-convex regions between wells. Our simulations confirm the role of the GL-network of valley floors as giving the deformation pathways for the strain cloud on ℍ during total-energy minimization, and show, as in [37], that these features of the energy landscape can help to better inform also other crystallographically-based approaches to crystal micromechanics, such as the phase-field models in [23, 52], and may usefully guide the analysis of their relations with fully variational formulations such as in the present setting.

Fig. 5

Snapshot of the bursty evolution of the strain field on ℍ ($\gamma = 0.24$) shown directly on a portion of the energy landscape during the shear-driven $s$-$h$ transformation. The GL-energy is $\sigma _{\text{R}}$ in Eq. (10) with $\beta = 1$, so the $s$-$h$ wells $i$, $\rho$, $\zeta$, $i+1$, …, all have the same depth. The avalanching strain values largely follow the energy valley floors indicated in dashed red, as in Fig. 3(b). See also Video V3 in the SM

As a concluding remark, we highlight explicitly the strain avalanches characterizing the intermittent deformation in this shear test, showing one such event in Fig. 6. It originates from a burst of coordinated basin-hopping activity of the strain values on ℍ, associated to microstructure and defect evolution in the lattice, producing a stress-relaxation drop. The computed sequence of such strain avalanches observed in our quasi-static simulation is shown in Video V4 of the SM, and the corresponding serrated stress-strain relations, as in Fig. 6, are observed experimentally in mechanically-induced reconstructive transformations in crystals [53]. Our modelling predicts this to occur through strain intermittence as in Video V4, in close analogy to the strain-intermittence experimentally observed during mechanically-induced martensitic transformation in shape-memory alloys [54]. Such behavior is expected at least at the early stages of the mechanically induced reconstructive transformation. For the latter, the developing plasticity effects, also evidenced by the model, would eventually impede microstructure formation and mobility, leading to the transition’s irreversibility (Footnote 4, [7]), unlike with the case of memory martensites.

Fig. 6

Highlight of a strain burst characterizing the shear-driven $s$-$h$ transformation in Fig. 4. This strain-evolution event in the body undergoing quasi-static shear ($\gamma = 0.20$) corresponds to a spike (in orange) in the intermittent basin-hopping activity of the local strains on ℍ, with the associated relaxation drop in the stress-strain relation (blue jagged line). The strain avalanche is computed by considering the norm of the strain difference at each cell for two consecutive values of $\gamma$ in the simulation, near $\gamma = 0.20$. See Video V4 in the SM for more details

Appendix

1.1 A.1 Variational Problem and Computational Details

In the numerical simulation the body is square-shaped containing an initially a defect-free square lattice ($z_{0} = i$), with the lower, fixed side, aligned with a dense square-lattice direction, in turn parallel to the vector $\mathbf{u}_{1}$ in the given ortho-normal basis. The simple shear boundary condition is $\mathbf{I} +\gamma \, \mathbf{u}_{1}\otimes \mathbf{u}_{2}$, imposed to the top side of the square body, with shear parameter $\gamma$ such that $\gamma =0$ at the ground state $z_{0} = i$ and $\gamma =1$ at the fully-sheared square ground state $i+1$. The corresponding shearing path in ℍ is the straight dashed blue line $i \rightarrow 1+i$ in panels (a), (b), (c) of Fig. 4. The remaining two sides of the square body are free.

For each $\gamma$ a local minimizer of the body’s total strain-energy functional is computed through the density (10), with $\beta =\beta _{\text{M}} =1$, and complying with the imposed boundary conditions given by the increasing $\gamma$. In the minimization algorithm the volumetric energy $\sigma _{\text{v}}(\det \mathbf{C})$ in (2) has the form

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

(12)

with modulus $\lambda$ such that $\lambda /\mu = 30$ to obtain quasi-incompressibility. Furthermore, the function $j(z) \equiv 1728J(z)$ used to compute $\sigma _{\text{R}}$ in (10), has, as a basic property, all integers in its Fourier coefficients [28, 29, 33, 34], and is given explicitly by

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

(13)

where $q = \exp (2\pi i\tau )$.

The variational problem derived from minimizing the body’s total energy functional is discretized through a classical FEM implemented through the free-source Finite-Element software FreeFEM++ [55], using piece-wise affine functions on a mesh for a total of about $3\times 10^{4}$ degrees of freedom. As both the energy density $\sigma _{R}$ and its related functional are highly non-convex, the variational problem shows an infinite number of metastable states in which minimization algorithms such as gradient-descent or Newton’s methods perform poorly. A suitable preconditioner was thus used via spectral equivalence. At each $\gamma$-step in the shear test the initial guess coincides with the metastable state associated with the previous value of the imposed loading. The algorithm then searches for metastable states which are close (in the sense of the preconditioned gradient-descent) to the previous equilibrium. The FEM discretization introduces in this non-convex problem a regularizing length scale related to lattice discreteness.

1.2 A.2 Caption to Supplementary Video V1

Animation showing the GL-periodic $\beta$-bifurcation diagram for the critical points of the $s$-$h$ potential $\sigma _{\text{R}}$ in (10).

A section of this diagram is shown in Fig. 3(a). The Poincaré disk model [6, 39, 46, 56] is used here for the 2D hyperbolic space. The purple square and green circles in the disk respectively correspond to square and hexagonal points, with fat [slim] rhombic points given by bold [thin] blue lines, while red curves represent rectangular points. The 3D-diagram corresponds to increasing $\beta$ from bottom to top. Vertical lines indicate stable (bold) or unstable (dotted) square and hexagonal critical points. Stability ranges for these $s$-$h$ extremals are as detailed in Fig. 3(a) and in the text. In particular, there is the $s$-$h$ coexistence interval $\frac{1}{2}<\beta <\frac{3}{2}$ where both the square and hexagonal points are local minimizers, with branches of rhombic saddle points bifurcating from these two maximally symmetric ones with features and multiplicities as described in Fig. 3(a): three transverse rhombic branches issue from each hexagonal bifurcation point, while two rhombic branches issue with a subcritical pitchfork from each square bifurcation point. For completeness here are also indicated the supercritical pitchforks from square to rectangular local minimizers of the potential $\sigma _{\text{R}}$ in (10) for $\beta \leq - \frac{3}{2}$.

1.3 A.3 Caption to Supplementary Video V2

Shearing of a homogeneous square body containing an initially defect-free square lattice for an $s$-$h$ phase-transforming crystal with strain energy $\sigma _{\text{R}}$ in (10) with $\beta =1$. The imposed loading is along a primary shear direction in the square lattice, aligned with both the parallel square-cell side and bottom body side. The shearing boundary condition is imposed through the constrained horizontal sides of the body, with the two remaining sides free. See also Figs. 2-3-4 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 ℍ visited by the strain of each lattice cell during loading. Defect nucleation and evolution accompany the phase transformation, intrinsically produced in this model through energy minimization. See also Fig. 4(d)-(e)-(f).

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

(b) Evolution of the 2D histogram (strain cloud) of the density of strain parameters on the Dedekind tessellation of ℍ during shearing. The histogram color-coding provides the percentage of body cells with strain located at each point of ℍ. The straight horizontal dashed-blue line between the two neighboring square configurations $i$ and $i+1$ is the image in ℍ 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$. Shading indicates the convexity domains around the $s$-$h$ minimizers in ℍ of the energy $\sigma _{\text{R}}$, while the red dashed lines mark the valley floor segments as in Fig. 3, which largely direct the strain-cloud evolution under the slow driving. See also Fig. 4(a)-(b)-(c) and Fig. 5.

(c) Color-code map used in (a) for the GL-energy basins of $\sigma _{\text{R}}$ on the Dedekind tessellation of ℍ. The ridges marking the basins’ boundaries are computed via a Eq. (11) in the main text.

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

(e) Bursty $s$-$h$ phase transformation in the shearing test. Jagged stress-strain body behavior (blue), with the underlying spikes (orange) showing the percentage of basin-hopping strain values during loading, as $\gamma$ grows. The body response is elastic to about $\gamma = 0.13$, where, after a first large stress drop, intermittent stress-relaxation continues for growing $\gamma$. See also Fig. 4(g) and Fig. 6.

(f) Evolution of the histogram for the values of the angle $\theta$ in the polar decomposition of the deformation gradient $\mathbf{F}$, indicating local lattice rotation accompanying the phase transformation process. Notice the different rotation angles of the $s$-$h$ phase-microstructure bands, resulting in an evolving bimodal distribution for $\theta$.

1.4 A.4 Caption to Supplementary Video V3

Quasi-static bursty evolution of the strain field on ℍ for growing parameter $\gamma$, shown directly on the energy landscape during the shear-driven $s$-$h$ transformation. The GL-energy is $\sigma _{\text{R}}$ in Eq. (10) with $\beta = 1$, so the $s$-$h$ wells $i$, $\rho$, $\zeta$, $i+1$, …, all have the same depth. Under the slow driving the strain values largely follow the energy valley floors indicated in dashed red in Fig. 3(b) and Fig. 5. See also panel (b) in Video V2 in the SM.

1.5 A.5 Caption to Supplementary Video V4

Strain avalanching during the shear-driven $s$-$h$ phase transformation.

(c)-(a) For reference, these two panels report respectively the evolution of the strain field in the shearing body, as in Video V2(a), and the associated jagged stress-strain relation, from Video V2(e).

(b) Intermittent strain avalanching within the shearing crystal. Each event corresponds to phase transformation and defect evolution in the deforming lattice. Avalanches are computed by considering the difference in strain norm at each lattice cell for two consecutive values of $\gamma$ in the simulation. The GL-energy is $\sigma _{\text{R}}$ in Eq. (10) with $\beta = 1$, as in Fig. 2. See also Fig. 6.