Quantifying the role of the lattice in metal–insulator phase transitions

Many materials exhibit phase transitions at which both the electronic properties and the crystal structure change. Some authors have argued that the change in electronic order is primary, with the lattice distortion a relatively minor side-effect, and others have argued that the lattice distortions play an essential role in the energetics of the transition. In this paper, we introduce a formalism that resolves this long-standing problem. The methodology works with any electronic structure method that produces solutions of the equation of state determining the electronic order parameter as a function of lattice distortion. We use the formalism to settle the question of the physics of the metal–insulator transitions in the rare-earth perovskite nickelates (RNiO3) and Ruddlesden–Popper calcium ruthenates (Ca2RuO4) in bulk, heterostructure, and epitaxially strained thin film forms, finding that electron-lattice coupling is key to stabilizing the insulating state in both classes of materials. For phase transitions involving concomitant structural and electronic changes, it has been challenging to determine which component provides the dominant contribution. Here, a general formalism capable of resolving this question is introduced and applied to metal–insulator transitions in two material families.

T he relative importance of electronic and lattice effects in driving phase transitions in quantum materials is the subject of long-standing interest and controversy. One issue of particular focus is the metal to insulator transition (MIT), which in many materials is accompanied by a change in atomic positions, in other words, by a lattice distortion 1 . Multiple studies have addressed various aspects of the interplay between electronic and lattice contributions to the metal-insulator transition, both theoretically and experimentally. Some works have argued that the lattice plays a crucial role 2 , while others have emphasize the electronic aspects 3 ; recently, explicit attempts at disentangling the two have appeared [4][5][6] . Information on the coupling of electronic and structural modes has been used to design and study new materials [7][8][9] and to help clarify their interconnected roles. Machine learning approaches are discovering features from comparison of MIT and non-MIT materials [10][11][12] . Nonetheless, the problem remains an open topic of exploration with a wide variety of theoretical and experimental approaches [2][3][4][5][6][7][8][9] . Yet, the relative importance of electronic and lattice effects in producing a transition has remained unclear, in part because of the lack of an appropriate theoretical framework and of computational schema for evaluating the needed quantities.
Metal insulator transition phenomena are properly analyzed by constructing an energy landscape in the space of lattice distortions and electronic order parameters as shown in Fig. 1 for NdNiO 3 . However, construction of the full energy landscape as a function of electronic order parameter and lattice configuration has not been feasible, as current quantum many-body methods do not allow the order parameters to be independently varied.
Here we show how to avoid this difficulty by building the needed energy landscapes from solutions of the equation of state (electronic order parameter as a function of lattice configuration). This information is available from modern quantum many-body methods at reasonable computational expense.
Our work was inspired by the pioneering work of Fisher and collaborators 55 , who combined elegant experimental measurements with an insightful free energy analysis to show that the observed second order nematic transition in an iron pnictide material has an intrinsically electronic origin. Our work may be viewed as a generalization of the ideas of Fisher et al. to first order transitions, where a more global view of the energy landscape is required, and as an extension of the insightful equation of state analysis of Peil, Hampel, Ederer, and Georges 6 to the computation and interpretation of the full free energy landscape. We make substantial use of concepts put forward in analyses of the metal-insulator transition in epitaxially strained Ca 2 RuO 4 5 and in nickelate heterostructures 4 . We apply the methodology developed here to two paradigmatic metal-insulator materials, the perovskite rare earth nickelates, and the Ruddlesden-Popper calcium ruthenates, settling a decades-old controversy by showing that in both cases the energetics associated with the lattice degrees of freedom are essential to the stabilization of the insulating phase.

Results
Framework. Following refs. [4][5][6] we consider a free energy F(ΔN, Q) that depends on an electronic order parameter, ΔN, and a lattice distortion, Q. The precise definitions of ΔN and Q will depend on the specifics of the system. We choose (ΔN, Q) = (0, 0) as the equilibrium configuration of the metallic phase. The insulating state then corresponds to a configuration (ΔN, Q) ≠ (0, 0) and the question of the existence of a purely electronic transition relates to the properties of F as a function of ΔN at Q = 0.
A second order transition corresponds to a linear instability of the (ΔN, Q) = (0, 0) state, in other words to a change in sign of the smallest eigenvalue of the Hessian matrix ∂ 2 F/∂ 2 ΔN, Q evaluated at ΔN = Q = 0. In general, the eigenvector associated with the negative eigenvalue will have components along both ΔN and Q, indicating that the electronic and lattice orders are coupled. A purely electronically driven transition would correspond to a change in sign of ∂ 2 F/∂ 2 ΔN, a lattice-driven transition would correspond to a change in sign of ∂ 2 F/∂ 2 Q, and in a mixed situation, neither single derivative changes sign but the lowest eigenvalue of the Hessian matrix does change sign. In the case of the nematic transition in pnictides, Fisher and co-workers were able 55 to experimentally determine ∂ 2 F(Q = 0)/∂ΔN 2 and show that it changed sign at a temperature only slightly lower than the observed nematic transition temperature, thereby establishing that a purely electronically driven nematic transition exists in these materials, with a transition temperature that is slightly enhanced by coupling to the lattice.
The transitions of interest in this paper are first order, characterized by the appearance of a new free energy minimum at which ΔN and Q are both different from zero (see e.g., Fig. 1). Study of first order transitions requires global knowledge of the free energy. One may say that the transition is electronically driven if F(ΔN, Q = 0) has an extremum at a ΔN ≠ 0 with ∂ 2 F(ΔN, 0)/∂ΔN 2 > 0 and with energy lower than F(0, 0). One may say that the transition is lattice-assisted if F(ΔN, 0) has a local minimum at a ΔN ≠ 0 but that the lattice coupling is required to make this minimum the ground state. Finally, if F(ΔN, Q = 0) has only one minimum, at ΔN = 0 but exhibits a stable minimum with ΔN ≠ 0 at a Q ≠ 0 then we may say that the transition is not electronically driven. The energy functional shown in Fig. 1 shows a transition that is not electronically driven because the global minima are at a Q, ΔN ≠ 0 but along the line Q = 0 the functional has only one minimum, at the origin.
The global knowledge of F required to assess the relative importance of electronic and lattice effects in stabilizing the insulating state has been challenging to obtain. Available quantum many-body methods can, with reasonable computational cost, obtain estimates of the optimal ΔN that minimizes F at fixed atomic positions (Q), but finding the dependence of F on ΔN at fixed Q or on Q at fixed ΔN is difficult, both because it is not straightforward to control ΔN independently of Q and because calculations of energies are very expensive.
However, the metal-insulator transitions of greatest current experimental interest share two features that greatly simplify an analysis. First, the lattice energetics is to good approximation harmonic 4-6 , as shown from calculations 4 and by the observation Fig. 1 Energy Landscape for NdNiO 3 as a function of the electronic and ionic lattice disproportionation. Energy landscape for NdNiO 3 as a function of electronic disproportionation ΔN and lattice distortion (bond disproportionation) Q as obtained from the energy functional methods introduced in this paper using data from ref. 4 . that phonon frequencies change only very slightly (~1%) across the metal-insulator transition. Second, and closely related to the first point, the coupling between the electronic order parameter and lattice modes is linear and is well determined by standard theoretical methods [4][5][6] , as are the lattice energetics. This means that we may write, schematically, 4-6 : Requiring stationarity of F with respect to variations in ΔN and Q yields the equations of state 4,6 and ∂F el ðΔNÞ Because g and K are known, Eq. 2 can be solved. Substituting the solution back into Eq. 1 produces: From this point of view, the key question is the magnitude of the lattice stabilization energy 1 8K gΔN À Á 2 : is it large enough to create a new extremum at a ΔN ≠ 0? Is it large enough to ensure that the extremum at ΔN = 0 is unstable? Answering this question requires knowledge of F el . Peil, Hampel, Ederer, and Georges 6 observed that existing many-body methods such as the dynamical mean field method enable computation of electronic properties at fixed lattice configuration Q, in effect solving Eq. 3. Here we take one further step: the ansatz for the electronic energy functional in Eq. 1 means that the computed ΔN(Q) determines dF el (ΔN)/dΔN; the derivative can then be integrated, thereby constructing F, F el itself. In practice, we perform the integration by fitting the ΔN(Q) results to a polynomial, which is analytically integrated. While this method requires a choice of polynomial form, we empirically find that the uncertainty thereby introduced is small, of the order of a few meV, comparable to the error in performing an explicit energy calculation. This approach is computationally inexpensive, relatively unaffected by noise, and the analytically specified functional forms provide additional insight, in particular enabling a straightforward examination of the global energy landscape.
The method is general: the electronic and lattice order parameters can be of scalar, vector, or tensor character, higherorder couplings can be included and any theoretical method that delivers an equation of state can be used in Eq. 1. The applicability of the model in its current form is, however, restricted to situations in which the orbital basis used to define the electronic order parameter can be clearly defined, and the order parameter energetics can be separated from the energetics of the rest of the material which can then be subsumed into the lattice stiffnesses K. We will discuss the limits of our method's applicability in more detail in the "Discussion" section, and in the Conclusion mention other materials to which it may be applicable.
We use this approach to analyze the metal-insulator transitions occurring in the rare-earth nickelates and their heterostructures, and in bulk and epitaxially strained Ca 2 RuO 4 . We find that in all of these situations the metal-insulator transition is, in the sense defined above, not electronically driven.
It is important to emphasize that the distinction between electronically driven and lattice driven transitions is not the same as the distinction between weak and strong electron correlations. A material may be said to be weakly correlated if density functional band theory (DFT) methods (or their hybrid and +U extensions) adequately describe the behavior and may be said to be strongly correlated if beyond DFT methods are required for the correct description of materials properties. While the weak/ strong correlation distinction is not the main focus of this paper we believe the available evidence indicates that beyond DFT methods are needed to properly calculate the electronic term F el of the materials we study here.
Rare-earth nickelates. The perovskite rare-earth nickelates have the chemical formula RNiO 3 ; when synthesized in bulk form the materials with R = Lu, Y, Sm, Nd, Pr exhibit a first order transition from a high-temperature metal to a low-temperature insulator while LaNiO 3 remains metallic down to the lowest measured temperature 56 . The material properties depend systematically on the choice of rare-earth ion, with LuNiO 3 exhibiting the highest metal insulator transition temperature T MIT = 600 K, SmNiO 3 having T MIT ≈ 400 K, NdNiO 3 having a T MIT ≈ 200 K and PrNiO 3 having a T MIT ≈ 130 K 17,57,58 . The materials may also be grown in heterostructure form 4,7,13,17,18,24,48,59,60 , with a small number of layers of RNiO 3 sandwiched between layers of other materials. The metal-insulator transition temperature in the heterostructures differs from that in the bulk.
Except for the La-based compound, where the symmetry is slightly different, the high temperature metallic state forms in a Pbnm structure which for present purposes may be regarded as a pseudo-cubic lattice of Ni ions with an O ion at the midpoint of each Ni-Ni bond. The insulating state is characterized by a twosublattice bond disproportionation, with one Ni sublattice having a short mean Ni-O bond length and the other by a long one, represented qualitatively in Fig. 2. The difference in mean Ni-O bond lengths defines the lattice mode Q appropriate to this transition.
In the low T insulating state, the two Ni sublattices differ in electronic configuration. This difference has been characterized in various ways in the literature 6,36,44,61-63 . We follow Peil and collaborators 4,6 and define the electronic change ΔN as the Fig. 2 Schematic of the electronic and structural disproportionation in the insulating state in rare-earth nickelates. NiO 6 octahedra disproportionate electronically and structurally in a three-dimensional checkerboard pattern in RNiO 3 metal-insulator transition compounds, with R a rare-earth. In the low-energy model of refs. 4,6 , the Ni ions corresponding to the larger octahedra have 1e − þ ΔN 2 electrons (dark blue) in their e g shell, while those corresponding to the smaller one (light blue) have 1e − À ΔN 2 . The difference in average bond lengths defines Q for the nickelate materials. Red circles correspond to oxygen atoms, while the green arrow is a measure of the oxygen displacement along the Ni-Ni direction.
difference in the occupation of the e g states obtained from a narrow window Wannier fit of the bands crossing the Fermi level, i.e., with N HF (N LF ) the total e g occupation on the octahedron with higher (lower) Ni occupation NiO 6 octahedra, while Q is defined as the average Ni-O bond length difference between the two, as shown in Fig. 2. In Fig. 3 we present results for three different bulk-form rare earth nickelates, (LuNiO 3 , SmNiO 3 , PrNiO 3 ) as a function of Q, digitized from the calculations of ΔN presented in ref. 6 . In this reference, the calculations performed are DFT+DMFT (DFT +Dynamical Mean Field Theory) calculations in the paramagnetic phase. In all compounds, two regimes are found, a low Q regime where ΔN is proportional to Q and a high Q regime where ΔN is large but weakly dependent on Q; the two regimes are separated by a discontinuity.
To determine the free energy model to which the results should be fit we observe that the two sublattice nature of the insulating state means that the free energy must be invariant under simultaneous inversion Q ↔ −Q, ΔN ↔ −ΔN, so in particular, F el (ΔN) must be even in ΔN. Motivated by the observed firstorder nature of the transition we assume that F el is a 6th order polynomial in ΔN: Solving Eq. 2 and eliminating Q from the full free energy then changes the quadratic term Note that the shift in the quadratic term involves the parameter g 2 /K, which varies slightly across the RNiO 3 series although as noted in ref. 6 the ratio g/K is almost material-independent. The K for the Lu, Sm, and Pr materials were reported in ref. 6 We fit the points in Fig. 3 to Eq. 7, using the g values presented in ref. 6 , obtaining the light solid curves shown in Fig. 3 and the fit parameters given Table 1. It is important to note that the DMFT equation of state shown in Fig. 3 has two branches, with the solution discontinuously changing from one to the other as Q is varied. In our fits we only use the data points shown as symbols in Fig. 3; these points lie outside the region of multistability (i.e., where two ΔN solutions exist within DMFT for the same Q). However, our theoretical energy function correctly reproduces behavior that was not part of the original fit, for example, the  6 ; Light lines: fits of the points to Eq. 7. Also shown is the Q(ΔN) line obtained from the Q equation of state Eq. 2. Note that the ratio of the lattice stiffness to the linear coupling K/g has the same value for all three materials. b Electronic energy, c total energy F tot = F of three RNiO 3 compounds as a function of electronic order parameter ΔN, minimized over structural order parameter (bond disproportionation) Q. d Total energy as function of lattice distortion Q, minimized over ΔN. Energies per 10 atom unit cell. Table 1 Parameters that characterize the electronic and total energy of the three RNiO 3 materials, as extracted by fitting the data from ref. 6 to Eq. 7.
The χ À1 0 ; β; γ parameters correspond to the quadratic, quartic, and sixth-order terms of the electronic energy, while g corresponds to the linear coupling between the electronic and lattice terms, and k is the lattice stiffness, as described in the main text.
locally stable large ΔN at Q = 0 solution for LuNiO 3 6 . Including points within the region of multistability does not lead to significant changes to our results. The free energy we obtain here may be thought of as the free energy that would be obtained from the two component Landau theory of ref. 6 if the metal-insulator variable were integrated out.
Panel b of Fig. 3 shows the purely electronic energy F el corresponding to the fit parameters given in Table 1. We see that for all of the materials the purely electronic theory is characterized by a locally and globally stable metallic minimum at Q = 0. Only for LuNiO 3 does the purely electronic theory even exhibit a second, insulating extremum and even in this case it is not energetically favored. Panel c presents the total free energy as a function of ΔN obtained by minimization over the lattice degrees of freedom Q. We see that inclusion of the lattice coupling is necessary to stabilize the ΔN, Q ≠ 0 solution: in other words, it is the electron-lattice coupling that drives the transition. Table 1 shows that as one moves across the RNiO 3 series from R = Lu (most insulating) to R = Pr (least insulating), the linear response electronic susceptibility of the metallic state decreases, consistent with the empirical observation that the Lu material is the strongest insulator and the Pr material is the weakest. We see also that the nonlinear terms in the energy are the same for the Lu and Sm compounds, which are the two strong insulators, but are different in the Pr compound, which is close to the metal-insulator transition. This shows that to understand composition-dependent changes one must consider more than just the variation in linear response susceptibilities. Results to be presented below for superlattices further confirm this point.
The present theory predicts that PrNiO 3 is undistorted and metallic, although close to the metal-insulator transition boundary, while the compound is empirically distorted and insulating, although again with a very low transition temperature. We believe this discrepancy arises because the transition in PrNiO 3 is from a paramagnetic metal to an antiferromagnetic insulator, in contrast to the Sm and Lu materials where the transition is between two paramagnetic phases. Inclusion of antiferromagnetism in the theory will yield a lower energy for the insulating phase of PrNiO 3 at the U, J values used here while a slightly different choice of interaction parameters would also push PrNiO 3 to the insulating side of the phase diagram.
We may take the analysis further by substituting the ΔN(Q) obtained from the solution of Eq. 7 into the full free energy expression to obtain the free energy as a function of the lattice coordinate Q shown in the panel d of Fig. 3. This is analogous to the energy that is usually calculated by electronic structure codes, which work at fixed atomic positions. The discontinuities arise from the different branches of the ΔN(Q) curves. The resulting free energy is, to a good approximation, the combination of two parabolic energy curves, one corresponding to the metallic state, and one to the insulating state. We see that in this theory the ΔN = 0, Q = 0 state is stable to lattice distortions for PrNiO 3 , marginally stable for SmNiO 3 , and unstable for LuNiO 3 .
The analysis presented here is based on a fit of the computed ΔN(Q) to a free energy model. We emphasize that this fit is in no way essential to our method; one could simply numerically integrate a suitably dense set of ΔN(Q) data. But for the procedure used here the question of the sensitivity of the results to the choice of free energy arises. To quantify the uncertainties we have refit the PrNiO 3 data shown in Fig. 3, now constraining the 4th and 6th order coefficients to have the same values as found in SmNiO 3 . Results are shown in Fig. 4. We see from panel a that this fit to the data is not quite as good, especially in the small to intermediate Q region. The fitted energy (panel b) is very similar to that shown in Fig. 3, with the value of ΔN at the minimum almost the same in the two cases and the energy of insulating solution in the constrained fit is now slightly lower than the energy of the metallic solution, placing paramagnetic PrNiO 3 slightly on the insulating side of the transition. These differences indicate that the systematic uncertainties in the method are not large.
Layered structures of NdNiO 3 /NdAlO 3 . An important modern direction in quantum materials science is heterostructuring: the ability to synthesize structures in which atomically thin planes of one material alternate with atomically thin planes of another. Previous experimental 13 and theoretical 4 work analyzed systems of the type shown schematically in Fig. 5 in which a few layers of metal-insulator transition compound NdNiO 3 were grown as a thin film sandwiched epitaxially between layers of the wide-gap insulator NdAlO 3 . The few-layer systems had significantly higher transition temperatures than did bulk NdNiO 3 . Theoretical DFT +DMFT results are available 4 for three systems: bulk NdNiO 3 and two heterostructured systems with a 4 unit cell repeat distance: one layer of NdNiO 3 followed by three of NdAlO 3 (monolayer), and two layers of NdNiO 3 followed by two of NdAlO 3 (bilayer). The ΔN(Q) results are shown in Fig. 6 Two competing effects were found 4 : the electronic confinement of the electrons in the NdNiO 3 by the nearby NdAlO 3 layers favors a higher electronic disproportionation ΔN, while the energy cost of imposing a lattice distortion on the adjacent Al-O octahedra is equivalent to an increase in the stiffness K of the NiO 6 lattice with a K = 31.6914eV/Å 2 , for the bulk material, K = 35.6269 eV/Å 2 for the bilayer and K = 40.4753 eV/Å 2 for the monolayer; all energies per 10 atom unit cell. The values for g are those quoted in ref. 4 ; we reproduce them here for convenience: g = 3.2212 eV/Å for the bulk, g = 3.3465 eV/Å for the bilayer, and g = 3.3780 eV/Å for the monolayer. Performing an analysis similar to that presented in section "Rare-earth nickelates" leads to the energy curves presented in Fig. 6 and to the fit parameters shown in Table 2. The original data with the fit are presented for convenience in Fig. 6 Turning first to the purely electronic correlation energy shown in the (b) panel of Fig. 6, we observe that neither NdNiO 3 nor the heterostructured materials show a purely electronically driven insulating state; however, hints of an inflection point, the precursor of the formation of a higher ΔN extremum, can be seen in bilayer and monolayer curves around ΔN = 1 reflecting the effect of quantum confinement in increasing electron correlation effects. Inclusion of the lattice energy produces a global, insulating minimum in all three materials. The monolayer and bilayer show a significantly larger energy difference between the insulating and the metallic states than does bulk NdNiO 3 , in agreement with their higher T MIT .
The origin of the differences between monolayer, bilayer, and bulk nickelates is surprising: the modest decrease in the magnitude of g 2 /K as we pass from bulk to bilayer to monolayer reflects the increase in lattice stiffness, opposing the order, and nearly counterbalances the modest decrease in χ À1 0 reflecting the increased correlation physics of the confined metallic state. The more important effect, however, is the change in the magnitude of the higher order (ΔN 4 , ΔN 6 ) terms shown in Table 2. The change can be seen in Fig. 6 from the variation of the critical Q below which the insulating state is not stable. In the calculations presented in the previous subsection the higher order terms were also found to be different in PrNiO 3 (proximal to the metal-insulator transition) than in LuNiO 3 and SmNiO 3 (farther from the transition point).
Metal-insulator transition in bulk Ca 2 RuO 4 . Ca 2 RuO 4 forms in a slightly distorted version of the n = 1 Ruddlesden-Popper structure, a layered structure consisting of RuO 2 planes separated by pairs of CaO layers. Each Ru is six-fold coordinated by oxygen, forming RuO 6 octahedra. The relevant electronic states (actually Ru-O antibonding states) may be thought of as Ru t 2g symmetry d-states (d xy , d xz , d yz ) with 4 t 2g electrons per Ru. Below 340 K the material undergoes a first-order transition from a high temperature metallic to a low temperature insulating phase. No symmetry is broken at the transition: the two phases share the same crystal symmetry but differ in the occupancies of the d orbitals and shape of the RuO 6 octahedra, and the relative occupancies of the d-levels. As sketched in Fig. 7, the metallic state is characterized by an approximately equal occupancy of the t 2g orbitals, while in the insulating state the xy orbitals are doubly occupied and the xz/yz orbitals each contain a single electron. The appropriate electronic order parameter is the occupancy difference between the xy and xz/yz orbitals No crystal symmetry protects the orbital occupancies in the high temperature metallic state, which is characterized by a ΔN = N H close to, but not exactly, zero. The change in electron occupancy across the metal-insulator transition occurs simultaneously with a decrease in the Ru-O apical bond length and an increase in the in-plane Ru-O bond length (a flattening of the octahedra). At the transition, there is a change in the lattice constants, and rearrangements of other internal coordinates. The harmonic nature of the lattice Hamiltonian means that most of the lattice modes may be integrated out, leaving an effective theory involving two structural degrees of freedom. The two lattice degrees of freedom are needed, because if the point symmetry of the Ru ion is lower than cubic, then both the unit cell volume change and the relative octahedral distortion couple linearly to the differential level occupancy ΔN.
Following ref. 5 we write the two relevant structural degrees of freedom as a variable Q 3 parametrizing a volume preserving change in the c-direction Ru-O bonds relative to the average inplane Ru-O bonds, and a change Q 0 in the octahedral volume and define the lattice variables such that in the high-T metallic state Q 3 = Q 0 = 0. In terms of changes δx, δy, δz to the three Ru-O bond lengths, defined as: We now build the free energy. We define the lattice distortion relative to the high-T metallic state as the vector Q ! ¼ ðQ 3 ; Q 0 Þ; the lattice restoring term K is a 2 × 2 matrix with entries determined in previous work 5 to be K 33 = 17.7, K 03 = 7.6, K 00 = 46.2 eV/Å 2 per formula unit. We write the linear combination of the Q that couples to the electronic disproportionation as F ! Á Q ! with F ! ¼ F 3 ð1; ÀλÞ. Reference 5 finds F 3 = 2.8 eV/Å and λ = 0.45.
Here ΔN H is the value of ΔN in the high temperature insulating phase; it is nearly, but not quite, zero.  For the purely electronic energy we chose the 4th order form: Odd powers occur because no symmetry is broken at the transition. Figure 8 shows previously published results 5 for ΔN as a function of the relevant combination of lattice distortions along with our fit to Eq. 11. As in the previous section we have not used points in the coexistence region for the fit. We find σ 0 = −0.009, a = 0.70, b = −1.718, c = 1.28.
The higher curve (blue) in Fig. 8, panel b, shows the purely electronic free energy resulting from the fit. We see that in the absence of lattice effects the metallic, undistorted state is energetically favored and that there is not even a metastable minimum corresponding to the insulating state, although there is a hint of an inflection point that is a precursor of a higher ΔN extremum. The lower trace (dark orange line) shows F, the full free energy after the lattice modes have been integrated out. We see that the inclusion of the lattice energetics qualitatively changes the free energy, strongly favoring large ΔN. Indeed we see that the energy does not have a minimum in the physically allowed range ΔN ≤ 1; rather it is minimized at the boundary ΔN = 1, as |ΔN ≤ 1| is the largest allowed value. We again conclude that the The χ À1 0 ; β; γ parameters correspond to the quadratic, quartic, and sixth-order terms of the electronic energy, while g corresponds to the linear coupling between the electronic and lattice terms, and k is the lattice stiffness, as described in the main text.  transition in this compound is driven by the lattice contribution to the energetics.
Epitaxially constrained Ca 2 RuO 4 . Ca 2 RuO 4 has also been grown expitaxially on insulating substrates. The epitaxial constraint means that the in-plane lattice parameters of Ca 2 RuO 4 are forced to be equal to the in-plane lattice parameters of the substrate material. A priori, this constraint does not fix the values of the internal coordinates including the the Ru-O in-plane bond length of relevance here. Density functional theory calculations, however 5 , indicate that in practice the material adapts to the epitaxial constraint by adjusting the in-plane Ru-O bond lengths, rather than by rotating the Ru-O 6 octahedra, so that the percentage change in the in-plane Ru-O bond lengths is fixed by the percentage change in the in-plane lattice parameters. The zdirection lattice constants and Ru-O apical bond lengths are still of course free to relax. This means that the theory of the previous section can be simply adapted to the epitaxial case by a change of variables. We write: and with δx 0 and δy 0 fixed by the epitaxial constraint. It is also convenient to define the zero of δz to be the value that δz would take at ΔN = N H . We obtain: with and The epitaxial constraint thus has three effects: it increases the net stiffness to lattice distortions by preventing relaxations of the lattice coordinates associated with in-plane bond lengths, making the insulating phase more difficult to obtain. Second, it reduces the electron lattice coupling (Eq. 17) which also makes it more difficult to obtain an insulating phase. Finally, it provides a term linearly proportional to ΔN, which for positive epitaxial strain favors the insulating phase. Panel a of Fig. 9 shows the change in apical Ru-O bond length as a function of epitaxial strain; panel b the occupancy difference ΔN. The first order transition is evident. Note that for an epitaxial strain matching the in-plane lattice constant of the high temperature metallic phase, insufficient lattice energy is available to stabilize the insulating phase. As the tensile strain is increased, the octahedral deformation, which couples linearly to ΔN, increases, and beyond a critical value drives a first order transition.
Experiments (refs. 64 Figure 8 shows that the theory, with no adjustable parameters, predicts that for the theoretically predicted coupling F 3 = 2.8 eV/Å the films grown on NdAlO 3 and LaAlO 3 should be metallic, while the films grown on NSAT and NdGaO 3 should be insulating. In the experiment, the LaAlO 3 -strained material has a transition at T ≈ 200 K from a high temperature metal to a low temperature weak insulator/bad metal, while the others are metallic and insulating as predicted by the theory. The qualitative agreement is good; the quantitative discrepancies may arise from a more subtle relation between Ru-O bond length and in plane lattice constant than was assumed in ref. 5 or from small systematic errors in the theory. The good correspondence between experiment and theory shows the utility of the methodology introduced here. Temperature Dependence. The methods presented here provide some insight into the physics of the temperature-driven metal-insulator transition, indicating an important area for future research. We see from Figs. 3 and 6 that for the nickelates the low T energy landscape is characterized by two minima, with energy differences per formula unit ranging from 80 meV (LuNiO 3 ) to 25 meV (NdNiO 3 ). In the nickelates, as the temperature increases, a first order metal-insulator transition occurs. The lattice distortion Q exhibits negligible temperature dependence over the entire insulating phase, indicating that the minima remain robust and the transition is driven by an entropic effect that lifts up the free energy of the large ΔN extremum relative to that of the ΔN = 0 extremum.
To investigate whether the entropic effect arises from the local correlations captured by the dynamical mean field approximation, for bulk NdNiO 3 we have extended the calculations shown in Fig. 6 to much higher temperatures by performing calculations with the same methodology and interaction parameters as in 4 but changing the electronic temperature. The parameters, as obtained from the fit for certain temperatures can be found in Table 3. Figure 10 shows the results: all of the parameters in the electronic free energy vary systematically with temperature, decreasing in magnitude as the temperature is increased from~290 to 1000 K; the result is a very weakly first order transition at ≈1080 K (the model actually passes very close to the tricritical point at which the second and fourth order coefficients of the theory simultaneously change sign). This behavior, while theoretically interesting, is inconsistent with the observed strongly first order character and much lower transition temperature behavior, suggesting that either the single site dynamical mean field approximation does not adequately describe the temperature dependence of the electronic contribution to the free energy, or that an entropic effect associated with the lattice degrees of freedom plays an important role. Indeed the neglect of spatial fluctuations generically means that single-site DMFT approximations overestimate transition temperatures 66,67 . It should also be noted that in this material antiferromagnetism plays an important role in the ordering, and that a direct computation 44 of the free energy F(ΔN, Q) using a different DMFT formalism places paramagnetic NdNiO 3 on the paramagnetic metal side of the transition and yields a stronger low T temperature dependence, and shows for LuNiO 3 a temperature dependence of F(ΔN, Q) in the range T < 500K very similar to that found here for NdNiO 3 . Ca 2 RuO 4 presents different issues. A substantial temperature dependence is observed 68 across the insulating phase. Since the electronic order parameter ΔN in this material is fixed at its largest possible value for a range of couplings, the temperature dependence must imply a temperature dependence of the electron-lattice coupling F 3 as proposed for related reasons in 5 . Within single-site dynamical mean field theory, there is no theoretical justification for a strongly temperature-dependent coupling parameter, but on the assumption that changes in electron-lattice coupling F 3 are a proxy for changes in temperature we show in Fig. 8, panel d, the energy computed for different values of F 3 . We see that a first order transition occurs at an F 3 ≈ 1.58 ≈ 0.6 of the F 3 we estimated from the low temperature calculation. The phase diagram presented in Fig. 8 panels e, and c, then indicates that films grown on substrates with a tensile epitaxial constraint exhibit significantly higher transition temperatures than found in bulk, roughly consistent with experiment.

Discussion
It is important to critically examine the assumptions that underly our results. The first assumption enabling the construction is that the lattice energy is quadratic in deviations from an equilibrium position at fixed value of the electron order parameter, in other words that the physics controlling the relevant atomic force constants arises from chemical bonds inside the solid that are only weakly affected by the transition from metal to insulator. This and the harmonic nature of the lattice response are found in Table 3 Parameters that characterize the temperature dependence of the electronic energy as extracted by fitting the data in Supplementary Note 1; electronic energies presented in Fig. 10. The χ À1 0 ; β; γ parameters correspond to the quadratic, quartic, and sixth-order terms of the electronic energy, as described in the main text. Fig. 10 Evolution of the electronic and total energy landscape of NdNiO 3 as a function of electronic temperature. a Electronic free energy as function of electronic disproportionation ΔN computed for NdNiO 3 at temperatures from 290 to 1160 K. b Total free energy computed for bulk NdNiO 3 for the same temperature range. c Total expanded view of total free energy for bulk NdNiO 3 at temperatures close to the transition. d Optimal bond disproportionation Q, as obtained by a linear transformation from the electronic disproportionation ΔN, using the lattice stiffness K and electron-lattice coupling term g: Q(T) = ΔNðTÞ g 2K . Data for fits provided in Supplementary Note 1.
density functional calculations (see e.g., the supplementary material of ref. 4 ) and are further supported by the observation that phonon frequencies do not change much across the transitions of interest. The second assumption is that the coupling between electronic and lattice degrees of freedom is linear, and may be determined a priori. This second assumption has been confirmed by explicit calculations (see, for example, the supplementary information in ref. 4 ). Further investigation of these assumptions, and identification of compounds in which they break down, is an important topic for future research. One important direction is a detailed comparison of results obtained along the lines indicated here to direct computations of energies and free energies 44,62,63,69 .
Given these assumptions, a many-body calculation of electronic configuration as a function of lattice distortion may be used to construct free energies, essentially by integrating the equation of state. We accomplish the integration by fitting the equation of state results to a polynomial which is integrated analytically. This procedure is convenient because it avoids ambiguities associated with the coexistence region of the first order transition, but is not necessary.
We applied the methodology to two currently interesting families of compounds, that exhibit first order metal-insulator transitions as the temperature is decreased: the rare earth perovskite nickelates RNiO 3 , and Ca 2 RuO 4 . The two materials were chosen because they exhibit qualitatively different electronlattice couplings. The transition in the RNiO 3 is a symmetrybreaking transition leading to two inequivalent Ni sites and an alternating pattern of long and short Ni-O bond lengths. The transition in Ca 2 RuO 4 is isosymmetric, with no broken crystal symmetries between the two states. In both materials, we find that if the lattice coordinate is set to the value appropriate to the high temperature metallic state, the resulting energy F(ΔN, Q = 0) has only one minimum, at ΔN = 0. It is only after the coupling to the lattice is included that the total functional F(ΔN, Q) acquires a minimum at a ΔN, Q ≠ 0. We, therefore, conclude that lattice effects are essential in driving the metal-insulator transition in both compounds, rather than being merely a minor consequence of a fundamentally electronically driven transition, thereby settling a long-standing controversy. It is important to note however that although the purely electronic theory does not produce a metal-insulator transition, the correlation contribution to the electronic energy is still important. DFT-level calculations would lead to a F el , which rises so steeply with ΔN that the only extremum would be at ΔN = Q = 0. A beyond-DFT theory is needed to obtain an electronic energy which, although it does not by itself have a minimum at a non-zero ΔN, is "soft" enough to enable the lattice energies to produce the ordered state.
It is important to consider the limitations of our conclusions. First, while the approach introduced here will work with any method that calculates the electronic order parameter as a function of lattice distortion ΔN(Q), in practice the information available in the literature comes from the density functional plus dynamical mean field methodology. While this method successfully produces results that are consistent with many experiments, its precise quantitative accuracy is unknown. The method does require an identification of correlated orbitals, it relies on the assumption that density functional theory gives an adequate account of the energetics of the uncorrelated orbitals, it requires specification of interaction parameters, and its solution of the many body problem involves a stringent locality assumption. Cross checking via other methods (as for example was done for Ca 2 RuO 4 in ref. 70 ) would be desirable. Further, there are several variants of the DMFT plus DFT methodology, differing in the choice of energy window and the representation of correlated states 71 ; the dependence of our conclusions on the implementation is an interesting question.
The specification of interaction parameters is of particular importance. As shown e.g., in refs. 3,4 , as the interaction parameters are increased in magnitude, a transition can be generated in the purely electronic theory. The precise statement made here is that with interaction parameters determined by reproducing physical properties including gaps, and effective masses, and structural distortion, the purely electronic theory does not exhibit a phase transition: electron-lattice coupling plays an essential role in stabilizing the insulating states.
It will also be noted that while the theory presented here reproduces trends and orders of magnitude very well, it has some quantitative deficiencies, for example predicting that PrNiO 3 is metallic but close to the phase boundary when in fact it is insulating but close to the phase boundary, with the lowest transition temperature of any bulk member of the RNiO 3 nickelate family. This deficiency is likely to be remedied by the inclusion of magnetism in the theory, since in PrNiO 3 and NdNiO 3 , but not in the other compounds studied in this paper, the metal-insulator transition is accompanied by a magnetic transition. Similarly, the phase diagram of epitaxially constrained Ca 2 RuO 4 films is qualitatively correct but exhibits quantitative discrepancies with experiment. We note that in this paper no attempt was made to adjust the parameters from the literature or to fine-tune the fits to obtain better agreement with the experiment. Small changes, reflecting moderate parameter uncertainties and modest systematic errors in the DFT +DMFT approach, would likely cure these discrepancies.

Conclusions
Electronic phase transitions in quantum materials are almost always accompanied by some form of lattice distortion, and in many specific cases, the question of whether the transition is driven by the lattice or by correlated electron energetics has been hotly debated. In theoretical terms, the answer is determined by the properties of a free energy functional F(ΔN, Q) that depends on suitably defined electronic (ΔN) and lattice (Q) modes. For second order transitions, the needed information can often be determined from symmetry arguments along with the values of a small number of parameters determined from the linear response of the symmetry-unbroken state. However, many transitions of current interest, in particular metal-insulator transitions, are first order, so that knowledge of the free energy landscape over a wide range is required. This information has been more difficult to obtain both because of the simple expense of computing an entire energy landscape and because these computations require fixing the electron order parameter at values that do not extremize the energy. In this paper, we proposed a method that constructs the needed free energy from the equation of state (electron order parameter as a function of lattice distortion). The method is computationally feasible and can be applied using any many-body method that obtains the equation of state. We demonstrate the power of the method by using it to resolve the long-standing question of the role of the lattice in the metal-insulator transitions in two representative classes of quantum materials, the perovskite rare earth nickelates, and the Ruddlesden-Popper calcium ruthenates. We find that in both materials the lattice degree of freedom is essential to the metal-insulator transition.
We also observe that while we presented results based on dynamical mean field calculations, the method is generic in nature and can easily be used with other electronic structure methods, including DFT, DFT+U 72,73 , Gutzwiller 74,75 , Auxiliary-Boson 76-79 -particularly with the advent of new codes that calculate energy and naturally handle symmetry breaking [80][81][82][83][84] , and quantum Monte Carlo 85 .
The methods can be applied to quantify the importance of lattice effects in other systems of current interest, including CaFeO 3 86,87 which exhibits charge ordering and lattice disproportionation, perovskite titanates and vanadates, where lattice distortions couple to orbital ordering [88][89][90] , VO 2 where a dimerization instability occurs 91 , V 2 O 3 where the metal-insulator transition couples to the volume of the material 92 , and the rare earth manganites where electronic charge, orbital and magnetic ordering are tightly coupled to complex lattice distortions [93][94][95] , two-dimensional correlated van der Waals dihalides and trihalides MX 2 , MX 3 , with M a transition metal ion and X a halogen [96][97][98][99] , as well as in charge density waves in transition metal dichalcogenides 100,101 .
The access to the free energy provided by our methods enables additional insights. Issues relating to temperature dependence were discussed in the temperature dependence section "Temperature Dependence". The energy landscape shown in Fig. 1 makes it clear that only one narrowly defined path connects the metallic and insulating minima; this information may be used to investigate the kinetics of order parameter nucleation if the material is supercooled or superheated across the phase boundary. Further, the differing roles of the electronic and lattice degrees of freedom in defining the energy landscape make possible an analysis of the kinetics of state evolution after excitation. For example, excitation might rapidly heat the electrons (leading to changes in F el ) but only slowly heat the lattice degrees of freedom (which would in any event respond more slowly). Extension of the theory presented here to include the antiferromagnetism that occurs in the rare earth nickelates may help understand recent ultrafast experiments exploring the relation between antiferromagnetism, lattice distortions, and metallicity in NdNiO 3 102,103 .

Methods
The new numerical ΔN(Q) data for NdNiO 3 used to build the electronic energies in this paper is obtained by using identical methodology and parameters as in our previous work (ref. 4 ). We've used Quantum Espresso 104 for density functional theory calculations, Wannier90 105 to build the low-energy Wannier model, and dynamical mean field theory calculations using the TRIQS implementation 106 , using the continuous time hybridization 107 solver. We've used the exact same numerical and interaction parameters as in (ref. 4 ), but with a different electronic temperature which can be obtained by changing the parameter β.

(Ca 2 RuO 4 ).
Data obtained either from new calculations or digitizing from previous work, is fit using Matlab and built into our equation of state (3), as appropriate for each material.

Data availability
The data used for this work and scripts used to fit the data can be found in the manuscripts referenced herein, both original and digitized from previous work, 4-6 , can be accessed directly on GitHub at: https://github.com/alexandrub53/EnergyLandscapes.

Code availability
The scripts used to fit the DFT+DMFT data and obtain electronic and total energies throughout the paper can be accessed directly on GitHub at: https://github.com/ alexandrub53/EnergyLandscapes.