Mechanical instability at finite temperature

Many physical systems including lattices near structural phase transitions, glasses, jammed solids, and bio-polymer gels have coordination numbers that place them at the edge of mechanical instability. Their properties are determined by an interplay between soft mechanical modes and thermal fluctuations. In this paper we investigate a simple square-lattice model with a $\phi^4$ potential between next-nearest-neighbor sites whose quadratic coefficient $\kappa$ can be tuned from positive negative. We show that its zero-temperature ground state for $\kappa<0$ is highly degenerate, and we use analytical techniques and simulation to explore its finite temperature properties. We show that a unique rhombic ground state is entropically favored at nonzero temperature at $\kappa<0$ and that the existence of a subextensive number of"floppy"modes whose frequencies vanish at $\kappa = 0$ leads to singular contributions to the free energy that render the square-to-rhombic transition first order and lead to power-law behavior of the shear modulus as a function of temperature. We expect our study to provide a general framework for the study of finite-temperature mechanical and phase behavior of other systems with a large number of floppy modes.


I. INTRODUCTION
Crystalline solids can undergo structural phase transitions in which there is a spontaneous change in the shape or internal geometry of their unit cells [1][2][3][4]. These transitions are signalled by the softening of certain elastic moduli or of phonon modes at a discrete set of points in the Brillouin zone. Lattices with coordination number z = z c = 2d in d spatial dimensions, which we will call Maxwell lattices [5] [6], exist at the edge of mechanical instability, and they are critical to the understanding of systems as diverse as engineering structures [7,8], diluted lattices near the rigidity threshold [9][10][11], jammed systems [12][13][14], biopolymer networks [15][16][17][18][19], and network glasses [20,21]. Hypercubic lattices in d dimensions and the kagome lattice and its generalization to higher dimensions with nearest-neighbor (NN) Hookean springs of spring constant k are a special type of Maxwell lattice whose phonon spectra have harmonic-level zero modes not at a discrete set of points but at all N (d−1)/d points on (d − 1)-dimensional hyperplanes oriented along symmetry directions and passing through the origin [22]. A question that arises naturally is whether these lattices can be viewed as critical lattices at the boundary between phases of different symmetry and, if so, what is the nature of the two phases and the phase transition between them.
Here we introduce and study, both analytically and with Monte-Carlo simulations, a square-lattice model (easily generalized to higher dimensions) in which nextnearest-neighbors (NNNs) are connected via an anharmonic potential consisting of a harmonic term with a spring constant κ tuned from positive to negative and a quartic stabilizing term. When κ > 0, the square lattice is stable even at zero temperature. When κ = 0, NNN springs contribute only at anharmonic order, and the harmonic phonon spectrum is identical to that of the NN-lattice. When κ < 0, the NNN potential has two minima, and the ground state of an individual plaquette is a rhombus that can have any orientation (Fig. 1c). Plaquettes in the same row (or column), however, are constrained to have the same orientation, but plaquettes in adjacent rows can either tilt in the same direction or in the mirror-image direction as shown in Figs. 1d-f, leading to 2 × 2 √ N equivalent ground states and a subextensive but divergent entropy of order √ N ln 2. The properties of this model, including the subextensive entropy at zero temperature, are very similar to those of colloidal particles confined to a low-height cells [23,24] and to the anti-ferromagnetic Ising model on a deformable triangular lattice [25]. In addition, the scaling of the shear modulus near the zero-temperature critical point is analogous to that observed in finite-temperature simulations of randomly diluted lattices near the rigidity-percolation threshold [26] and to finite-temperature scaling near the jamming transition [27], suggesting that generalizations of our model and approach may provide useful insight into the thermal properties of other systems near the Maxwell rigidity limit.

II. RESULTS
Strong fluctuations arising from the large number of zero modes lead to interesting physics at T > 0 in this square lattice model. We show the following: • Among all the equal-energy zigzagging configurations at κ < 0, the uniformly sheared rhombic lattice, shown in Fig. 1d, has the lowest free energy: the ground-state degeneracy is broken by shows the uniformly sheared rhombic lattice, which we show to be the preferred configuration at small T in the thermodynamic limit. (e) is a randomly zigzagging configuration, and (f) is the ordered maximally zigzagging configuration, which has a unit cell consisting of two particles.
• Thermal fluctuations lead to a negative coefficient of thermal expansion and corrections to the shear rigidity [35][36][37][38] that enable the square lattice state to remain thermodynamically stable in a region of the phase diagram at κ ≤ 0 (Fig. 2).
• Fluctuations drive the transition from the rhombic to the square phase first order and lead to the phase diagram shown in Fig. 2 in which the temperature of the transition approaches zero as κ → 0 − .
• The low-T shear modulus G (Fig. 2) in both phases is proportional to |κ| at low temperature, and there is a critical regime with T > const. × κ 3/2 in which G ∼ T 3/2 . In addition, there is a region in the square phase (mostly metastable with respect to the rhombic phase) with T < const. × |κ| 3/2 in which G ∼ (T /|κ|) 2 . This behavior near the T = 0 critical point is analogous to that found in the randomly diluted triangular lattice [26] near the central-force rigidity threshold. Interestingly, the critical regime in our model is fundamentally a consequence of nonlinearity as is the case for dynamical scaling near the jamming transition [27]. These predictions are supported by our Monte-Carlo simulations and by direct calculations of entropic contributions to the free energy phonon fluctuations in different arrangements of kinks.

III. DISCUSSION
Our model of the square-to-rhomic transition is very similar to a model, studied by Brazovskii [39], for the transition from an isotropic fluid to a crystal and later applied to the Rayleigh-Bénard instability [40] and to the nematic-to-smectic-C transitions in liquid crystals [41,42]. In all of these systems, there is a subextensive but infinite manifold of zero modes [a (d − 1 = 2)dimensional spherical shell in the Brazovskii case, a (d − 1 = 1)-dimensional circle in the Rayleigh-Bénard case, and two (d − 2 = 1)-dimensional circles in the liquid crystal case] in the disordered phase leading to a singular contribution to the free energy. We use the Brazovskii theory to calculate the temperature of the firstorder square-to-rhomic transition as a function of κ < 0 and negative thermal expansion in the square phase.
In applying the Brazovskii approach to our problem, we develop an expansion of the free energy that maintains rotational invariance of elastic distortions in the target space. Previous treatments [1] of structural transitions tend to mix up nonlinear terms arising from nonlinearities in the strain tensor required to ensure rotational invariance and nonlinear terms in the elastic potential itself. This advance should be useful for the calculation of renormalized free energies and critical exponents in standard structural phase transitions.
The number and nature of zero modes of the critical NN square lattice has a direct impact on the properties of the lattices with κ > 0 and κ < 0, and it is instructive to review their origin. A powerful index theorem [43] relates the number of zero modes of a frame consisting of N points and N B ≡ 1 2 zN bonds, where z is the average coordination number, via the relation where S is the number of independent states of self-stress in which bonds are under tension or compression and in which the net force on each point is zero. In his seminal 1864 paper [5], Maxwell considered the case with S = 0 that yields the Maxwell relation for the critical coordination number at which N 0 is equal to the number n(d) (= d for periodic and d(d + 1)/2 for free boundary conditions) of zero modes of rigid translation and rotation: (3.2) In the limit of large N , z N c → z ∞ c = 2d. There are many small unit cell periodic Maxwell lattices with N 0 = S. The NN square and kagome lattices in two dimensions and, the cubic and pyrochlore lattices in three dimensions are special examples of these lattices that have sample-spanning straight lines of bonds that support states of self-stress under periodic boundary conditions. They, therefore, have of order N (d−1)/d states of self-stress and the same number of zero modes, which are indicators of buckling instabilities of the lines when subjected to compression. Geometrical distortions of theses lattices that remove straight lines, as is the case with the twisted kagome lattice [44], remove states of selfstress and associated zero modes. When subjected to free rather than periodic boundary conditions, these distorted lattices have a deficiency of order N (d−1)/d bonds and as a result the same number of surface zero modes (there are no bulk zero mode other than the trivial ones of uniform translation), which can have a topological character [45] or be described in the long-wavelength limit by a conformal field theory [44]. Unlike the infinitesimal zero modes of hypercubic lattices, those of the kagome and pyrochlore do not translate into finite zero modes of the lattices when finite sections are cut from a lattice under periodic boundary conditions. Thus, it is not yet clear whether the ground state of the latter lattices will are highly degenerate or not. Nevertheless, Brasovskii theory should provide a sound description of thermal properties of these lattices in the vicinity of the point T = 0, κ = 0.

IV. MODEL AND ORDER-BY-DISORDER IN THE LOW-SYMMETRY PHASE
The model we consider is a square lattice with two different types of springs -those connecting nearest neighbors and those connecting next-nearest neighbors, as shown in Figure 1(a). The NN springs are Hookian, with potential where k > 0. The NNN springs are introduced with an anharmonic potential where κ can be either positive or negative and g, introduced for stability, is always positive. The Hamiltonian of the whole lattice is thus, where R i is the positions of the node i and a is the lattice constant. In what follows, we will use the reduced variables to measure the strength of couplings in V NNN . For κ > 0, V NNN (x) has a unique minimum at x = 0, and the ground state of H is the square lattice with lattice spacing a. All elastic moduli of this state are nonzero, and it is stable with respect to thermal fluctuations, though, as we shall see, it does undergo thermal contraction at nonzero temperature. When κ < 0, V NNN (x) has two minima at x = ± 6κ/g, corresponding to stretch and compression, respectively. This change in length of NNN springs is resisted by the NN springs, and in minimum energy configurations one NNN bond in each plaquette will stretch and the other will contract. The alternative of having both stretch or contract would cost too much NN energy. A reasonable assumption, which is checked by our direct calculation, is that the stretching and contraction will occur symmetrically about the center so that the resulting equilibrium shape is a rhombus rather than a more general quadrilateral (see Supplementary Information Sec. I). The shape of a rhombus is uniquely specified by the lengths d 1 and d 2 of its diagonals (which are perpendicular to each other), whose equilibrium values are obtained by minimizing the sum over plaquettes of The grounds state of the entire lattice when κ < 0 must correspond to a tiling of the plane by identical rhombi each of whose vertices are four-fold coordinated. It is clear that zigzag arrangements of rows (or columns) of rhombi in which adjacent rows tilt in either the same or opposite directions constitute a set of ground states. A derivation showing that this is the complete set can be found in Ref. [25], which considered packing of isosceles triangles, which make up half of each rhombus. The ground state energy per site, ǫ 0 , is simply V PL evaluated at the equilibrium values of d 1 and d 2 .
Each ground-state configuration of a system with N x vertical columns and N y horizontal rows has K = 0, · · · , N y horizontal zigzags or K = 0, · · · , N x vertical zigzags. Thus, the ground state entropy diverges in the thermodynamic limit, though it is sub-extensive and proportional to N x + N y in a system of N = N x N y particles. Such ground state configurations are found in other systems, most notably the zigzagging phases seen in suspensions of confined colloidal particles [23,24]. The confined colloidal system has a phase diagram that depends only on the planar density and the height of confinement of the colloids. For sufficiently large heights, the colloids form a phase of two stacked square lattices. In a neighboring region of the phase diagram, explored in simulations of Refs. [46,47], this square lattice symmetry is broken through a weakly discontinuous transition and a rhombic phase is observed. This region of the phase diagram of confined colloids thus provides a physical realization of the Hamiltonian (4.3).
At low but nonzero temperatures, the degeneracy of the ground state is broken by thermal fluctuations through the order-by-disorder mechanism [25,28,[30][31][32][33][34]. This splitting of degeneracy due to small phonon fluctuations around the ground state may be calculated using the dynamical matrix in the harmonic approximation to the Hamiltonian (4.3). For each ground state configuration R i , we write the deformation as R i → R i + u i , and expand to quadratic order in u, H = 1 2 ij u i D ij u j . The Fourier transform of D ij , D q , is block-diagonal and the phonon free energy, which is purely entropic, for that configuration is (4.6) where w p is the free-energy per site in units of k B T . In general, F p depends not only on K, but also on the particular sequence of zigzags. We numerically calculated this free energy for all periodically zigzagged configurations with up to 10 sites per unit cell. We found that the lowest-free-energy state is the uniformly sheared state [ Fig. 1(c)] with K = 0, and the highest-free-energy state is the maximally zigzagged sate [ Fig. 1(e)] with K = N y and two sites per unit cell. Both of these energies are extensive in the number of sites N , and we define Figure 3(a) displays our calculation of ∆w 0 (τ ), which vanishes, as expected, at τ = 0 and also at large τ at which the rhombus collapses to a line. Figure 3(b) plots ∆F/∆F 0 as a function of φ = K/N y for different values of τ . By construction, this function must vanish at φ = 0 and be equal to one at φ = 1. All of the points lie approximately on a straight line of slope one. Thus, we can approximate F p (N x , N y , K) by Note that for each K, this energy is extensive in N = N x N y as long as φ = 0.
These calculations were carried out in the thermodynamics limit, N → 0, in which the sum over q is replaced by the continuum limit by an integral. In order to compare these results with the Monte Carlo results of the next section, it is necessary to study finite size effects. The first observation is that for any finite N x , the system will be effectively one-dimensional for a sufficiently large N y , and as a result, we would expect the number of zigzags to fluctuate. To proceed, we continue to use the continuum limit to evaluate Eq. (4.6), and we use Eq. where E 0 (N x , N y ) denotes the potential energy which is independent of K, and the full free energy is where e 0 and f 0 are the potential energy and phonon free energy per site of the uniformly sheared state. Thus, as expected, when N x ∆w 0 ≫ 1, zigzag configurations make only a very small, subextensive contribution to the free energy. On the other hand, in the opposite limit, they make an extensive contribution of N x N y k B T ∆w 0 to the energy. Therefore, at a given τ , the zigzag configurations are favored when the system is small, and in thermodynamic limit, the rhombic configuration is always favored. Our Monte Carlo simulation verified this [see inset of Fig. 3(b)].

V. SIMULATION
We simulate the system using a Monte Carlo (MC) algorithm inside a periodic box whose shape and size The phonon contribution to the free energy for various zigzagging configurations. (a) shows the free energy difference ∆w0 between the maximally zigzagging configuration (Fig. 1e) and the uniformly sheared square lattice configuration (Fig. 1c). For sufficiently large τ , the lattice collapses onto a line, at which point the free energy difference goes to zero. (b) shows the free energy of the phonons as a ratio ∆F/∆F0 for τ = −0.02, −0.1 and −1 (λ = 10), where ∆F and ∆F0 are defined in Eq. (4.7) and ∆F is evaluated for all possible configurations with unit cells with at most 10 sites as a function of the zigzag fraction φ = K/Nx. The lattice without zigzags is entropically favored for any any value of τ in our calculations. ∆F/∆F0 is well approximated by the line φ, which corresponds to non-interacting zigzags. This interaction results in the dispersion of values of ∆F/∆F0. The inset of (b) shows Monte Carlo simulation results for the zigzagging fraction φ plotted against the linear system size √ N .
are allowed to change in order to maintain zero pressure. In this version of the Metropolis algorithm, also used in Refs. [25,48], for each MC step a particle is picked at random and a random trial displacement is performed. The trial displacement is initially uniformly distributed within a radius of 0.1a, but throughout the simulation the radius is adjusted to keep the acceptance probability between 0.35 and 0.45. Given the initial configuration energy E i and the trial configuration energy E j , the trial configuration is accepted with probability [1 + exp(E i − E j )/T ] −1 , i.e., Glauber dynamics is used. After initializing the system using the square lattice configuration with lattice constant a, the simulation is first run at a high temperature and is then annealed to the final low temperature. For each intermediate temperature, an equilibration cycle in a sample of N sites consists of at least N 4 × 10 2 MC steps. To accommodate areal and shear distortions in the different phases we encounter, the simulation box area and shape are changed using a similar acceptance algorithm, with the trial deformation adjusted to keep the acceptance probability between 0.35 and 0.45, such that the simulation box retains the shape of a parallelogram [49]. The simulation is thus performed at zero pressure, and a range of temperatures measured in units of ka 2 for up to N = 3600 sites. We use these simulations to investigate the phase diagram corresponding to the Hamiltonian (4.3) [50] and to investigate the properties of the phases we encounter, such as ground state degeneracy, order-by-disorder and negative thermal expansion. As all simulations involve a finite lattice and are run for a finite time, we took care to make sure that the system is sufficiently large to capture the thermodynamic behavior and that the simulation time is sufficiently long for the system to relax to equilibrium. To capture the subtlety of the order by disorder effect for a finite system, we simulated the model for a range of sizes and times and calculated the average fraction of zigzags n in equilibrium [inset of Fig. 3(b)]. While for small systems, φ ≈ 1 2 , as in a disordered zigzagging configuration, for large systems, φ approaches 0, suggesting that the system prefers the configuration of a uniformly sheared square lattice. Thus, we find good agreement with theoretical results from Sec. IV.
We calculated the shape of the phase boundary shown in Fig. 2, the behavior of the order parameter across the phase boundary shown in Fig. 4a, and the negative of the thermal expansion coefficient (L − L 0 )/L 0 shown in Fig. 4(b), where L is the length at T > 0 and L 0 that at T = 0. The phase boundary is obtained by calculating the heat capacity of the system as a function of temperature at fixed λ and τ . The location of the peak of the heat capacity corresponds to the location of the phase transition in the thermodynamic limit, and in our simulations, the locations of the peak converge to the values seen in Fig. 2. The order parameter values are calculated locally, i.e., t is calculated for each plaquette and each configuration via the angle between the two adjacent nearest-neighbor bonds and then averaged over all plaquettes in the system and over all 100 configurations. In this approach, t is independent of the particular zigzag configuration, and its evaluation does not exhibit the long relaxation process to the uniformly sheared square lattice. The behavior of the order parameter in Fig. 4 is consistent with a weakly-discontinuous transition.

VI. ANALYTIC THEORY AND THE PHASE DIAGRAM
The special feature of our model is its large but subextensive number of soft modes living on the q x and q y axes in the first Brillouin zone, as shown in Fig. 1(b) in the limit κ → 0) [22,51]. As we discuss below, these floppy modes provide a divergent fluctuation correction to the rigidity of the square lattice and render the transition from the square to the rhombic phase first order. Thus, this model is analogous to the one introduced by Brazovskii [39] for a liquid-solid transition in which mode frequencies of the form ω = ∆ + (q − q c ) 2 /m vanish on a (d − 1)-dimensional hypersphere when ∆ → 0 and contribute terms to the free energy singular in ∆.
To study the square-to-rhombic transition, we take the T = 0 square lattice, with site i at position R 0,i , as the reference state. We then represent positions in the distorted lattices as the sum of a part arising from uniform strain characterized by a deformation tensor Λ and deviations u ′ i from that uniform strain: The deviations u ′ i are constrained to satisfy periodic boundary conditions, and their average u ′ i is constrained to be zero. The former condition ensures that the sum over all bond stretches arising from these deviations vanishes for every configuration. Without loss of generality, we take Λ yx to be zero, leaving three independent parameters to parameterize the three independent strains. As we detail in Supplementary Information Sec. IV, the strain parameter characterizing pure shear with axes along the x and y axes of the reference lattice is of order t 2 and can be ignored near τ = 0, and we set Note that Λ is invertible even though it is not symmetric. t is the order parameter that distinguishes the rhombic phase from the square phase. Thermal fluctuations lead to s < 0 in both phases. Expanding the Hamiltonian [Eq. (4.3)] in to second order u ′ i about the homogeneously deformed state R i = Λ · R 0,i , we obtain where H 0 is the energy of the uniformly deformed state, and is the d × d dimensional (d = 2 being spatial dimension) dynamical matrix with scalar v q and second rank tensor M q determined by the potentials. There is no term linear in u ′ i in Eq. (6.3) because of the periodicity constraint (see Supplementary Information Sec. II).
Integrating out the fluctuations u ′ from the Hamiltonian (6.3), we obtain the free energy of the deformed state [52] where ǫ0 is the full nonlinear strain. Thus the one-loop free energy of Eq. (6.5) is a function of the nonlinear, rather than the linear strain, so that rotational invariance in the target space is guaranteed and there is a clean distinction between nonlinear terms in linearized deformations arising from nonlinearities in ǫ0 in from nonlinear terms in the expansion in powers of ǫ0.
To analyze the transition between the square and the rhombic phases at low temperature we expand F as a series in ǫ0, by expanding the transformed dynamical matrix asD = D 0 + A(ǫ0), where D 0 = D| ǫ0=0 is the dynamical matrix of the undeformed state. The free energy is then where G 0 ≡ D −1 0 is the phonon Green's function in the undeformed state, V = N a 2 and f is the free energy density. The expansion of F at small ǫ0 follows from this.
Close to the transition, F is dominated by fluctuations coming from the floppy modes as we discussed above. As κ → 0, the frequency of these floppy modes vanishes as ω ∼ √ κ, and the corresponding phonon Green's function diverges, leading to divergent fluctuation corrections to the coefficients of ǫ0 in Eq. (C4) as detailed in Supplementary Information Sec. III. Keeping leading order terms as τ → 0, we can identify the two phases through the equations of state, is a unitless reduced temperature. Eq. (6.8) has three solution for t: t = 0 corresponding to the square phase, and two solutions for t = 0 corresponding to the two orientations of the rhombic phase. There is only a single solution for s, with s < 0, from which we conclude that both phases exhibit negative thermal expansion. The elastic rigidity and thus the stability of the two phases is determined by the second derivatives of F with respect to s and t. In particular, the reduced shear modulus (G/k where G is the shear modulus) is (6.10) To obtain these leading order equations we (i) assume low T , so only terms singular in τ as τ → 0, such as t/ √ τ , in the integral of T 2 ln DetD(Λ), and (ii) assume that |s| ∼ t 2 ≪ τ ≪ 1, (6.11) the validity of which will be verified below.
As observed in the simulation (Fig. 2), thermal fluctuations at T > 0 stabilize the square relative to the rhombic phase even for κ < 0. To understand this phenomenon within the analytic approach, we use a self-consistentfield approximation in which τ is replaced by with its renormalized value r in the phonon Green's function G 0 and thus in the denominators on the right hand sides of Eqs. (6.7), (6.8), and (6.10). In this approximation, the shear rigidity of the square (t = 0) and rhombic (t = 0) satisfy where we used the equation of state, Eq. (6.8) to eliminate t 2 from Eq. (6.10). In the square phase, Eq. (6.12) has a solution r > 0, implying local stability, everywhere except atT = 0, τ < 0 and in the uninteresting limit, z → −∞. This local stability implies that the transition to the rhombic phase must be first order. In the rhombic phase, solutions r > 0 only exist for τ < τ c1 (z < z c1 ), where τ c1 = − 3 2 (T λ) 2/3 (6.13) The solutions to Eq.(6.12) can conveniently be expressed as scaling functions in the two phases: , (6.14) where ν = s, r for the square and rhombic phases, respectively. The scaling functions h s (|z|) and h r (|z|) depicted in Fig. 5 have the following limits 16) where z c1 = −3/2. The z 3 regimes of h s is in the metastable regime where the rhombic phase is stable. These results yield the scaling phase diagram of Fig. 2. The phase boundary of this discontinuous transition occurs along the coexistence line (i.e., equal free-energy line) of the two phases. Following Brazovskii, we have already calculated the limit of metastability of the rhombic phase, i.e., for the value of κ = τ c1 [Eq. (6.13)] at the local free energy minimum where that phase first appears. We then calculate the free energy difference between the two phases, which is evaluated through the following integral for a given τ , where we have substituted r(t) for the integral measure.
Here r 0 and r 1 are the values of r at the minima of F corresponding to the square and the rhombic phases determined by Eq. (6.12). Along the path of this integral, Eq. (6.10) is valid, but the equations of state (6.7) and (6.8) and the equation [Eq. (6.12)] for r in the rhombic phase are not satisfied, because they only apply to equilibrium states. The phase boundary corresponds to the curve of τ = τ c2 along which ∆F | τ =τc2 = 0. As shown in Supplementary Information Sec. IV, an asymptotic solution valid at low τ , can be obtained by expanding the equation around τ = τ c1 , assuming that τ 1 and τ 2 are of the same order of magnitude (verified below). This yields forT ≪ 1 where c ≃ 0.216 is a constant. This transition line is shown in Fig. 2. Excellent agreement between theory and simulation is obtained without any fitting parameter.
Along the phase boundary, r 0 , r 1 ∼T 2/3 > 0, so that both phases are locally stable. The order parameter for the transition, t, jumps from 0 to t c2 = 3.4T 1/3 λ −1/6 (6. 19) at the transition. As T → 0 this discontinuity vanishes, consistent with the continuous nature of the transition at T = 0. A good agreement between t values in theory and in simulation is shown in Fig. 4(a). From Eq. (6.7), the negative thermal expansion coefficient in both the square and rhombic phases, is determined by the equation of state, Eq. (6.7). Equation (6.12) for r then implies the following behavior for s in different regions of the phase diagram: In the critical region 0 < |τ | <T 2/3 of the square phase, Deep in the square and rhombic phases, where 0 < T 2/3 ≪ |τ |, Finally along the coexistence curve in both phases, s ∼ −T / |τ | ∼T 2/3 . These results agree well with simulation measurement of negative thermal expansion, as shown in Fig. 4b. In this lattice the negative thermal expansion behavior results from strong transverse fluctuations associated with soft modes. These solutions for s and t verifies that our assumptions in Eq. (6.11) are satisfied, provided that λ ≫ 1.

VII. REVIEW
We have presented an analysis of a model based on the square lattice with NN harmonic and NNN anharmonic springs that can be tuned at zero temperature from a stable square lattice through the mechanically unstable NN square lattice to a highly degenerate zigzag state by changing the coefficient κ of the harmonic term in the NNN spring from positive through zero to negative. Using analytic theory, including a generalization of the Brazovskii theory for the liquid-to-crystal transition, we investigated the phase diagram and mechanical properties of this model at T > 0. The degeneracy of the zero-T zigzag state is broken by an order-by-disorder effect, thermal fluctuations drive the square-to-rhombic phase transition first order, and the elastic modulus of the square phase a crossover from being proportional to κ for κ ≫ k(T g/k 2 ) 2/3 > 0 to T 2/3 for |κ| ≪ k(T g/k 2 ) 2/3 to T 2 for κ ≪ −k(T g/k 2 ) 2/3 < 0 as function of the scaling variable z = τ /(T λ) 2/3 =∼ (κ/k)/(T g/k 2 ) 2/3 . This behavior arises because the spectrum of the NN square lattice with N sites exhibits √ N zero modes on a onedimensional manifold in the Brillouin Zone. Other lattices such as the 2D kagome lattice, the 3D simple cubic lattice, and the 3D pyrochlore and β-cristobalite [53] lattices have similar spectra, and it is our expectation that generalizations of our model to these lattices will exhibit similar behavior. It is also likely that our model can inform us about more physically realistic models in which interactions lead to spectra with a large set of modes with small but not zero frequency.
Acknowledgments -A.S. gatefully acknowledges discussions with P. A. Rikvold, Gregory Brown, Shengnan Huang and Andrea J. Liu The minimum of V PL in Eq. (4.5) in the main text can be found analytically by assuming a rhombic plaquette and solving the set of equations given by ∂V PL /∂d 1 = 0 and ∂V PL /∂d 2 = 0 for τ and λ. These equations are linear, and may be written in terms of the plaquette side b and the inner angle α, where d 1,2 = √ 2b √ 1 ± sin α. The solutions have the form It was verified numerically that the rhombic plaquette with side length b minimizes the potential energy for the range of parameters considered in this work relative to plaquettes with sides of unequal length.

Appendix B: Expansion of lattice Hamiltonian at deformed reference states
For a generic lattice with pair-wise potentials, we can write the Hamiltonian as where b labels bonds {i, j}, V b is the interaction potential of the bond, and are the bond vectors in the reference and target states. We consider a macroscopic deformation Λ (corresponding to ǫ 0 in the continuum theory) and fluctuations u ′ , so the deformation of an arbitrary site can be written as Thus for a bond, The change of bond length can then be expanded for small The terms we kept are exact to O(ũ 2 ). The expansion of the potential of a bond is then where V bΛ , V ′ bΛ , V ′′ bΛ are the potential and its derivatives at the macroscopic deformation value We can then sort terms in powers ofũ ′ Therefore the total Hamiltonian of the lattice is where H 0 (Λ) is the energy for reference state with the uniform deformation (Λ), and is the additional potential energy coming from fluctuations around the uniformly deformed reference state. Here we express it in momentum space, v 0 is the area of the unit cell, and and B represent the sum over all bonds in one unit cell. From this we can write the dynamical matrix D as a sum of two parts D q (Λ) = v q (Λ)I + Λ · M q (Λ)Λ T .
We can then calculate the expression for D q (Λ) for the special case of the square lattice model with NN harmonic springs of spring constant k and NNN springs with the potential where g starts from O(ǫ 0 ) and M starts from O(1). This confirms the Ward identity in this problem.
To analyze the transition we then expand F as a series of ǫ 0 . Because the dynamical matrix can be expanded as where D 0 = D| ǫ0=0 is the dynamical matrix of the undeformed state. The free energy is then where G 0 ≡ D −1 0 is the phonon Green's function in the undeformed state. The expansion of F at small ǫ 0 thus follows from this, The above free energy can be calculated by performing integrals in momentum space. Because we are interested in characterizing the square-rhombic phase transition, where τ ≡ κ/k is small, we can expand our results at small τ . In this limit, floppy modes of frequency √ τ lie along q x and q y axes, and to lowest order in τ , the Green's function can be approximated by Integrals involving this Green's functions can be evaluated following the calculation discussed in Ref. [40] (add the kagome lattice ref here which contain more detail). We also verified that to leading order in small τ , integrals done using this approximation agree with exact results.