Mechanical instability at finite temperature
Xiaoming Mao, Anton Souslov, Carlos I. Mendoza, T. C. Lubensky
aa r X i v : . [ c ond - m a t . s o f t ] J u l Mechanical instability at finite temperature
Xiaoming Mao, Anton Souslov, Carlos I. Mendoza, and T. C. Lubensky Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332, USA Instituto de Investigaciones en Materiales, Universidad Nacional Aut´onoma de M´exico,Apdo. Postal 70-360, 04510 M´exico, D.F., Mexico Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA (Dated: October 12, 2018)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 insta-bility. Their properties are determined by an interplay between soft mechanical modes and thermalfluctuations. In this paper we investigate a simple square-lattice model with a φ potential betweennext-nearest-neighbor sites whose quadratic coefficient κ can be tuned from positive negative. Weshow that its zero-temperature ground state for κ < κ < κ = 0 leads to singularcontributions to the free energy that render the square-to-rhombic transition first order and leadto power-law behavior of the shear modulus as a function of temperature. We expect our study toprovide a general framework for the study of finite-temperature mechanical and phase behavior ofother systems with a large number of floppy modes. I. INTRODUCTION
Crystalline solids can undergo structural phase transi-tions in which there is a spontaneous change in the shapeor internal geometry of their unit cells [1–4]. These tran-sitions are signalled by the softening of certain elasticmoduli or of phonon modes at a discrete set of pointsin the Brillouin zone. Lattices with coordination num-ber z = z c = 2 d in d spatial dimensions, which we willcall Maxwell lattices [5][6], exist at the edge of mechani-cal instability, and they are critical to the understandingof systems as diverse as engineering structures [7, 8], di-luted lattices near the rigidity threshold [9–11], jammedsystems [12–14], biopolymer networks [15–19], and net-work glasses [20, 21]. Hypercubic lattices in d dimensionsand the kagome lattice and its generalization to higherdimensions with nearest-neighbor (NN) Hookean springsof spring constant k are a special type of Maxwell latticewhose phonon spectra have harmonic-level zero modesnot at a discrete set of points but at all N ( d − /d pointson ( d − κ tuned from positive to negative and aquartic stabilizing term. When κ >
0, the square lat-tice is stable even at zero temperature. When κ = 0, NNN springs contribute only at anharmonic order, andthe harmonic phonon spectrum is identical to that of theNN-lattice. When κ <
0, the NNN potential has twominima, and the ground state of an individual plaque-tte is a rhombus that can have any orientation (Fig. 1c).Plaquettes in the same row (or column), however, areconstrained to have the same orientation, but plaquettesin adjacent rows can either tilt in the same direction or inthe mirror-image direction as shown in Figs. 1d-f, leadingto 2 × √ N equivalent ground states and a subextensivebut divergent entropy of order √ N ln 2. The propertiesof this model, including the subextensive entropy at zerotemperature, are very similar to those of colloidal par-ticles confined to a low-height cells [23, 24] and to theanti-ferromagnetic Ising model on a deformable trian-gular lattice [25]. In addition, the scaling of the shearmodulus near the zero-temperature critical point is anal-ogous to that observed in finite-temperature simulationsof randomly diluted lattices near the rigidity-percolationthreshold [26] and to finite-temperature scaling near thejamming transition [27], suggesting that generalizationsof our model and approach may provide useful insightinto the thermal properties of other systems near theMaxwell rigidity limit. II. RESULTS
Strong fluctuations arising from the large number ofzero modes lead to interesting physics at
T > • Among all the equal-energy zigzagging configura-tions at κ <
0, the uniformly sheared rhombiclattice, shown in Fig. 1d, has the lowest free en-ergy: the ground-state degeneracy is broken by
FIG. 1. (a) The square lattice model with NN (blue thick)bonds and NNN (brown thin) bonds. White disks showinga shift of the second row is one example of a floppy modesof the lattice with no NNN bond. (b) Density plot of thephonon spectrum of the NN square lattice showing lines ofzero modes (darker color corresponds to lower frequency). (c)The T = 0 ground states of a plaquette when κ <
0, with theundeformed reference state shown in gray. (d-f) Examples of T = 0 ground states of the whole lattice when κ <
0. (d)shows the uniformly sheared rhombic lattice, which we showto be the preferred configuration at small T in the thermody-namic limit. (e) is a randomly zigzagging configuration, and(f) is the ordered maximally zigzagging configuration, whichhas a unit cell consisting of two particles. thermal fluctuations through an order-by-disordereffect [25, 28–34]. • Thermal fluctuations lead to a negative coefficientof thermal expansion and corrections to the shearrigidity [35–38] that enable the square lattice stateto remain thermodynamically stable in a region ofthe phase diagram at κ ≤ • Fluctuations drive the transition from the rhombicto the square phase first order and lead to the phasediagram shown in Fig. 2 in which the temperatureof the transition approaches zero as κ → − . • The low- T shear modulus G (Fig. 2) in both phasesis proportional to | κ | at low temperature, and thereis a critical regime with T > const. × κ / in which G ∼ T / . In addition, there is a region in thesquare phase (mostly metastable with respect tothe rhombic phase) with T < const. × | κ | / inwhich G ∼ ( T / | κ | ) . This behavior near the T = 0critical point is analogous to that found in therandomly diluted triangular lattice [26] near thecentral-force rigidity threshold. Interestingly, thecritical regime in our model is fundamentally a con-sequence of nonlinearity as is the case for dynamicalscaling near the jamming transition [27]. FIG. 2. (a) An example of the phase diagram of the modelsquare lattice at k = 1, g = 10, a = 1. The black solid line(the red dots) shows the boundary obtained from analytictheory (Monte Carlo simulation) between the square phaseon the right and the rhombic phase on the left. The squarephase is stabilized by thermal fluctuations even for κ < G of the square lattice exhibit differentscaling regimes (separated by black dashed lines) determinedby Eq. (6.16). (b) and (c) show Monte Carlo snapshots ofthe square and rhombic phases respectively. A small numberof zigzags exist in the rhombic phase Monte Carlo snapshots,resulting from finite size effects as discussed in Sec. IV. These predictions are supported by our Monte-Carlosimulations and by direct calculations of entropic contri-butions to the free energy phonon fluctuations in differentarrangements of kinks.
III. DISCUSSION
Our model of the square-to-rhomic transition is verysimilar to a model, studied by Brazovskii [39], for thetransition from an isotropic fluid to a crystal and laterapplied to the Rayleigh-B´enard instability [40] and tothe nematic-to-smectic- C transitions in liquid crystals[41, 42]. In all of these systems, there is a subexten-sive but infinite manifold of zero modes [a ( d − d − d − κ < κ > κ <
0, and it is instructiveto review their origin. A powerful index theorem [43]relates the number of zero modes of a frame consisting of N points and N B ≡ zN bonds, where z is the averagecoordination number, via the relation N = dN − N B + S, (3.1)where S is the number of independent states of self-stressin which bonds are under tension or compression and inwhich the net force on each point is zero. In his seminal1864 paper [5], Maxwell considered the case with S = 0that yields the Maxwell relation for the critical coordi-nation number at which N is equal to the number n ( d )(= d for periodic and d ( d + 1) / z Nc = 2 d − n ( d ) N . (3.2)In the limit of large N , z Nc → z ∞ c = 2 d .There are many small unit cell periodic Maxwell lat-tices with N = S . The NN square and kagome latticesin two dimensions and, the cubic and pyrochlore latticesin three dimensions are special examples of these lat-tices that have sample-spanning straight lines of bondsthat support states of self-stress under periodic bound-ary conditions. They, therefore, have of order N ( d − /d states of self-stress and the same number of zero modes,which are indicators of buckling instabilities of the lineswhen subjected to compression. Geometrical distortionsof theses lattices that remove straight lines, as is the casewith the twisted kagome lattice [44], remove states of self-stress and associated zero modes. When subjected to freerather than periodic boundary conditions, these distortedlattices have a deficiency of order N ( d − /d bonds and asa result the same number of surface zero modes (thereare no bulk zero mode other than the trivial ones of uni-form translation), which can have a topological charac-ter [45] or be described in the long-wavelength limit by aconformal field theory [44]. Unlike the infinitesimal zeromodes of hypercubic lattices, those of the kagome andpyrochlore do not translate into finite zero modes of thelattices when finite sections are cut from a lattice underperiodic boundary conditions. Thus, it is not yet clearwhether the ground state of the latter lattices will arehighly degenerate or not. Nevertheless, Brasovskii theoryshould provide a sound description of thermal propertiesof these lattices in the vicinity of the point T = 0, κ = 0. IV. MODEL AND ORDER-BY-DISORDER INTHE LOW-SYMMETRY PHASE
The model we consider is a square lattice with two dif-ferent types of springs – those connecting nearest neigh-bors and those connecting next-nearest neighbors, asshown in Figure 1(a). The NN springs are Hookian, withpotential V NN ( x ) = k x (4.1)where k >
0. The NNN springs are introduced with ananharmonic potential V NNN ( x ) = κ x + g x , (4.2)where κ can be either positive or negative and g , intro-duced for stability, is always positive. The Hamiltonianof the whole lattice is thus, H = X h i,j i∈ NN V NN (cid:16) | ~R i − ~R j | − a (cid:17) + X h i,j i∈ NNN V NNN (cid:16) | ~R i − ~R j | − √ a (cid:17) , (4.3)where ~R i is the positions of the node i and a is the lat-tice constant. In what follows, we will use the reducedvariables τ ≡ κ/k and λ ≡ ga /k (4.4)to measure the strength of couplings in V NNN .For κ > V NNN ( x ) has a unique minimum at x = 0,and the ground state of H is the square lattice withlattice spacing a . All elastic moduli of this state arenonzero, and it is stable with respect to thermal fluc-tuations, though, as we shall see, it does undergo ther-mal contraction at nonzero temperature. When κ < V NNN ( x ) has two minima at x = ± p κ/g , correspondingto stretch and compression, respectively. This change inlength of NNN springs is resisted by the NN springs, andin minimum energy configurations one NNN bond in eachplaquette will stretch and the other will contract. Thealternative of having both stretch or contract would costtoo much NN energy. A reasonable assumption, whichis checked by our direct calculation, is that the stretch-ing and contraction will occur symmetrically about thecenter so that the resulting equilibrium shape is a rhom-bus rather than a more general quadrilateral (see Sup-plementary Information Sec. I). The shape of a rhombusis uniquely specified by the lengths d and d of its di-agonals (which are perpendicular to each other), whoseequilibrium values are obtained by minimizing the sumover plaquettes of V PL = 2 V NN ( 12 q d + d − a ) + V NNN ( d − √ a ) + V NNN ( d − √ a ) . (4.5)The grounds state of the entire lattice when κ < ǫ , is simply V PL evaluatedat the equilibrium values of d and d .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 verticalzigzags. Thus, the ground state entropy diverges in thethermodynamic limit, though it is sub-extensive and pro-portional to N x + N y in a system of N = N x N y particles.Such ground state configurations are found in other sys-tems, most notably the zigzagging phases seen in suspen-sions of confined colloidal particles [23, 24]. The confinedcolloidal system has a phase diagram that depends onlyon the planar density and the height of confinement of thecolloids. For sufficiently large heights, the colloids forma phase of two stacked square lattices. In a neighbor-ing region of the phase diagram, explored in simulationsof Refs. [46, 47], this square lattice symmetry is brokenthrough a weakly discontinuous transition and a rhombicphase is observed. This region of the phase diagram ofconfined colloids thus provides a physical realization ofthe Hamiltonian (4.3).At low but nonzero temperatures, the degeneracyof the ground state is broken by thermal fluctuationsthrough the order-by-disorder mechanism [25, 28, 30–34].This splitting of degeneracy due to small phonon fluctu-ations around the ground state may be calculated usingthe dynamical matrix in the harmonic approximation tothe Hamiltonian (4.3). For each ground state configura-tion ~R i , we write the deformation as ~R i → ~R i + ~u i , andexpand to quadratic order in ~u , H = P h ij i ~u i D ij ~u j .The Fourier transform of D ij , D q , is block-diagonal andthe phonon free energy, which is purely entropic, for thatconfiguration is F p ( N x , N y , K ) = 12 k B T X q ln det D q ≡ N k B T w p , (4.6)where w p is the free-energy per site in units of k B T . Ingeneral, F p depends not only on K , but also on the par-ticular sequence of zigzags. We numerically calculatedthis free energy for all periodically zigzagged configura-tions with up to 10 sites per unit cell. We found thatthe lowest-free-energy state is the uniformly sheared state[Fig. 1(c)] with K = 0, and the highest-free-energy stateis 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∆ F ( N x , N y ) = F p ( N x , N y , N y ) − F p ( N x , N y , ≡ N x N y k b T ∆ w (4.7)∆ F ( N x , N y , K ) = F p ( N x , N y , K ) − F p ( N x , N y , . Figure 3(a) displays our calculation of ∆ w ( τ ), whichvanishes, as expected, at τ = 0 and also at large τ atwhich the rhombus collapses to a line. Figure 3(b) plots∆ F/ ∆ F as a function of φ = K/N y for different values of τ . By construction, this function must vanish at φ = 0and be equal to one at φ = 1. All of the points lieapproximately on a straight line of slope one. Thus, wecan approximate F p ( N x , N y , K ) by F p ( N x , N y , K ) = F p ( N x , N y ,
0) + N x Kk B T ∆ w . (4.8)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 thermody-namics limit, N →
0, in which the sum over q is replacedby the continuum limit by an integral. In order to com-pare these results with the Monte Carlo results of thenext section, it is necessary to study finite size effects.The first observation is that for any finite N x , the sys-tem will be effectively one-dimensional for a sufficientlylarge N y , and as a result, we would expect the numberof zigzags to fluctuate. To proceed, we continue to usethe continuum limit to evaluate Eq. (4.6), and we useEq. (4.8) for F p ( N x , N y , K ). The partition function forthis energy is Z = N y X K =0 e − β [ F p ( N x ,N y , E ( N x ,N y )] (cid:18) N y K (cid:19) e − N x K ∆ w = e − β [ F p ( N x ,N y , E ( N x ,N y )] (cid:0) e − N x ∆ w (cid:1) N y (4.9)where E ( N x , N y ) denotes the potential energy which isindependent of K , and the full free energy is F ( N x , N y , T ) = − T ln Z = N x N y ( f + e ) − N y k B T ln(1 + e − N x ∆ w ) , (4.10)where e and f are the potential energy and phonon freeenergy per site of the uniformly sheared state. Thus, asexpected, when N x ∆ w ≫
1, zigzag configurations makeonly a very small, subextensive contribution to the freeenergy. On the other hand, in the opposite limit, theymake an extensive contribution of N x N y k B T ∆ w to theenergy. Therefore, at a given τ , the zigzag configurationsare favored when the system is small, and in thermody-namic limit, the rhombic configuration is always favored.Our Monte Carlo simulation verified this [see inset ofFig. 3(b)]. V. SIMULATION
We simulate the system using a Monte Carlo (MC)algorithm inside a periodic box whose shape and size | | τ / F F ∆ ∆ φ / B F Nk T ∆ ab φ FIG. 3. The phonon contribution to the free energy for var-ious zigzagging configurations. (a) shows the free energy dif-ference ∆ w between the maximally zigzagging configuration(Fig. 1e) and the uniformly sheared square lattice configura-tion (Fig. 1c). For sufficiently large τ , the lattice collapsesonto a line, at which point the free energy difference goes tozero. (b) shows the free energy of the phonons as a ratio∆ F/ ∆ F for τ = − . − . − λ = 10), where ∆ F and ∆ F are defined in Eq. (4.7) and ∆ F is evaluated for allpossible configurations with unit cells with at most 10 sitesas a function of the zigzag fraction φ = K/N x . The latticewithout zigzags is entropically favored for any any value of τ in our calculations. ∆ F/ ∆ F is well approximated by the line φ , which corresponds to non-interacting zigzags. This interac-tion results in the dispersion of values of ∆ F/ ∆ F . The insetof (b) shows Monte Carlo simulation results for the zigzaggingfraction φ plotted against the linear system size √ N . are allowed to change in order to maintain zero pres-sure. In this version of the Metropolis algorithm, alsoused in Refs. [25, 48], for each MC step a particle ispicked at random and a random trial displacement isperformed. The trial displacement is initially uniformlydistributed within a radius of 0 . a , but throughout thesimulation the radius is adjusted to keep the acceptanceprobability between 0 .
35 and 0 .
45. Given the initial con-figuration energy E i and the trial configuration energy E j , the trial configuration is accepted with probability[1 + exp( E i − E j ) /T ] − , i.e., Glauber dynamics is used.After initializing the system using the square lattice con-figuration with lattice constant a , the simulation is firstrun at a high temperature and is then annealed to the fi-nal low temperature. For each intermediate temperature,an equilibration cycle in a sample of N sites consists ofat least N × MC steps. To accommodate areal andshear distortions in the different phases we encounter,the simulation box area and shape are changed using asimilar acceptance algorithm, with the trial deformationadjusted to keep the acceptance probability between 0 .
35 and 0 .
45, such that the simulation box retains the shapeof a parallelogram [49]. The simulation is thus performedat zero pressure, and a range of temperatures measuredin units of ka for up to N = 3600 sites.We use these simulations to investigate the phase di-agram corresponding to the Hamiltonian (4.3) [50] andto investigate the properties of the phases we encounter,such as ground state degeneracy, order-by-disorder andnegative thermal expansion. As all simulations involvea finite lattice and are run for a finite time, we tookcare to make sure that the system is sufficiently large tocapture the thermodynamic behavior and that the simu-lation time is sufficiently long for the system to relax toequilibrium. To capture the subtlety of the order by dis-order effect for a finite system, we simulated the modelfor a range of sizes and times and calculated the averagefraction of zigzags n in equilibrium [inset of Fig. 3(b)].While for small systems, φ ≈ , as in a disordered zigzag-ging configuration, for large systems, φ approaches 0,suggesting that the system prefers the configuration ofa uniformly sheared square lattice. Thus, we find goodagreement with theoretical results from Sec. IV.We calculated the shape of the phase boundary shownin Fig. 2, the behavior of the order parameter across thephase boundary shown in Fig. 4a, and the negative ofthe thermal expansion coefficient ( L − L ) /L shown inFig. 4(b), where L is the length at T > L that at T = 0. The phase boundary is obtained by calculatingthe heat capacity of the system as a function of temper-ature at fixed λ and τ . The location of the peak of theheat capacity corresponds to the location of the phasetransition in the thermodynamic limit, and in our simu-lations, the locations of the peak converge to the valuesseen in Fig. 2. The order parameter values are calcu-lated locally, i.e., t is calculated for each plaquette andeach configuration via the angle between the two adja-cent nearest-neighbor bonds and then averaged over allplaquettes in the system and over all 100 configurations.In this approach, t is independent of the particular zigzagconfiguration, and its evaluation does not exhibit the longrelaxation process to the uniformly sheared square lat-tice. The behavior of the order parameter in Fig. 4 isconsistent with a weakly-discontinuous transition. VI. ANALYTIC THEORY AND THE PHASEDIAGRAM
The special feature of our model is its large but subex-tensive number of soft modes living on the q x and q y axesin the first Brillouin zone, as shown in Fig. 1(b) in thelimit κ →
0) [22, 51]. As we discuss below, these floppymodes provide a divergent fluctuation correction to therigidity of the square lattice and render the transitionfrom the square to the rhombic phase first order. Thus,this model is analogous to the one introduced by Bra-zovskii [39] for a liquid-solid transition in which modefrequencies of the form ω = ∆ + ( q − q c ) /m vanish on -0.20 -0.15 -0.10 -0.05 0.000.00.10.20.30.40.5 ab T κ t − − − − − − () / LLL − FIG. 4. (a) Order parameter t calculated from theory (lines)and simulation (data points), at λ = 10, and from right toleft, T / ( ka ) = 0 . , . , . , . , . L − L ) /L as a function of T at κ = 0 , g = 10(red triangles: MC data; red upper line: theory), and at κ =0 . , g = 0 (black circles: MC data; black lower line: theory). a ( d − → T = 0 square lattice, with site i at position ~R ,i , as thereference state. We then represent positions in the dis-torted lattices as the sum of a part arising from uniformstrain characterized by a deformation tensor Λ and devi-ations ~u ′ i from that uniform strain: ~R i = Λ · ~R ,i + ~u ′ i , (6.1)The deviations ~u ′ i are constrained to satisfy periodicboundary conditions, and their average h ~u ′ i i is con-strained to be zero. The former condition ensures thatthe sum over all bond stretches arising from these devi-ations vanishes for every configuration. Without loss ofgenerality, we take Λ yx to be zero, leaving three indepen-dent parameters to parameterize the three independentstrains. As we detail in Supplementary Information Sec.IV, the strain parameter characterizing pure shear withaxes along the x and y axes of the reference lattice is oforder t and can be ignored near τ = 0, and we set Λ = (cid:18) s t s (cid:19) . (6.2)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 leadto s < ~u ′ i about the homogeneously deformed state h ~R i i = Λ · ~R ,i , we obtain H = H ( Λ ) + 1 V X q ~u ′ q · D q ( Λ ) · ~u ′− q + O (( ~u ′ ) ) , (6.3)where H is the energy of the uniformly deformed state,and D q ( Λ ) = v q ( Λ T Λ ) I + Λ · M q ( Λ T Λ ) · Λ T (6.4)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 linearin ~u ′ i in Eq. (6.3) because of the periodicity constraint(see Supplementary Information Sec. II).Integrating out the fluctuations ~u ′ from the Hamilto-nian (6.3), we obtain the free energy of the deformedstate [52] F ( Λ ) = H ( Λ T Λ ) + T D ( Λ ) , (6.5)where ˜ D = v q + M q Λ T Λ depends only on Λ T Λ = 1+2 ǫ ,where ǫ is the full nonlinear strain. Thus the one-loopfree energy of Eq. (6.5) is a function of the nonlinear,rather than the linear strain, so that rotational invariancein the target space is guaranteed and there is a clean dis-tinction between nonlinear terms in linearized deforma-tions arising from nonlinearities in ǫ in from nonlinearterms in the expansion in powers of ǫ .To analyze the transition between the square and therhombic phases at low temperature we expand F as aseries in ǫ , by expanding the transformed dynamicalmatrix as ˜ D = D + A ( ǫ ), where D = D | ǫ =0 is thedynamical matrix of the undeformed state. The free en-ergy is then F ( ǫ ) = H ( ǫ ) + T D [ I + G · A ( ǫ )] ≡ V f ( ǫ ) , (6.6)where G ≡ D − is the phonon Green’s function in theundeformed state, V = N a and f is the free energydensity. The expansion of F at small ǫ follows fromthis.Close to the transition, F is dominated by fluctuationscoming from the floppy modes as we discussed above.As κ →
0, the frequency of these floppy modes vanishesas ω ∼ √ κ , and the corresponding phonon Green’s func-tion diverges, leading to divergent fluctuation correctionsto the coefficients of ǫ τ →
0, we can identify the two phases throughthe equations of state, ∂f ( ǫ ) ∂s ≃ ka ˜ T √ τ + s ! = 0 (6.7) ∂f ( ǫ ) ∂t ≃ ka t τ + λ ˜ T √ τ + 112 λt ! = 0 , (6.8)where ˜ T = πT ka (6.9)is a unitless reduced temperature. Eq. (6.8) has three so-lution for t : t = 0 corresponding to the square phase, andtwo solutions for t = 0 corresponding to the two orienta-tions of the rhombic phase. There is only a single solu-tion for s , with s <
0, from which we conclude that bothphases exhibit negative thermal expansion . The elasticrigidity and thus the stability of the two phases is de-termined 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 r = 1 k ∂ f ( ǫ ) ∂t ≃ τ + ˜ T √ τ + 14 λt . (6.10)To obtain these leading order equations we (i) assumelow T , so only terms singular in τ as τ →
0, such as t/ √ τ , in the integral of T ln Det ˜ D ( Λ ), and (ii) assumethat | s | ∼ t ≪ τ ≪ , (6.11)the validity of which will be verified below.As observed in the simulation (Fig. 2), thermal fluctua-tions at T > κ <
0. To understand this phenomenonwithin the analytic approach, we use a self-consistent-field approximation in which τ is replaced by with itsrenormalized value r in the phonon Green’s function G and thus in the denominators on the right hand sides ofEqs. (6.7), (6.8), and (6.10). In this approximation, theshear rigidity of the square ( t = 0) and rhombic ( t = 0)satisfy r = ( τ + ˜ T λ √ r square − τ − ˜ T λ √ r rhombic , (6.12)where we used the equation of state, Eq. (6.8) to elimi-nate t from Eq. (6.10). In the square phase, Eq. (6.12)has a solution r >
0, implying local stability, everywhereexcept at ˜ T = 0 , τ < z → −∞ . This local stability implies that the transitionto the rhombic phase must be first order. In the rhombicphase, solutions r > τ < τ c ( z < z c ),where τ c = −
32 ( ˜
T λ ) / (6.13) -4 -2 0 2 4 z h h s h r z c1 z c2 FIG. 5. (color online) Plots of h s ( | z | ) (black) and h r ( | z | ) (red)as a function of z . Note the singular behavior of h r ( | z | in thevicinity of z c and the large difference between h s and h r atthe first order transition at z = z c . The solutions to Eq.(6.12) can conveniently be ex-pressed as scaling functions in the two phases: | τ | r = h ν ( | z | ) , z = τ ( ˜ T λ ) / , (6.14)where ν = s, r for the square and rhombic phases, re-spectively. The scaling functions h s ( | z | ) and h r ( | z | ) de-picted in Fig. 5 have the following limits h s ( | z | ) ∼ z → + ∞| z | z → | z | z → −∞ (6.15) h r ( | z | ) ∼ ( (3 / − p | z − z c | z → z − c / z → z → −∞ , (6.16)where z c = − /
2. The z regimes of h s is in themetastable regime where the rhombic phase is stable.These results yield the scaling phase diagram of Fig. 2.The phase boundary of this discontinuous transitionoccurs along the coexistence line (i.e., equal free-energyline) of the two phases. Following Brazovskii, we have al-ready calculated the limit of metastability of the rhombicphase, i.e., for the value of κ = τ c [Eq. (6.13)] at the localfree energy minimum where that phase first appears. Wethen calculate the free energy difference between the twophases, which is evaluated through the following integralfor a given τ ,∆ F = Z t rhombic dF ( ǫ ) dt dt = Z r r dF ( ǫ ) dt dtdr dr, (6.17)where we have substituted r ( t ) for the integral measure.Here r and r are the values of r at the minima of F corresponding to the square and the rhombic phases de-termined 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 rhombicphase are not satisfied, because they only apply to equi-librium states. The phase boundary corresponds to thecurve of τ = τ c along which ∆ F | τ = τ c = 0. As shown inSupplementary Information Sec. IV, an asymptotic so-lution valid at low τ , can be obtained by expanding theequation around τ = τ c , assuming that τ and τ are ofthe same order of magnitude (verified below). This yields τ c = − (cid:18)
32 + c (cid:19) ( ˜ T λ ) / = − .
716 ˜ T / , (6.18)for ˜ T ≪ c ≃ .
216 is a constant. This transi-tion line is shown in Fig. 2. Excellent agreement betweentheory and simulation is obtained without any fitting pa-rameter.Along the phase boundary, r , r ∼ ˜ T / >
0, so thatboth phases are locally stable. The order parameter forthe transition, t , jumps from 0 to t c = 3 . T / λ − / (6.19)at the transition. As T → T = 0. A good agreement between t values in theoryand in simulation is shown in Fig. 4(a).From Eq. (6.7), the negative thermal expansion coeffi-cient in both the square and rhombic phases, s = − ˜ T √ r , (6.20)is determined by the equation of state, Eq. (6.7). Equa-tion (6.12) for r then implies the following behavior for s in different regions of the phase diagram: In the criticalregion 0 < | τ | < ˜ T / of the square phase, s ≃ − ˜ T / λ − / . (6.21)Deep in the square and rhombic phases, where 0 < ˜ T / ≪ | τ | , s = ( − ˜ T / p | τ | square − ˜ T / p | τ | rhombic . (6.22) Finally along the coexistence curve in both phases, s ∼− ˜ T / p | τ | ∼ ˜ T / . These results agree well with sim-ulation measurement of negative thermal expansion, asshown in Fig. 4b. In this lattice the negative thermalexpansion behavior results from strong transverse fluctu-ations associated with soft modes.These solutions for s and t verifies that our assump-tions in Eq. (6.11) are satisfied, provided that λ ≫ VII. REVIEW
We have presented an analysis of a model based onthe square lattice with NN harmonic and NNN anhar-monic springs that can be tuned at zero temperaturefrom a stable square lattice through the mechanicallyunstable NN square lattice to a highly degenerate zigzagstate by changing the coefficient κ of the harmonic termin the NNN spring from positive through zero to neg-ative. Using analytic theory, including a generalizationof the Brazovskii theory for the liquid-to-crystal transi-tion, we investigated the phase diagram and mechanicalproperties of this model at T >
0. The degeneracy ofthe zero- T zigzag state is broken by an order-by-disordereffect, thermal fluctuations drive the square-to-rhombicphase transition first order, and the elastic modulus ofthe square phase a crossover from being proportional to κ for κ ≫ k ( T g/k ) / > T / for | κ | ≪ k ( T g/k ) / to T for κ ≪ − k ( T g/k ) / < z = τ / ( ˜ T λ ) / = ∼ ( κ/k ) / ( T g/k ) / . Thisbehavior arises because the spectrum of the NN squarelattice with N sites exhibits √ N zero modes on a one-dimensional manifold in the Brillouin Zone. Other lat-tices such as the 2D kagome lattice, the 3D simple cubiclattice, and the 3D pyrochlore and β -cristobalite [53] lat-tices have similar spectra, and it is our expectation thatgeneralizations of our model to these lattices will exhibitsimilar behavior. It is also likely that our model can in-form us about more physically realistic models in whichinteractions lead to spectra with a large set of modes withsmall but not zero frequency. Acknowledgments –
A.S. gatefully acknowledges dis-cussions with P. A. Rikvold, Gregory Brown, ShengnanHuang and Andrea J. Liu. This work was supported inpart by the NSF under grants DMR-1104707 and DMR-1120901 (TCL), DMR-1207026 (AS), grant DGAPA IN-110613 (CIM) and by the Georgia Institute of Technology(AS). [1] Folk, R., Iro, H. & Schwabl, F. Critical elastic phasetranstions.
Z. Physik B , 69–81 (1976).[2] Cowley, R. A. Structural phase-transitions i: Landautheory. Advances in Physics , 1–110 (1980).[3] Bruce, A. Structural phase transitions ii: Static criticalbehaviour. Advances in Physics , 111–217 (1980).[4] Fujimoto, M. The Physics of Structural Phase Transi- tions (Springer, 2005), 2nd edn.[5] Maxwell, J. C. On the calculation of the equilibrium andstiffness of frames.
Philos. Mag. , 294 (1864).[6] The term isostatic is often incorrectly used to describeany system with z = z c . Finite Isostatic lattices have z = z Nc = z c − d ( d + 1) /N and no states of self stress- See [43] for example. There is no universally accepted definition of isostatic in lattices with periodic boundaryconditions, but one such as the square lattice with manystates of self-stress is surely not isostatic.[7] Heyman, J. The Science of Structural Engineering (Cen-gage Learning, Stamford CT, 2005).[8] Kassimali, A.
Structural Analysis (Cengage Learning,Stamford CT, 2005).[9] Feng, S. & Sen, P. N. Percolation on elastic networks:New exponent and threshold.
Phys. Rev. Lett. , 216–219 (1984).[10] Feng, S., Sen, P. N., Halperin, B. I. & Lobb, C. J. Per-colation on two-dimensional elastic networks with rota-tionally invariant bond-bending forces. Phys. Rev. B ,5386–5389 (1984).[11] Jacobs, D. J. & Thorpe, M. F. Generic rigidity percola-tion: The pebble game. Phys. Rev. Lett. , 4051–4054(1995).[12] Liu, A. J. & Nagel, S. R. Jamming is not just cool anymore. Nature , 21 (1998).[13] Wyart, M. On the rigidity of amorphous solids.
Ann.Phys. Fr , 1–96 (2005).[14] Liu, A. J. & Nagel, S. R. The Jamming Transition andthe Marginally Jammed Solid , vol. 1 of
Annual Review ofCondensed Matter Physics , 347–369 (2010).[15] Elson, E. L. Cellular mechanism as an indicator of cskstructure and function.
Annu. Rev. Biophys. Biophys.Chem , 397–430 (1988).[16] Kasza, K. et al. The cell as a material.
Curr. Opin. CellBiol. , 101–107 (2007).[17] Alberts, B. et al. Molecular Biology of the Cell (Garland,New York, 2008), 4th edn.[18] Janmey, P. et al. Resemblance of actin-binding actin gelsto covanlently cross-linked networks.
Nature (1990).[19] Broedersz, C. P. & MacKintosh, F. C. Modeling semi-flexible polymer networks. arXiv:1404.4332 (2014).[20] Phillips, J. C. Topology of covalent non-crystalline solids.2. medium-range order in chalcogenide alloys and a-si(ge).
J. Non-Cryst. Solids , 37–77 (1981).[21] Thorpe, M. Continuous deformations in random net-works. Journal of Non-Crystalline Solids , 355 – 370(1983).[22] Souslov, A., Liu, A. J. & Lubensky, T. C. Elasticity andresponse in nearly isostatic periodic lattices. PhysicalReview Letters , 205503 (2009).[23] Pieranski, P., Strzelecki, L. & Pansu, B. Thin colloidalcrystals.
Physical Review Letters , 900–903 (1983).[24] Han, Y. et al. Geometric frustration in buckled colloidalmonolayers.
Nature , 898–903 (2008).[25] Shokef, Y., Souslov, A. & Lubensky, T. C. Order by dis-order in the antiferromagnetic ising model on an elastictriangular lattice.
Proceedings of the National Academyof Sciences , 11804–11809 (2011).[26] Dennison, M., Sheinman, M., Storm, C. & MacKin-tosh, F. C. Fluctuation-stabilized marginal networks andanomalous entropic elasticity.
Physical Review Letters , 095503 (2013).[27] Ikeda, A., Berthier, L. & Biroli, G. Dynamic criticalityat the jamming transition.
Journal of Chemical Physics , 12a507 (2013).[28] Villain, J., Bidaux, R., Carton, J.-P. & Conte, R. Orderas an effect of disorder.
J. Phys. France , 1263–1272(1980).[29] Shender, E. Anitferromagnetic garnets with fluctuation-ally interacting sublattices. Sov. Phys. JETP , 178– 184 (1982).[30] Henley, C. L. Ordering by disorder: Ground state selec-tion in fcc vector antiferromagnets. Journal of AppliedPhysics , 3962–3964 (1987).[31] Henley, C. L. Ordering due to disorder in a frustratedvector antiferromagnet. Phys. Rev. Lett. , 2056–2059(1989).[32] Chubukov, A. Order from disorder in a kagom antiferro-magnet. Phys. Rev. Lett. , 832–835 (1992).[33] Reimers, J. N. & Berlinsky, A. J. Order by disorder inthe classical heisenberg kagom antiferromagnet. Phys.Rev. B , 9539–9554 (1993).[34] Bergman, Doron, Alicea, Jason, Gull, Emanuel, Trebst,Simon & Balents, Leon. Order-by-disorder and spiralspin-liquid in frustrated diamond-lattice antiferromag-nets. Nat Phys , 487–491 (2007).[35] Rubinstein, M., Leibler, L. & Bastide, J. Giant fluctu-ations of cross-linked positions in gels. Physical ReviewLetters , 405–407 (1992).[36] Barriere, B. Elatic moduli of 2d grafted tethered mem-branes. Journal De Physique I , 389–398 (1995).[37] Plischke, M. & Jo´os, B. Entropic elasticity of dilutedcentral force networks. Physical review letters , 4907(1998).[38] Tessier, F., Boal, D. H. & Discher, D. E. Networks withfourfold connectivity in two dimensions. Physical ReviewE , 011903 (2003).[39] Brazovskii, S. A. Phase-transition of an isotropic systemto an inhomogeneous state. Zh Eksp Teor Fiz , 175–185 (1975).[40] Swift, J. & Hohenberg, P. C. Hydrodynamic fluctuationsat convective instability. Physical Review A , 319–328(1977).[41] Chen, J. H. & Lubensky, T. C. Landau-ginzburg mean-field theory for nematic to smectic-c and nematic tosmectic-a phase- transitions. Physical Review A ,1202–1207 (1976).[42] Swift, J. Fluctuations near nematic-smectic-c phase-transition. Physical Review A , 2274–2277 (1976).[43] Calladine, C. R. Buckminster fuller ”tensegrity” struc-tures and clerk maxwell ruels for the construction of stiffframes. International Journal of Solids and Structures , 161–172 (1978).[44] Sun, K., Souslov, A., Mao, X. & Lubensky, T. C. Sur-face phonons, elastic response, and conformal invariancein twisted kagome lattices. Proceedings of the NationalAcademy of Sciences , 12369–12374 (2012).[45] Kane, C. & Lubensky, T. Topological boundary modesin isostatic lattices.
Nature Physics , 39–45 (2014).[46] Schmidt, M. & L¨owen, H. Freezing between two andthree dimensions. Physical Review Letters , 4552–4555(1996).[47] Schmidt, M. & L¨owen, H. Phase diagram of hard spheresconfined between two parallel plates. Physical Review E , 7228–7241 (1997).[48] Shokef, Y. & Lubensky, T. C. Stripes, zigzags, and slowdynamics in buckled hard spheres. Phys. Rev. Lett. ,048303 (2009).[49] Frenkel, D. & Smit, B.
Understanding Molecular Simu-lations (Academic Press, San Diego, 2001).[50] Binder, K. Critical properties from monte carlo coarsegraining and renormalization.
Physical Review Letters , 693–696 (1981).[51] Mao, X., Xu, N. & Lubensky, T. C. Soft modes and elasticity of nearly isostatic lattices: Randomness anddissipation. Phys. Rev. Lett. , 085504 (2010).[52] This form is similar to Eq. 4.6, except that here we ex-pand around uniformly deformed lattice of continuously varying Λ, rather than about a zigzagging state.[53] Hammonds, K. D., Dove, M. T., Giddy, A. P., Heine,V. & Winkler, B. Rigid-unit phonon modes and struc-tural phase transitions in framework silicates.
AmericanMineralogist , 1057–1079 (1996). Appendix A: Plaquette Ground State
The minimum of V PL in Eq. (4.5) in the main text can be found analytically by assuming a rhombic plaquette andsolving the set of equations given by ∂V PL /∂d = 0 and ∂V PL /∂d = 0 for τ and λ . These equations are linear, andmay be written in terms of the plaquette side b and the inner angle α , where d , = √ b √ ± sin α . The solutionshave the form τ ( b, α ) = ( b − (cid:2) b (cid:0) α + b cos α − b cos α (cid:1) − (cid:1) b (cid:0) b cos α − (cid:1) (cid:0) − b cos α + b cos α (cid:1) ,λ ( b, α ) = 3( b − b (cid:0) b cos α − (cid:1) (cid:0) − b cos α + b cos α (cid:1) . (A1)It was verified numerically that the rhombic plaquette with side length b minimizes the potential energy for the rangeof 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 H = X b V b (cid:16) | ~R b | − | ~R b | (cid:17) (B1)where b labels bonds { i, j } , V b is the interaction potential of the bond, and ~R b = ~R i − ~R j ~R b = ~R i − ~R j (B2)are the bond vectors in the reference and target states.We consider a macroscopic deformation Λ (corresponding to ǫ in the continuum theory) and fluctuations ~u ′ , sothe deformation of an arbitrary site can be written as ~R i → ~R i = Λ · ~R i + ~u ′ i . (B3)Thus for a bond, ~R b → ~R b = Λ · ~R b + ~u ′ b . (B4)The change of bond length can then be expanded for small ~u ′ b | ~R b | − | ~R b | = | Λ · ~R b | − | ~R b | + ˜ u ′ b k + (˜ u ′ b ⊥ ) | Λ · ~R b | + · · · (B5)where ˜ u ′ b k = ~u ′ b · ˆ t b (˜ u ′ b ⊥ ) = ~u ′ b · ( I − ˆ t b ˆ t b ) · ~u ′ b (B6)with ˆ t b = Λ · ~R b | Λ · ~R b | . (B7)1The terms we kept are exact to O (˜ u ).The expansion of the potential of a bond is then V b = V b Λ + ˜ u ′ b k + (˜ u ′ b ⊥ ) | Λ · ~R b | ! V ′ b Λ + ˜ u ′ b k + (˜ u b ⊥ ) | Λ · ~R b | ! V ′′ b Λ + . . . (B8)where V b Λ , V ′ b Λ , V ′′ b Λ are the potential and its derivatives at the macroscopic deformation value Λ · ~R b . (B9)We can then sort terms in powers of ˜ u ′ V b = V b Λ + V ′ b Λ ˜ u ′ b k + V ′ b Λ | Λ · ~R b | (˜ u ′ b ⊥ ) + V ′′ b Λ u ′ b k ) ! + O ((˜ u ′ ) ) . (B10)Therefore the total Hamiltonian of the lattice is H (Λ , ~u ′ ) = H (Λ) + H (2) (Λ , ~u ′ ) + O ((˜ u ′ ) ) , (B11)where H (Λ) is the energy for reference state with the uniform deformation (Λ), and H (2) (Λ , ~u ′ ) = 1 V X q ~u ′ q · D q (Λ) · ~u ′− q (B12)is the additional potential energy coming from fluctuations around the uniformly deformed reference state. Here weexpress it in momentum space, v is the area of the unit cell, and D q (Λ) = X B − cos( q · R B )] A B (B13)with A B = V ′ B Λ | Λ · ~R B | I + V ′′ B Λ − V ′ B Λ | Λ · ~R B | ! Λ · ~R B ~R B | Λ · ~R B | ! · Λ T (B14)and P B represent the sum over all bonds in one unit cell. From this we can write the dynamical matrix D as a sumof two parts D q (Λ) = v q (Λ) I + Λ · M q (Λ) Λ T . (B15)We can then calculate the expression for D q (Λ) for the special case of the square lattice model with NN harmonicsprings of spring constant k and NNN springs with the potential V NNN (∆ R ) = κ R + g
4! ∆ R . (B16) Appendix C: Lattice Free Energy
From the Hamiltonian derived in the previous section, one can calculate the free energy of a state with a uniformdeformation Λ by integrating out small fluctuations ~u ′ around this state. This lead to F (Λ) = − T ln Z D ~u ′ e − H (Λ ,~u ′ ) /T = H (Λ) + 12 T ln det [ D q (Λ)] . (C1)Using the form of D q (Λ) and the equality det( I + AB ) = det( I + BA ) we find that F only depends on the rotationallyinvariant combination of the uniform deformation ǫ ≡ (Λ T Λ − I ) / F ( ǫ ) = H ( ǫ ) + 12 T ln det (cid:2) v ( ǫ ) I + ( I + 2 ǫ ) M ( ǫ ) (cid:3) , (C2)2where g starts from O ( ǫ ) 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 ǫ . Because the dynamical matrix can be expanded as D = D + A ( ǫ ) (C3)where D = D | ǫ =0 is the dynamical matrix of the undeformed state. The free energy is then F ( ǫ ) = H ( ǫ ) + T D (cid:2) I + G · A ( ǫ ) (cid:3) , (C4)where G ≡ D − is the phonon Green’s function in the undeformed state. The expansion of F at small ǫ thus followsfrom this, F ( ǫ ) = H ( ǫ ) + 12 T ln det D + 12 T ln det (cid:2) I + G · A (cid:3) ≃ H ( ǫ ) + 12 T ln det D + 12 T Tr (cid:2) G · A − G · A · G · A + 13 G · A · G · A · G · A (cid:3) (C5)The above free energy can be calculated by performing integrals in momentum space. Because we are interested incharacterizing the square-rhombic phase transition, where τ ≡ κ/k is small, we can expand our results at small τ . Inthis limit, floppy modes of frequency √ τ lie along q x and q y axes, and to lowest order in τ , the Green’s function canbe approximated by D = (cid:18) kq x + κq y kq y + κq x (cid:19) . (C6)Integrals involving this Green’s functions can be evaluated following the calculation discussed in Ref.[40] (add thekagome lattice ref here which contain more detail). We also verified that to leading order in small τ , integrals doneusing this approximation agree with exact results. Appendix D: Identify the transition
The square and the rhombic phases can be written in terms of a uniform deformation Λ = (cid:18) s t s + w (cid:19) , (D1)where s denote a hydrostatic expansion, t measures simple shear, and w denote an anisotropic expansion in the y direction. The shear t is the order parameter of the transition: t = 0 in the square phase and t = 0 in the rhombicphase. The free energy can then be written in terms of deformations { t, s, w } .At T = 0, one can solve for { t, s, w } by minimizing H . For κ > { t = 0 , s = 0 , w = 0 } corresponding to the square lattice. For κ < τ ≡ κ/k and λ ≡ ga /k , t = ± p | τ | /λ + O ( τ / ) ,s = O ( τ ) ,w = − | τ | /λ (D2)corresponding to a shear, associated with a small vertical shrink to keep the length of the vertical bonds unchanged.Because t ∝ | τ | / this is a continuous transition at T = 0.At T > F as given in Eq. (C5). In particular, we canidentify the transition from the equations of state ∂F ( ǫ ) ∂t = 0 , ∂F ( ǫ ) ∂s = 0 , ∂F ( ǫ ) ∂w = 0 . (D3)These equations of state involves integrals that diverge at small τ , originating from floppy modes on q x and q y axes.To make the series expansion convergent, we assume, based on the T = 0 solution and the assumption λ ≫
1, that(these assumptions will be verified later for the transition) s ∼ w ∼ t ≪ τ ≪ . (D4)3In addition, we assume that the temperature is low so that in calculating the integral from T ln det D we only need tokeep terms that are singular as τ →
0, which is the leading order contribution from thermal fluctuations. Thereforeto leading order the equations of state become ∂F∂t (cid:12)(cid:12)(cid:12) s,w = kt " λ ˜ T √ τ + τλ + t ! + ˜ T √ τ + s + t w ∂F∂s (cid:12)(cid:12)(cid:12) t,w = k T √ τ + 2 s + t w ! ∂F∂w (cid:12)(cid:12)(cid:12) t,s = k ˜ T √ τ + s + t w ! , (D5)where we have defined the reduced temperature ˜ T ≡ πT ka , (D6)which is dimensionless.Combining these equations we get w = − t / w and Eq.(D5) reduce to ˜ T √ τ + s = 0 t " τλ + ˜ T √ τ + t = 0 . (D8)From these equations, it is clear that s has only one solution, but t has three solutions t = 0 corresponding to thesquare phase, and t = ± s − τλ − ˜ T √ τ , (D9)which exist for τ < − (cid:16) ˜ T λ/ (cid:17) / , corresponding to the rhombic phase, which only exist when τ is below thetransition. (When T = 0 this reduces to Eq.(D2)).As discussed in the text, we then take the self-consistency approximation, replacing τ in the denominators by thecorresponding fluctuation corrected shear rigidity r ≡ ∂ F ( ǫ ) ∂t = τ + λ ˜ T √ r + λt . (D10)The correction diverges near the transition, meaning that the square phase is always locally stable, and that thetransition between the square and the rhombic phase becomes a first order transition.Therefore the phase boundary between the square and the rhombic phases are determined by the equal free energyline. In what follows we evaluate the free energy difference ∆ F between the two phases.In the square phase, t = 0, and the fluctuation-corrected shear rigidity satisfies the equation r = τ + λ ˜ T √ r . (D11)In the rhombic phase, t is one of the nonzero solutions to Eq.(D9), and12 r + λ ˜ T √ r = − τ. (D12)4These two equations determine the fluctuation-corrected rigidity in the two phases (in the following discussion we use r for the square phase and r for the rhombic phase) for a given τ . In particular, the rhombic-phase equation (D12)only has a solution for τ < τ c = − (cid:16) ˜ T λ (cid:17) / . (D13)The free energy difference then follows from the integral, as discussed in the text,∆ F = Z t rhombic ∂F∂t (cid:12)(cid:12)(cid:12) ∂ s F =0 ,∂ w F =0 dt = Z r r ∂F∂t (cid:12)(cid:12)(cid:12) ∂ s F =0 ,∂ w F =0 ∂t∂r (cid:12)(cid:12)(cid:12) ∂ s F =0 ,∂ w F =0 dr. (D14)To evaluate this integral, we use ∂F∂t (cid:12)(cid:12)(cid:12) ∂ s F =0 ,∂ w F =0 = ktλ " ˜ T √ r + τλ + t , (D15)and Eq.(D10) (with τ → r on the denominator) to eliminate t in the above equation, and also calculate ∂t∂r (cid:12)(cid:12)(cid:12) ∂ s F =0 ,∂ w F =0 = 2 t (cid:20) λ + ˜ T r − / (cid:21) . (D16)The integration thus gives ∆ F = T " √ r r λ ˜ T − λ ˜ T r − τ √ r + 2 rτ λ ˜ T r r . (D17)By plugging in r , r we can solve for the value of τ where ∆ F = 0. To leading order, we get τ c ≃ ( − − c ) (cid:16) ˜ T λ (cid:17) / , (D18)as the boundary of the first order transition between the square and the rhombic phase, where c ≃ . τ c isslightly lower than τ c1