Strain tuning of highly frustrated magnets: Order and disorder in the distorted kagome Heisenberg antiferromagnet

Strain applied to a condensed-matter system can be used to engineer its excitation spectrum via artificial gauge fields, or it may tune the system through transitions between different phases. Here we demonstrate that strain tuning of the ground state of otherwise highly degenerate frustrated systems can induce novel phases, both ordered and disordered. For the classical Heisenberg antiferromagnet on the kagome lattice, we show that weak triaxial strain reduces the degeneracies of the system, leading to a classical spin liquid with non-coplanar configurations, while stronger strain drives the system into a highly unconventional state which displays signatures of both spin-glass behavior and magnetic long-range order. We provide experimentally testable predictions for the magnetic structure factor, characterize the ground-state degeneracies and the excitation spectrum, and analyze the influence of sample shape and boundaries. Our work opens the way to strain engineering of highly frustrated magnets.

Strain applied to a condensed-matter system can be used to engineer its excitation spectrum via artificial gauge fields, or it may tune the system through transitions between different phases. Here we demonstrate that strain tuning of the ground state of otherwise highly degenerate frustrated systems can induce novel phases, both ordered and disordered. For the classical Heisenberg antiferromagnet on the kagome lattice, we show that weak triaxial strain reduces the degeneracies of the system, leading to a classical spin liquid with non-coplanar configurations, while stronger strain drives the system into a highly unconventional state which displays signatures of both spin-glass behavior and magnetic long-range order. We provide experimentally testable predictions for the magnetic structure factor, characterize the ground-state degeneracies and the excitation spectrum, and analyze the influence of sample shape and boundaries. Our work opens the way to strain engineering of highly frustrated magnets.

I. INTRODUCTION
The manipulation of many-body systems by external stimuli is widely used in both the search for novel phenomena and the realization of applications. Among the possible tools, pressure and the resulting lattice deformations are particularly appealing, as they do not introduce disorder -as opposed to chemical substitution -and they can either preserve or modify in a controlled fashion the lattice symmetries of the underlying system. Recent experimental progress in applying uniaxial or otherwise inhomogeneous forces has led to the notion of straintronics, where specific strain patterns enable to engineer states and functionalities of novel materials [1][2][3].
Examples of strain manipulation include the mechanical switching of nano-electronic graphene devices [4], the occurrence of strain-induced Landau levels in graphene, with a spacing corresponding to ultra-large magnetic fields [5][6][7], the creation of artificial gauge fields for ultracold atoms and photonic crystals [8,9], the proposals to realize Landau levels for emergent charge-neutral excitations in solids, such as magnons or Majorana spinons in quantum antiferromagnets [10,11], and the modification of multi-component superconducting states [12,13]. In essentially all of these examples, one starts from a unique microscopic state which is modified by strain, changing either its static properties or its excitation spectrum.
In this paper, we extend the concept of strain tuning to highly degenerate many-body systems: Here, applying strain can be expected to have a singular effect, i.e., even small strain modifies the system's properties in a qualitative fashion. We choose to discuss the effect of inhomogeneous strain applied to highly frustrated antiferromagnets where degenerate ground-state manifolds result from competing interactions. We show that suitably chosen strain patterns can be used to induce particular forms of magnetic order as well as novel spin-liquid regimes, thus opening a new arena for strain-based engineering of states of matter. Specifically we consider the classical Heisenberg antiferromagnet on the kagome lattice. Its highly degenerate ground state features Coulombic spin correlations, and it displays remarkably complex order-by-disorder phenomena at finite low temperature [14][15][16][17][18][19]. Its spin-1/2 cousin is a prime candidate to realize a quantum spin-liquid ground state [20][21][22][23][24][25][26][27][28][29]. Here we focus on the effect of triaxial strain which partially preserves discrete lattice symmetries of the kagome lattice, Fig. 1, while partially relieving strong geometric frustration. Increasing strain lifts the classical degeneracies of the unstrained system, first deforming the spin liquid into a non-coplanar one with pronounced spin correlations at Q = 0. Larger strain induces a transition into a highly unusual state, being connected with the inability to independently minimize the energy on every triangle. This state displays characteristics of a spin glass, but at the same time its magnetic structure factor shows sharp peaks corresponding to Q = 0 long-range order. While details of the ground-state configurations and lowenergy excitation modes depend on the sample shape and boundaries, the gross features of the strained magnetic state appear robust. We connect our findings to known results for homogeneous uniaxial strain applied to the kagome Heisenberg model, and we comment on the role of quantum effects.
The remainder of the paper is organized as follows: In Sec. II we introduce the inhomogeneously strained Heisenberg model and discuss the re-writing of its Hamiltonian as sum of complete squares. Sec. III describes the results for spin configurations and the spin structure factor, obtained from minimizing the classical energy. Sec. IV then discusses the complex energy landscape at finite strain, implying glassy features, and the properties of the low-energy excitations. A summary and discussion of open questions closes the paper. (a-d): Distorted kagome-lattice antiferromagnet, with displacements from triaxial strain, Eq. (3), for β = 1 and |C| = 0.025/a0; the undistorted lattice is shown in light gray. Longer (shorter) bonds correspond to weaker (stronger) exchange couplings Jij as indicated by the color code. Arrows indicate the force applied to the sample. Panels (a,b) show the system with c-type boundaries, panels (c,d) with ic-type boundaries. Left (a,c) and right (b,d) panels correspond to positive and negative strain, respectively. The linear system size is N = 6. (e) Schematic phase diagram of the distorted kagome-lattice antiferromagnet in the limit T → 0, for details see text.
As an aside, we note that the effect of strain on a kagome-lattice tight-binding model has been recently studied in Ref. 30, with focus on single-particle pseudomagnetic fields. Also, spontaneous (instead of imposed) distortions which relieve frustration in highly frustrated magnets have been discussed in earlier papers [31][32][33].

A. Kagome Heisenberg model
We consider a nearest-neighbor antiferromagnetic Heisenberg model on the kagome lattice, formed by corner-sharing triangles, with spatially varying couplings: We exclusively focus on the classical case and treat the S i as unit vectors. The homogeneous system, J ij ≡ J, features a ground-state manifold which includes both coplanar and non-coplanar states [14]. Finite-temperature fluctuations tend to select coplanar states via an orderby-disorder mechanism [14][15][16][17], and the low-temperature regime displays weak long-range spin order corresponding to a √ 3 × √ 3 ordering pattern [18,19].

B. Triaxial strain
Strain engineering goes back to the discussion of electron-phonon coupling in carbon nanotubes [34] where strain-induced modulations of hopping matrix elements can emulate a vector potential for electrons [5]. Central to our work is the modification of magnetic exchange couplings due to strain. In the distorted lattice, each magnetic ion is characterized by a displacement vector U i . This results in exchange couplings between neighboring ions, entering the Hamiltonian (1), which we assume to follow where a 0 is the reference bond length and For instance, the hopping matrix elements t of graphene display a bond-length dependence with β t ≈ 2 . . . 3 [3]; for exchange couplings following J = t 2 /U where U is an on-site Coulomb repulsion this would mean β ≈ 4 . . . 6. Note that Eq. (2) represents a linear approximation to the full (typically exponential) bond-length dependence of the exchange constant; assuming an exponential dependence yields qualitatively similar results as shown in the Appendix. Most numerical results are shown for β = 1.
In the following, we focus on triaxial strain where the displacement vector is given by [5,36] U (x, y) =C 2xy, withC encoding the distortion amplitude, and we employ The dimensionless parameter C =Cβa 0 specifies of the modulation strength of the J ij . The distortions described by Eq. (3) increase linearly with increasing distance from the sample center. Structural stability then requires to consider finite-sized samples. Combining Eqs. (2) and (3) we define -for fixed sample size -a maximum strain C max beyond which the outermost couplings become formally negative due to the linearization in Eq. (2). For β = 1 this means that the longest (i.e. weakest) bond takes twice its original length at maximum strain. This maximum strain is inversely proportional to the linear system size, C max ∝ 1/N , therefore the thermodynamic limit N → ∞ cannot be taken at fixed strain. As our results show, it is instead meaningful to consider the thermodynamic limit at fixed C/C max , i.e., N → ∞ with (CN ) fixed. Moreover, we will sometimes refer to the combined limit β → ∞ andC → 0, keeping C =Cβa 0 fixed, which reduces non-linearities in the strain dependence of couplings [35]. In fact, the results at fixed C depend only weakly on β for β > 20.
In order to partially preserve the discrete lattice symmetries, we primarily consider samples of triangular shape [35], Fig. 1. Here, discrete rotation and mirror symmetries exist w.r.t. the sample center. For such samples, positive and negative strain,C ≷ 0, correspond to qualitatively distinct distortion patterns, and we will display results for both.
Since finite-size properties will depend on the structure of the edges, we consider different types of edges: Kagome-lattice frustration is best preserved for edges with complete kagome triangles (dubbed c-type), Fig. 1(a,b). As a representative for different edges, we choose those with all outer spins removed such that the outward triangles are incomplete (dubbed ic-type), Fig. 1(c,d). In both cases, we denote the linear sample size by N where N counts the number of complete triangles along a sample edge. Then, the total number of spins is N s = (3/2)N (N + 1). Depending on (N mod 3) the center of the sample is either formed by an elementary triangle or hexagon.
For C > 0 there are three weakest bonds located in the sample corners. For c-type edges we find C + max = √ 3/(4N − 3) and for ic-type edges C + max = √ 3/(4N − 9), both valid for arbitrary β. For C < 0 there are now six weakest bonds in the corner triangles. C max is given by a lengthy expression which is not particularly enlightening. However, in the limit β → ∞ it simplifies to |C − max | = √ 3/(2N − 3) for both c-type and ic-type edges.

C. Constraint satisfiability and critical strain
The homogeneous Heisenberg Hamiltonian on the kagome lattice can be written as sum of complete squares, and this rewriting can be generalized to inhomogeneous couplings [37]: For each triangle α with spins ijk we can define γ iα = (J ij J ik /J jk ) 1/2 such that the Hamiltonian With small triaxial strain applied, the γ iα will weakly deviate from their unstrained reference value unity, such that the minimization constraint L α = 0 can be fulfilled for all triangles of a finite sample. In contrast, for larger strain the γ iα do no longer fulfill the triangle inequality for triangles α near the sample corners or edges, depending on the sign of C. This change defines a critical value of strain, C ± crit , where ± correspond to positive and negative strain, respectively.
While these considerations strictly apply to samples with c-type edges, samples with ic-type edges contain bonds not belonging to triangles, rendering the system less frustrated. Hence, the nature of the ground-state Critical value of triaxial strain, plotted as C ± crit /C ± max , as function of inverse linear system size, 1/N , for different β for (a,b) positive strain and (c,d) negative strain, both for samples with (a,c) c-type boundaries and (b,d) ictype boundaries. Ccrit/Cmax approaches the value 0.4 in the limit β → ∞, N → ∞ in all cases, for details see text.
manifold depends on the type of sample edges. However, this difference turns out to be of minor importance for the magnetic structure factor for sufficiently large samples.
Numerical results for the ratio of critical and maximum strain, C crit /C max , for triangular samples are shown in Fig. 2. This quantity displays a mild dependence on system size N , but a stronger dependence on the magnetoelastic coupling β. For |C| < |C ± crit | the strained system is strongly frustrated, with L α = 0 ∀α defining a degenerate manifold of liquid-like states [37,38], while for |C| > |C ± crit | the constraint L α = 0 cannot be fulfilled for all triangles. Then, the condition α L 2 α → min induces tendencies to magnetic order, as we will see in the next section.
The spatial profile of constraint satisfiability can be visualized by plotting the quantity | L α | for each triangle α in the ground state, this is in Fig. 3. With increasing positive strain, the constraint is first violated in the sample corners. For samples with c-type edges, we are able to find an analytic expression in the limit of β → ∞ which reads C + crit = √ 3/(10N − 9). Hence, in this limit we have which again tends to 0.4 for N → ∞. For C < 0, the constraint can be satisfied in the corner triangles for any strain up to C max . However, the constraint gets first violated for the triangles in the middle of the boundaries. While we have not been able to obtain a closed-form expression, our numerical evaluation shows that C − crit /C − max → 0.4 for β → ∞ and N → ∞ for both c-type and ic-type edges, as in the case of positive strain, Fig. 2.
For |C| > |C crit | the strained system contains multiple triangles where the constraint L α = 0 cannot be satisfied. If we define the number of these triangles as M us , we can consider its ratio with the total number of triangles M tot (which is N s /3 for c-type edges). Apparently r ≡ M us /M tot is zero (non-zero) for |C| < |C crit | (|C| > |C crit |), respectively. In the limit of large system size, N → ∞, r approaches a finite value which remains smaller than unity for |C| = |C max | as triangles near the center of the sample remain weakly distorted even in this limit. Numerical results for r are shown in Fig. 4.
Parenthetically, we note that for |C| < |C crit | all constraints L α = 0 can be satisfied, but the ground states cannot be mapped to an origami analog as discussed in Ref. 38, because the geometric condition of Eq. (2) in that paper is in general not fulfilled by the couplings of the strained system.

III. NUMERICAL RESULTS: CONFIGURATIONS AND SPIN STRUCTURE FACTOR
We now turn to our core numerical results, obtained for finite-size triaxially strained kagome-lattice Heisenberg systems using system sizes up to N = 24.

A. Iteration scheme
We use an iterative scheme to find spin configurations corresponding to local minima of the total energy in configuration space. For given values of N , β, and C which determine the Hamiltonian we start from a random initial spin configuration and iteratively minimize the total energy by aligning each spin according to its mean field, supplemented by appropriate random mixing. The iteration is aborted once the average energy per bond, ε, changes less than a threshold ε conv in one step. For ε conv = 10 −8 J this happens after typically 10 4 . . . 10 5 iteration steps. The iteration is repeated for N init = 10 5 different initial conditions. For the state with the globally lowest energy, E min , we denote by ε min its energy per bond, ε min = E min /N b , where N b is the number of bonds. Given the SU(2) spin symmetry of the underlying Hamiltonian, two of the resulting spin configurations are considered equivalent if they match (within a numerical threshold) up to global SU(2) rotations.
Convergence tends to be slow for large systems due to the glassy nature of the energy landscape, see Sec. IV below. Computation time therefore limits our ability to reach larger system sizes, and most calculations are restricted to N ≤ 20.

B. Ground-state spin configurations
For any finite strain, C = 0, we find that non-coplanar spin configuration are energetically preferred over coplanar ones: We have verified this tendency by comparing the ground-state energies between those for the SU(2)symmetric model and models with varying degree of easyplane anisotropy, obtained by reducing the prefactor of the S z S z coupling. The easy-plane models yield a consistently higher ground-state energy, except at zero strain where the ground-state energy does not depend on the anisotropy.
Representative ground-state configurations near maximum strain are shown in Fig. 5. While the spin configurations are naturally inhomogeneous, a clear tendency towards local three-sublattice 120 • order is visible near the sample center; for the unstrained kagome lattice such order is known as Q = 0 order [16]. In contrast, near the sample corners (for C > 0) or edges (for C < 0), triangles with ferrimagnetic-like configurations (↑↑↓) prevail. This can be rationalized by noting that, in these strongly distorted regions, the spatial distribution of coupling constants, Fig. 1, corresponds to locally uniaxial strain.
In fact, uniaxially strained kagome antiferromagnets have been investigated before [39][40][41][42] and display regimes of ferrimagnetism. Denoting the couplings along one direction by J and along the two others by J , the situation J J corresponds to chains weakly coupled via middle spins, while the case J J realizes a square lattice with spin-decorated bonds. In the classical limit, the ground state is a collinear ferrimagnet for J/J < 1/2, while for J/J > 1/2 there is an infinite family of degenerate canted ferrimagnetic ground states. While the former state is stable also for quantum spins S = 1/2, the latter is most likely replaced by a spiral state for large J/J [41], and a spin liquid is present in the intermediate regime. The J/J < 1/2 ferrimagnet is of obvious relevance to our triaxially strained system, and we will get back to this below.

C. Spin structure factor
We have analyzed the strain-induced states quantitatively by determining their static spin structure factor, defined as where . . . denotes an average over the ground-state manifold. For a state with magnetic LRO at wavevec-tor Q, the value of m 2 = S( Q)/N s corresponds to the squared order parameter m in the thermodynamic limit. In practice, we use N init different initial conditions to find local minima in the energy landscape as described in Sec. III A. For each converged configuration, we determine the average energy per bond, ε = E/N b . The averaging in Eq. (4) is then performed over those of the N init converged states whose energy per bond ε falls in the window [ε min , ε min + ∆ε]. The finite energy window ∆ε accounts for inaccuracies in convergence and may be interpreted in terms of a finite temperature; it is chosen sufficiently small, ∆ε = 10 −6 J, unless noted otherwise.
For c-type edges this procedure ensures averaging over a representative set of configurations from the continuously degenerate ground-state manifold. For ic-type edges the sampling is over the ground states plus a small number of low-energy excited states to improve statistics due to the glassy energy landscape which is characterized by many local minima close to the ground-state energy, see Sec. IV below. In calculating S( q) we have varied ∆ε for selected parameter sets and found that choosing smaller ∆ε changes the values of S( q) by less than 5% for the system sizes used.
In general, the calculated spin structure factor displays qualitatively similar behavior for positive and negative strain, for different values of β, and for different sample edges, as illustrated in Figs. 6, 7, and 8. At small strain, Figs. 6(a,e) and 7(a,e), the structure factor has the broad shape familiar from the classical kagome Heisenberg model [18,43], with pinch points located at reciprocal wavevectors Q = Γ , the centers of higher Brillouin zones, characteristic of the U(1) spin liquid. These pinch points remain sharp under strain, but gain weight with increasing strain, Fig. 6(b,f). At larger strain, Fig. 6(c,d,h), pronounced peaks at Γ emerge, which grow with increasing |C| and correspond to three-sublattice Q = 0 order. Differences between positive and negative strain appear minor. Individual differences, e.g. between Figs. 7(b) and (f), can be attributed to the different values of C/C crit which can be read off from Fig. 2.

D. Q = 0 order
Finite-size scaling for the height of the peaks in S( q) at q = Γ and fixed C/C max is demonstrated in Figs. 9(a,b) for β = 1. The data clearly show that S( Q)/N s scales to zero as N → ∞ for small |C|, but tends to a finite value at larger |C|. This signals the existence of a magnetically ordered state at large |C|, with the transition being located at C = C crit within our accuracy, Fig. 9(c,d). Consistent with this, we find that the width of the peaks in S( q) scales to zero as N → ∞ for |C| > |C crit |. Fig. 9 also illustrates non-monotonic N dependencies which can be traced back to commensurability effects of inhomogeneous spin arrangements near the sample boundaries. A finite-size-scaling comparison between samples with c-type and ic-type edges for positive strain  Fig. 6, but now for β = 100. The results are qualitatively similar to that shown in Fig. 6 for β = 1, with the quantitative differences reflecting the β dependence of Ccrit.
is shown in Fig. 10, here for β = 100. While finite-size systems with ic-type edges tend to have a larger order parameter and a smaller non-monotonic N dependence, the extrapolated data show only minor differences, partially related to the different C/C crit . Together, this underlines that the Q = 0 order is a robust bulk phenomenon, with boundary effects being subleading.
To rationalize the appearance of the Q = 0 order, we recall that the triaxially strained kagome flake realizes bond configurations at the corners (at the mid-edges) for C > 0 (C < 0) which correspondi to a uniaxially strained systen with J/J < 1/2. The latter displays ferrimagnetic order [39][40][41][42] with magnetic Bragg peaks at two of the Γ momenta, namely those perpendicular to the J bond direction. Therefore, the six Γ peaks present for |C| > |C crit | in the triaxially strained sample can be thought of as a superposition of three domains of ferrimagnetic configurations. We note, however, that this picture is oversimplified, as (i) the entire sample contributes to the Bragg peaks and (ii) the local spin configurations deviate significantly from the collinear ferrimagnet.
A particular situation arises for samples with ic-type edges and C < 0: Here, the sample edges tend to have configurations of the form ↑↑↓↓, Fig. 5(d), because -compared to c-type edges -the outer triangles with strong bonds have been removed, producing pairs of effectively FIG. 8. Spin structure factor S( q) as in Fig. 6, but now for samples with ic-type edges and β = 100. For positive strain (top) the results are qualitatively similar to that for c-type edges, Fig. 7, apart from a slightly smaller value of Ccrit. In contrast, for negative strain incommensurate correlations dominate for |C| > |C − crit |; those become commensurate only for larger systems, as shown in Fig. 11 below. strongly coupled spins along the edge. The ↑↑↓↓ configurations in turn lead to structure-factor peaks at incommensurate locations, see Fig. 8(f-h). As these configurations are restricted to the boundary row, their influence on the structure factor diminishes with increasing N , such that for sufficiently large samples commensurate peaks in S( q) are restored, Fig. 11. For strongly distorted lattices, it makes a difference whether the lattice coordinates used to calculate S( q) in Eq. (4) are taken as the unstrained R i or the strained R i + U i ; a neutron-scattering experiment would probe the latter. For simplicity, Figs. 6-11 have been calculated with unstrained coordinates. A comparison of the structure factors calculated with both unstrained and strained coordinates is shown in Fig. 12 for a large-strain case and a realistic value of β = 3: One sees that the straininduced change in lattice geometry leads to a broadening of the Bragg peaks. For larger β the geometric lattice distortion at a given C/C max is smaller and hence the difference between the two cases diminishes. For completeness we have also considered kagome flakes of shapes different from triangular, in particular hexagonal and circular. For sufficiently large sam- ples, the structure-factor results are similar to those for triangular-shaped samples. We illustrate this in Fig. 13 which shows spin configuration and structure factor for a hexagonal-shaped system under triaxial strain. We note that, for such samples, positive and negative strain are equivalent by symmetry.

IV. GLASSY ENERGY LANDSCAPE
Our iterative minimization scheme keeps track of a large set of local minima in the energy landscape. For the problem at hand, we hand found the latter to be surprisingly complex.

A. Local minima and their energy distribution
For very small systems, N 5, we find that essentially all iterations converge to states with the same lowest energy. Inspecting the spin configurations for |C| > |C crit | indicates the existence of few discrete degenerate ground states for ic-type boundaries, whereas systems with ctype boundaries appear to display continuously degenerate ground states (except at C = C max ), as repeated iterations find inequivalent states with identical energies.
The situation is drastically different for larger systems. While convergence to the same lowest energy is still common for |C| < |C crit |, the iteration scheme finds local minima with widely distributed energies for |C| > |C crit |. Sample distributions for the energy per bond, ε − ε min , are shown in Fig. 14: The distributions appear effectively continuous for large systems, with a width reaching up to 10 −4 J. Naturally, the distribution width increases with increasing strain C, while for C → 0 the width tends to zero.
These results signal a complex energy landscape for |C| > |C crit |, with abundant local minima which are separated by energy barriers. This behavior is well known for spin glasses where randomness and frustration conspire to produce complex low-energy states without long-range order [44]. Remarkably, we find a glassy energy landscape here in a system free of quenched disorder (but with inhomogeneous couplings), and our structure-factor results show that this glassy behavior coexists with signatures of magnetic long-range order.

B. Ground-state degeneracy
We made an attempt to estimate the ground-state degeneracy -up to global SU(2) rotations -by monitoring the number of different [45] converged spin configurations, N diff , whose energy equals the minimum energy within a small window ∆ε = 10 −9 J. For samples with c-type boundaries we find that N diff scales with N init for all values of C and system sizes N , indicating that the ground states are continuously degenerate. This is consistent with the finding of non-trivial zero modes reported below. Histograms of converged energies per bond, ε, representing local minima of the energy landscape. Data are shown for different values of C/Cmax and different system sizes N for samples with c-type boundaries and β = 1. The horizontal axis shows the energy relative to the global minimum, ∆ε = ε − εmin, on a logarithmic axis, the vertical axis the frequency out of 10 4 random initial conditions.
In contrast, for samples with ic-type boundaries we observe that N diff tends to saturate with increasing N init , at least for strain values |C crit | < |C| < |C max |. The actual number of ground states depends on both C and N , with non-monotonic variations, see Table I. We point out, however, that counting true ground states is a numerically expensive task for large systems due to the glassy nature of the energy landscape, and a compromise between runtime, convergence accuracy, and selection window ∆ε is required. Based on the available data, we are not able to determine how the ground-state degeneracy of ic-type samples scales with system size; the results clearly point toward a non-extensive number.

C. Low-energy modes
To further characterize the states of the strained kagome Heisenberg magnet, we determine the quadratic energy cost of fluctuations around ground-state configurations. To this end, we construct the Hessian matrix for a given minimum-energy state and determine its eigen-values and eigenvectors.
The Hessian is constructed as described in Ref. 46: For a spin configuration {s i } we choose an orthonormal local basis (s i , u i , v i ) at every lattice site and parameterize fluctuations ass i = 1 − 2 i s i + ui u i + vi v i with i = ( ui , vi ) which takes into account the spin normalization condition. The quadratic energy cost of fluctuations around a ground state is given by E = T M where M is the Hessian matrix with dimension 2N s × 2N s . Diagonalizing the Hessian matrix gives us the eigenvalues λ j and the corresponding eigenvectors.
The spectrum of the Hessian always contains three trivial zero modes (Goldstone modes) due to the SU(2) symmetry of the underlying Hamiltonian. Interpreting all eigenvalues below 10 −6 J as zero modes, we find that samples with c-type boundaries generically display additional, i.e., non-trivial zero modes implying a continuous ground-state degeneracy, except for C = C + max . The zeromode count is shown Table II. The number of zero modes increases with linear system size N for |C| < |C crit |, while it appears to saturate for |C| > |C crit |. We have analyzed the zero-mode eigenvectors by calculating their inverse participation ratio and by inspecting their spatial distribution (not shown) and concluded that the non-trivial modes primarily live near the c-type sample boundaries for any non-zero strain. Boundary-induced zero modes are in fact consistent with the analysis of bond-disordered kagome antiferromagnets in Ref. 37 which concluded that no zero modes should exist in the bulk for inhomogeneous distributions of magnetic couplings.
In contrast, samples with ic-type boundaries do not feature non-trivial zero modes, except at C = C − max where the corner spins are disconnected and can be trivially rotated. This implies that the ground states display discrete degeneracies only (apart from global rotations), again consistent with Ref. 37.
The character of the finite-energy spectrum can be analyzed via the cumulative distribution function F (λ) = 1/(2N s ) j Θ(λ − λ j ) of the Hessian eigenvalues λ j for a particular local minimum. Plots of the cumulative distribution function are shown in Figs. 15 and 16. Independent of the edges, the spectrum appears gapless for large N : existing gaps in the spectrum get filled with increasing N , indicating that these are finite-size effects. Interestingly, for c-type (ic-type) edges the density of low-energy modes increases (decreases) with increasing |C|. This can be rationalized by considering that, with increasing |C|, for c-type edges zero modes are converted  6  18  18  6  3  18  18  6  9  8  24  24  12  3  24  24  9  9  10  30  30  10  3  30  30 10 9 into low-E modes, whereas for ic-type edges all modes are shifted to higher energy.

V. SUMMARY AND OUTLOOK
Non-uniform strain can be used to drive highly frustrated magnets into novel states: We have demonstrated this for triaxial strain applied to the classical kagome Heisenberg antiferromagnet: While this model system, in the absence of strain, is in a highly degenerate classical spin-liquid state, weak strain partially lifts the degeneracies. The system enters a non-coplanar spin-liquid state with pronounced strain-driven short-range spin correlations. Larger strain drives a phase transition into a state with Q = 0 long-range order, and we have connected this to the tendency towards ferrimagnetism in the uniaxially strained kagome antiferromagnet. Most interestingly, the inhomogeneously ordered state at large strain displays a rugged energy landscape akin to that of a spin glass.
Our results demonstrate an intriguing coexistence of magnetic long-range order and a glassy energy landscape in a classical non-random spin system. This calls for a deeper understanding of its dynamic properties, not only at the linear-response level, but also concerning relaxation and quenches. In this context, the role of the zero modes present for c-type boundaries needs particular attention. At finite temperatures, strain effects will compete with thermal order by disorder, which may drive novel types of phase transitions.
A notoriously difficult question is that for quantum effects at T = 0. Perhaps most interesting is the physics of the strained S = 1/2 kagome Heisenberg antiferromagnet. Provided that the unstrained system is a topological Z 2 spin liquid, it features a gap to Z 2 vortex excitations (visons) and hence can be expected to be stable at least against small strain. However, numerics indicates that the unstrained system is sensitive to small perturbations, related to its proximity to one or more transitions between different ground-state phases [25][26][27]29]. Hence, moderate strain is likely sufficient to modify the quantum ground state: Given the rugged energy landscape, we believe the strained kagome quantum antiferromagnet presents a fascinating platform to study quantum glassiness and aspects of many-body localization.
Our work suggests to consider strain engineering of degenerate states on a more general level and, in particular, prompts generalizations to other strain patterns as well as highly frustrated magnets on other lattices, such as pyrochlore or hyperkagome. While this is left for future work, we speculate here that large strain generically induces ordering tendencies, and it will be extremely interesting to study the emergence of corresponding ordering transitions, both at zero and finite temperatures.
FIG. 17. Spin structure factor S( q) as in Fig. 6, but here for strained samples with c-type edges, β = 3, and exponential length dependence (A1) of coupling constants J. The results are similar to that obtained using the linearized length dependence, Eq. (2). The results presented in the body of the paper employ the simplified geometry dependence of the magnetic exchange couplings given in Eq. (2). In a real material this dependence is not linear and will in general also depend on bond angles. The angle dependence arises from anisotropic orbitals as well as from superexchange paths via intermediate ions and is clearly material-dependent.
Here we illustrate the robustness of our findings by as-suming an exponential bond-length dependence instead of the linearized one, i.e., we use J ij = J exp −β(| δ ij |/a 0 − 1) . (A1) As these couplings are always positive, our previous definition of C max ceases to be well-defined. Therefore, we now use CN =CN βa 0 as size-independent dimensionless measure of the strain effect on the coupling constants.
We have performed numerical simulations using Eq. (A1) instead of Eq. (2) and find the results to be qualitatively unchanged: Upon increasing the strain, both the non-coplanar spin liquid and the glassy Q = 0 ordered state appear in an essentially unchanged fashion. This is illustrated in Figs. 17 and 18, showing the spin structure factor and the finite-size scaling for its peak height; these figures can be compared to Figs. 6, 7, 9, and 10(a,c).
The robustness can be rationalized as follows: Compared to its linearized version, the exponential coupling dependence (A1) leads to somewhat larger couplings for long (i.e. weak bonds) and to significantly larger couplings for short (i.e. strong bonds), while undistorted bonds remain unchanged. As a result, elementary triangles whose linearized couplings (2) strongly violate the triangle inequality γ iα + γ jα > γ kα , see Sec. II C, and thus cannot satisfy the constraint L α = 0, continue to do so for an exponential dependence (A1). As unsatisfied triangles force the emergence of the ordered glassy phase, its appearance and character remains unchanged. Switching from linearized to exponential coupling dependence shifts the transition location C crit by about 15% (recall, e.g., C + max N = √ 3/4 ≈ 0.433 for N → ∞). We conclude that the linearization of the couplings' length dependence, Eq. (2), is a permissable approximation in the regime of interest. We note that the same linearization approximation is frequently used in the literature on strained graphene where it has been shown to remain reasonably accurate for sample deformations up to 10-15% [36,47].