Unraveling the Mott-Peierls intrigue in Vanadium dioxide

Vanadium dioxide is one of the most studied strongly correlated materials. Nonetheless, the intertwining between electronic correlation and lattice effects has precluded a comprehensive description of the rutile metal to monoclinic insulator transition, in turn triggering a longstanding “the chicken or the egg” debate about which comes first, the Mott localisation or the Peierls distortion. Here, we show that this problem is in fact ill-posed: the electronic correlations and the lattice vibrations conspire to stabilise the monoclinic insulator, and so they must be treated on equal footing not to miss relevant pieces of the VO2 physics. Specifically, we design a minimal model for VO2 that includes all the important physical ingredients: the electronic correlations, the multi-orbital character, and the two components antiferrodistortive mode that condenses in the monoclinic insulator. We solve this model by dynamical mean-field theory within the adiabatic Born-Oppenheimer approximation. Consistently with the first-order character of the metal-insulator transition, the Born-Oppenheimer potential has a rich landscape, with minima corresponding to the undistorted phase and to the four equivalent distorted ones, and which translates into an equally rich thermodynamics that we uncover by the Monte Carlo method. Remarkably, we find that a distorted metal phase intrudes between the low-temperature distorted insulator and high-temperature undistorted metal, which sheds new light on the debated experimental evidence of a monoclinic metallic phase.

The crystal structure of rutile VO 2 is formed by equally spaced apart Vanadium atoms sitting at the centre of edge-sharing oxygen octahedra that form linear chains along the R c-axis, which we shall denote as c R , see Fig. 1. The tetragonal crystal field splits the 3d-manifold into two higher e g and three lower t 2g levels. In the oxidation state V 4+ , the single valence electron of Vanadium can, therefore, occupy any of the three t 2g orbitals, which are in turn distinguished into a singlet a 1g (or d || ) and a doublet e π g (or d π * ), having, respectively, bonding and nonbonding character along the c R -axis. The M1 phase is instead characterised by an anti-ferroelectric displacement of each Vanadium away from the centre of the octahedra, see Fig. 1, so that the above-mentioned chains, from being straight in the R phase, turn zigzag and dimerise [36,37].
A simple portrait of the transition in VO 2 was proposed in 1971 by Goodenough [38]. According to his proposal, the basal-plane component of the anti-ferroelectric distortion raises the energy of e π g with respect to the a 1g [39]. In addition, the c R component of the distortion, which drives the chain dimerisation, opens a hybridisation gap between bonding and anti-bonding combinations of the a 1g . For large enough crystal field splitting and hybridisation gap, the bonding combination of the a 1g fills completely, while the anti-bonding as well as the e π g get empty, hence the insulating behaviour. The Goodenough's mechanism for the metal-insulator transition in VO 2 thus relies on a single-particle description: the Peierls instability of the quasi-one-dimensional a 1g band that becomes half-filled after the growth of the crystal field drained the e π g orbital.
However, Zylbersztejn and Mott [40] soon after argued that the role of electronic correlation cannot be neglected as in the Goodenough's scenario. Indeed, a tiny ∼ 0.2% substitution of V with Cr changes the lowtemperature insulator from the M1 crystal structure to a new monoclinic phase, named M2, where dimerised and zigzag chains alternate [34,41]. The M2 phase can be also stabilised under hydrostatic pressure or uniaxial stress [31,32,[42][43][44][45]. In addition, a triclinic (T) phase with intermediate structural properties [41] was shown to intrude between M1 and M2. The zigzag undimerised chains in M2 are still insulating and display magnetic properties akin those of a spin-1/2 antiferromagnetic Heisenberg chain [41,43,46]. This likeness can be rationalised only invoking sizeable electronic correlations. Given the low concentration of substitutional Chromium or the small value of uniaxial stress required to stabilise arXiv:1906.10632v1 [cond-mat.str-el] 25 Jun 2019 M2, it is reasonable to conclude that M1 must be as correlated as M2 [47][48][49][50].
We believe that, even though electronic correlations are likely necessary, they are nonetheless not sufficient to explain the phase diagram of VO 2 . It is known that a strong enough repulsion may drive a Mott transition in a three-band Hubbard model at the density of one electron per site [51]. Therefore, it is well possible that the insulating phase of VO 2 is driven by correlations alone, and that the structural distortion below T c is just the best way the Mott insulator can freeze the residual spin and orbital degrees of freedom to get rid of their entropy. However, should that be the case, VO 2 would most likely remain insulating even above T c , which is not the case, all the more so because k B T c is more than one order of magnitude smaller than the optical gap in the M1 phase [52]. For the same reason, we must exclude a transition merely driven by the larger electronic entropy of the metal.
We are thus inclined to believe that the structural distortion is also necessary to stabilise the insulating phase in VO 2 , but, once again, not sufficient in view of the behaviour of the M2 phase, and of the bad metal character of the R phase [53][54][55]. It is therefore quite likely that Goodenough's scenario is after all correct, though it requires an active contribution from electronic correlations.
Indeed, different DFT-based calculations, which should properly account for the effects of the lattice dis-tortion on the electronic structure, though within an independent-particle scheme, do not agree one with another, and none explains at once all experiments. For instance, straight LDA or GGA methods do not find any gap opening in M1 and M2 phases [56,57]. Such gap is instead recovered by GW [58][59][60] or LDA+U [61][62][63], in all its variants. However, GW does not give easy access to the total energy, and therefore it does not explain why low temperatures should favour the M1 distorted phase against the rutile undistorted one. In turns, LDA+U or GGA+U calculations, known to overemphasise the exchange splitting, predict the existence of local moments even in the rutile phase [61][62][63], not observed in experiments [64]. Relatively recent calculations based on HSE hybrid functionals bring even worst results: both rutile and M1 phases are predicted to be magnetically ordered insulators, with the former lower in energy [65,66], even though earlier calculations were claimed to be more in accordance with experiments [25]. In turn, mBJ exchange potentials seem to predict the proper conducting behaviour of the R and M1 phases, as well as their lack of magnetism [67], which is erroneously predicted to occur also in the M2 phase [63]. This suggests that suppression of magnetic moments is somehow the rule of mBJ functionals applied to VO 2 , which only by chance is the correct result for R and M1 phases. Finally, calculations based on PBE0 hybrid functionals properly account for the magnetic and electronic properties of M1 and M2 phases, but predict ferromagnetism in the rutile structure, at odds with experiments [33], as well as the existence of a never observed ferromagnetic and insulating monoclinic phase, dubbed M0 [68], also predicted by PBEsol functionals [69].
One might expect that combining ab-initio techniques with many-body tools, e.g., DFT with dynamical meanfield theory (DMFT) [70], should work better and finally provide uncontroversial results in accordance with experiments. Unfortunately, different calculations by state-ofthe-art DFT+DMFT methods do not even agree about an unanimous view of the M1 monoclinic phase. Specifically, M1 has been regarded from time to time as a correlation-assisted Peierls insulator [24,71], or, vice versa, as a Peierls-assisted Mott insulator [72], or, finally, as a genuine Mott insulator [26,73,74].
In view of the above controversial results, we think it is worth desisting from describing VO 2 straight from first principles, and rather focusing on a minimal model, which can include all the ingredients that are, by now, widely accepted to be essential. As we mentioned, electron-electron correlations must play an important role and thus need to be included and handled in a truly many-body scheme. At the meantime, the coupling of the electrons to the lattice is equally important and must be included as well. We earlier mentioned that the monoclinic distortion in the M1 phase actually entails two different antiferrodistortive components: the basal-plane displacement of V from the octahedron centre, resulting into a zigzag shape of the formerly straight chains, and the out-of-plane displacement that produces the chain dimerisation. The two phenomena may actually occur separately, as indeed proposed by Goodenough [38], who argued that, generically, the basal-plane distortion should appear at higher temperatures than dimerisation. Indeed, time-resolved spectroscopy measurements during a photoinduced monoclinic-to-rutile transition have shown that dimerisation melts on earlier time-scales than the basal-plane displacement [37,75,76], which therefore must be distinct and actually more robust than the former. We must mention, however, that this conclusion does not agree with other experiments [77][78][79][80]. More convincing evidence is offered by the monoclinic metal that intrudes, under equilibrium conditions, between rutile metal and monoclinic insulator at ambient pressure [81][82][83][84], and nor just above a critical pressure as originally believed [85]. This phase might correspond to a crystal structure where dimerisation is almost melted unlike the zigzag distortion [69,83], so that e π g are still above the a 1g , though the dimerisation is too weak to stabilise at that temperature an hybridisation gap within the a 1g band [27]. Even the disappearance prior to the metalinsulator transition [86] of the so-called singlet peak, which is associated to dimerisation and observed in optics, can be regarded as a consequence of the melting of dimerisation preceding the complete monoclinic-torutile transformation. All the above experimental facts point to the need to treat separately the basal-plane displacement and the out-of-plane one. Finally, the importance of the basal plane antiferrodistortive mode suggests the last ingredient to be considered: the multi-orbital physics. This aspect was originally emphasised by Goodenough [38] and successively confirmed by many optical measurements [52,87,88].
To summarise, we shall consider a microscopic model which includes and treats on equal footing the following relevant features: 1. the electron-electron correlations and the coupling to the lattice distortion [46,53,; 2. the existence of two different antiferrodistortive components, each playing its own distinctive role [37,38,75]; 3. the multi-orbital physics [38,52,87,88].
with the minimal requirement of capturing, at least at a qualitative level, the following aspects of the VO 2 physics: A. the existence of an undistorted paramagnetic metal and a monoclinic distorted insulator [42,[110][111][112][113]; B. the first-order character of the transition between them [18,[114][115][116][117][118][119][120][121][122][123][124]; C. the possible existence of an intermediate monoclinic metal [81][82][83][125][126][127][128][129][130][131][132][133]; Many models have been already put forth to describe VO 2 . However, most of them focus either on the role of the electron-electron correlations, or on that of the electron-lattice coupling [27,29,[134][135][136][137][138][139][140][141][142][143][144][145][146][147], and thus do not allow accessing in a single framework the whole VO 2 phase diagram, e.g., the points A., B. and C. above. There are actually some exceptions where electron-electron and electron-lattice interactions have been considered on equal footing [148][149][150]. In particular, the model studied in [149] includes explicitly all ingredients listed above, for instance, the two distinct effects of the monoclinic distortion parametrised, however, just by a single phonon mode. Nonetheless, the mean-field treatment of the electron-electron interaction, despite its strength being comparable to the conduction bandwidth, yields not surprisingly to the formation of local moments in the rutile metal, not in accordance with magnetic measurements [64]. This negative result, highlighted by the same authors of Ref. [149], solicits for a more rigorous treatment of the interaction. This is actually the scope of the present work, which is organised as follows. In Sec. II we introduce a simple model that includes the three ingredients previously outlined, which we believe should capture the main physics of Vanadium dioxide. In Sec. III we discuss the dynamical mean-field theory (DMFT) approach to the model Hamiltonian, and presents in Sec. III A its ground state phase diagram. In Sec. IV we discuss the insulator-metal transition that occurs in our model upon raising the temperature. In Sec. IV A we discuss the case in which such transition is driven solely by the electronic entropy, hence neglecting the lattice contribution to entropy, whereas in Sec. IV B the opposite case. We will show that the latter situation is rather suggestive, since it foresees different transition temperatures of the two antiferrodistortive components, as predicted by Goodenough [38]. In turn, this result might explain the evidence supporting the existence of a monoclinic metal phase. Finally, Sec. V is devoted to concluding remarks.

II. THE MODEL
As we mentioned, the orbitals that are relevant to describe the physics of VO 2 are the Vanadium 3d − t 2g ones, comprising the a 1g singlet and e π g doublet, which host a single conduction electron. We believe that in this circumstance the doublet nature of the e π g is not truly essential; what really matters is the distinction between a 1g and e π g based on their bonding character with the ligands and response to atomic displacement. Therefore, in order to simplify our modelling without spoiling the important physics, we shall associate the e π g doublet with just a single orbital [135,151], which, together with the other orbital mimicking the a 1g singlet, give rise to two bands, band 1 ↔ a 1g and band 2 ↔ e π g , which accommodate one electron per site, i.e., they are quarter filled.
The other ingredient that is necessary to properly de-scribe VO 2 is the electron-electron Coulomb interaction. However, since the main role that Coulomb repulsion is believed to play is to suppress charge fluctuations on V 4+ , we shall ignore the long range tail and replace Coulomb repulsion with a short-range interaction. Finally, we need to include the coupling to the lattice. For simplicity, we shall focus our attention only on the rutile and monoclinic M1 phases, as such ignoring the M2 phase, which is actually regarded by some as just a metastable modification of the M1 structure [43,44,145]. Under this assumption, we can model the lattice antiferrodistortion through a two-component zone boundary mode at momentum Q, with displacement X = (X 1 , X 2 ) and classical potential energy Φ X 1 , X 2 . The X 1 and X 2 components model, respectively, the dimerising outof-plane displacement and the band-splitting basal-plane one, see Fig. 1 [134,152].
The model Hamiltonian is thus written as the sum of three terms: H el is the purely electronic component reading: where n a,k is the occupation number at momentum k of the band a = 1, 2, n i the electron number operator at site i, µ the chemical potential used to enforce the quarter filling condition and, finally, U is the on-site Hubbard repulsion.
With the aim to reduce the number of independent Hamiltonian parameters, we assume that the densityof-states (DOS) D 1 ( ) and D 2 ( ), of the band 1 and 2, respectively, have same bandwidth and centre of gravity, which we shall take as the zero of energy. In addition, we consider both DOS symmetric with respect to their centre, and such that 1k = − 1k+Q , where Q is the wave-vector of the antiferrodistortive mode X. This assumption actually overestimates the dimerisation strength, since it entails that any X 1 = 0 is able to open a hybridisation gap in the middle of band 1, which, we remark, does not coincide with the chemical potential unless band 2 is pushed above it. This implies that a finite hybridisation gap within band 1 does not stabilise an insulator so long as band 2 still crosses the Fermi energy. Therefore our simplified modelling does not spoil the important feature that a distorted insulating phase may occur only above a critical threshold of the Hamiltonian parameters, although it affects the value of that threshold, whose precise determination is however behind the scope of the present model-study.
In order to emphasise the bonding character of the a 1g , band 1, along the c R axis, as opposed to the more isotropic e π g , band 2, we choose the following forms of the two corresponding DOS's: with ∈ [−D, D] and N a normalisation factor. We take b > a/D 2 > 0 so that D 1 ( ) has a double-peak structure evocative of a one-dimensional DOS [71,149,153]. Hereafter, we take the half bandwidth D = 1 as our energy unit, and fix aD 3 = 1.9 and bD 5 = 2.1. The resulting DOS's are shown in Fig. 2 (a) and (b). There we note the two-peak structure of the band 1 DOS, mimicking the Van Hove singularities of a quasi one-dimensional band structure, in contrast to the structureless band 2 DOS. We highlight that the electron-electron interaction in Eq. (2) only includes the monopole Slater integral U > 0, and not higher order multipoles responsible of Hund's rules. This approximation, that makes the analysis more transparent, requires some justification. The Coulomb interaction of a single Vanadium projected onto the t 2g manifold, which effectively behaves as an l = 1 atomic shell, can be written in terms of two Slater integrals as: where n, S and L are the total occupation, spin, and angular momentum, respectively. Common values of the parameters are U t2g 4 eV and J t2g 0.68 eV U t2g /6 [24]. Denoting as E 0 (n) the lowest energy with n electrons in the t 2g shell, the effective Hubbard U for V 4+ can be defined through: to be compared with the VO 2 bandwidth of about 2.6 eV [56]. In units of the half-bandwidth, U 1.5, the value we shall use hereafter [154,155]. We observe that the Coulomb exchange J t2g has no effect on the configurations with n = 0, 1, while it splits those with n = 2 in three multiplets, with (S, L) = (0, 0), (1, 1), (0, 2), which are spread out over an energy J t2g , about a quarter of the full bandwidth. Such small value is not expected to qualitatively alter the physical behaviour, see, e.g., [156], which justifies our neglect of the exchange splitting in the model Hamiltonian (2).
We model the potential energy Φ X 1 , X 2 using a Landau functional for improper ferroelectrics [56,75,157] expanded up to the sixth order in the lattice displacements: where N is the number of sites and the couplings α to γ are all positive. The terms proportional to α, i.e. the harmonic part of the potential, and that proportional to γ have full rotational symmetry in the X 1 -X 2 plane. On the contrary, β 1 favours a lattice distortion only along one of the two components, whereas β 2 a distortion with In the specific case of VO 2 , β 2 > β 1 , and thus it is preferable to equally displace both modes [75] rather than just one of them. Finally, we write the electron-lattice coupling as: where c 1kσ creates an electron at momentum k in orbital 1 with spin σ, and we recall that, by construction, 1k = − 1k+Q . The dimerisation induced by the outof-plane displacement X 1 is controlled by the coupling constant g, while δ parametrises the strength of the crystal field splitting generated by the basal-plane displacement X 2 . The quadratic coupling in X 2 is intentional and has a physical explanation. Indeed, X 2 corresponds to the Vanadium displacement parallel to the diagonal of the rutile basal plane away from the centre of the Oxygen octahedron. As a result, the hybridisation between the e π g and the Oxygen ligands closer to the new Vanadium position increases, whereas the hybridisation with the further Oxygens diminishes. At linear order in the V-displacement X 2 , the two opposite variations of the hybridisation cancel each other, but, at second order, they add up to a net rise in energy of the e π g level, hence the expression in Eq. (7). The Hamiltonian Eq. (1) is invariant under the transformations X 1/2 → −X 1/2 , reflecting a Z 2 ×Z 2 (also known as K 4 or "Vierergruppe") symmetry.
Despite the great simplification effort, the model Hamiltonian Eq. (1) has still several parameters to be fixed. We emphasise that our main goal is to reproduce qualitatively the physics of VO 2 , without any ambition of getting also a quantitative agreement. Nonetheless, just to be sure not to explore a Hamiltonian parameter space completely detached from the real VO 2 compound, we choose parameters in line with the existing literature. We already mentioned our choice of U = 1.5, in units of the half-bandwidth, which is in line with the value used in realistic calculations [24,48,149,[158][159][160]. The other parameters involve the phonon variables. We shall choose g = 0.4, δ = 0.2, α = 0.155, β 1 = 1.75 · 10 −3 , β 2 = 2β 1 and γ = 6.722 · 10 −4 , which are compatible with various estimates of the electron-phonon coupling [135,161], of the lattice energy change across the rutile-to-monoclinic transition [162], as well as with direct experimental fits of those coupling constants [75,163].

III. DMFT SOLUTION
We solve the model Hamiltonian Eq. (1) by means of dynamical mean-field theory (DMFT) [164] within the adiabatic Born-Oppenheimer approximation. This approach will allow us to treat correlation effects nonperturbatively beyond an independent-particle description. Within DMFT, the lattice problem at a fixed displacement X = (X 1 , X 2 ) is mapped onto a quantum impurity model with an effective bath subject to a self-consistency condition. We solve the DMFT equations using as impurity solver exact diagonalization at zero and finite temperature [165,166], which requires discretisation of the effective bath in a finite number N b of levels. In this work, we take N b = 8, though we did test the validity of our results with respect to N b . We calculate the total electronic energy, or the free-energy at finite temperature, which renormalizes the Born-Oppenheimer adiabatic potential of the displacement Φ(X 1 , X 2 ) → Φ eff (X 1 , X 2 ) through: We shall restrict our analysis to the paramagnetic sector forcing spin SU (2) symmetry. However, we did check that magnetic solutions are higher in energy. We first present results at zero-temperature T = 0, and then move to those at T > 0.

A. Ground state phase diagram
In Fig. 3a we show the adiabatic potential Φ eff (X 1 , X 2 ) in (8) calculated by DMFT at U = 1.5. The energy landscape shows five minima. A local minima is located at the origin X 1 = X 2 = 0, and corresponds to an undistorted metal that we identify with the R phase of Vanadium dioxide. Four degenerate global minima are instead located at X 1 ±1.5 and X 2 ±2.1, which are related  to each other by the Z 2 ×Z 2 symmetry and represent the four equivalent lattice distortions. We find that these global minima describe an insulating phase, and thus realize a two-band version of the Goodenough scenario [38] for the M1 phase, in qualitative agreement with ab-initio calculations of VO 2 [145,146]. A detailed discussion of the electronic properties of all minima is postponed to the next Sec. III A 1.
In figures 3b and 3c we instead show the evolution of the adiabatic potential Φ eff (X 1 , X 2 ) along some specific lines, as indicated in Fig. 3a. We note that along the horizontal and vertical cuts, marked by a diamond and a circle in Fig. 3a, respectively, the energy landscape shows a saddle point, i.e., a minimum along the cut direction, but maximum in the perpendicular one. Within our model description, the effect of a uniaxial tensile strain would be taken into account by adding to the Hamiltonian Eq. (1) terms like: −F 1 X 2 1 or −F 2 X 2 2 (F 1 , F 2 > 0), depending on the direction of the applied stress [167][168][169]. In presence of such terms, the saddle points observed in Fig. 3a along the lines X 1 = 0 or X 2 = 0 may turn into additional minima of the energy landscape [147], which can possibly describe the occurrence of the M2 phase in the framework of the same model Hamiltonian.
In order to understand what is the role of the Hubbard interaction U in stabilising the insulating solution, we studied the evolution of Φ eff (X 1 , X 2 ) for several values of U , along the line in the X 1 -X 2 plane connecting the rutile local minimum with one of the monoclinic global minima (the diagonal cut in Fig. 3a marked by a dia- mond symbol). Our results are reported in Fig. 4. We note that already at U = 0 the energy has two minima. One is at the origin and corresponds to the undistorted metal. The other is located at finite X, and thus represents a distorted phase that must evidently be also insulating in order to be a local energy minimum. Therefore at small U 0.2, the stable phase is the undistorted metal at X = 0 in Fig. 4, while the local minimum at X = 0 (monoclinic insulator) is metastable. However, for larger U 0.2, the situation is reversed: the distorted insulator becomes the global minimum, while the undistorted metal a local one, entailing the typical scenario of a first-order metal-insulator transition driven by interaction. The above results show that electron-electron interaction is crucial to stabilise the distorted insulator, though the active contribution of the lattice is equally essential. Indeed, the interaction strength, U 1.5 the half-bandwidth, is too small to drive on its own the metal-insulator transition [156]. In other words, the picture that emerges from Fig. 4, with the interaction and the coupling to the lattice both necessary to stabilise the insulator, fully confirm our expectation in Sec. I.

Spectral functions
Further insights into the properties of the metalinsulator transition can be gained by looking at the spectral functions: where a = 1, 2 and G loc,aa is the local interacting Green's function obtained within the DMFT solution of the model. In Fig. 5 we show A a (ω) at the different minima in 3a, with ω measured with respect to the chemical potential. We note that already in the absence of interaction, U = 0, the different shapes of the DOS's, see Fig. 2, lead to a larger occupation of band 1 than band 2. Such population unbalance is increased by U > 0, which effectively enhances the crystal field, leading to an even larger occupation of band 1 at expenses of 2 [151,170,171]. This is evident in the spectral function of the undistorted metal at U = 1.5, reported in Fig. 5(a) and Fig. 5(b), where the occupied ω ≤ 0 part of A 1 (ω) overwhelms that of A 2 (ω) more than in the U = 0 case of Fig. 2.
We also note in the figures 5(a) and 5(b) side peaks that correspond to the Hubbard bands. The scenario is radically different in the insulating solution, see Fig. 5(c) and Fig. 5(d). Here we observe the formation of a hybridisation gap opening at the chemical potential inside the band 1. Two coherent-like features flank the gap. The band 2 is instead pushed above the Fermi energy, and therefore is empty. We still observe Hubbard sidebands in A 1 (ω), as well as signatures of the upper Hubbard band in A 2 (ω), though rather spiky because of the bath discretisation.
We note that in the insulating solution the lowest gap corresponds to transferring one electron from band 1 to band 2, i.e., from a 1g to e π g in the VO 2 language, and has a magnitude of about E gap ∼ 0.8 eV, for a realistic value of the half-bandwidth of 1.3 eV [56]. This value of the gap is not too far from the experimental one, E ex gap ∼ 0.6 − 0.7 eV [35,52,118]. Therefore, our simplified modelling yields results that are not only qualitatively correct but, rather unexpectedly, also quantitatively not far off the actual ones. The band 1 → band 1 transition, i.e., a 1g → a 1g , though being slightly higher in energy, has a much steeper absorption edge since it involves the two coherent peaks in Fig. 5(c). This result is in loose agreement with XAS linear dichroism experiments [86,172] that are able to distinguish the two absorption processes.
In order to assess the degree of electronic correlations, we calculate the quasiparticle residue of each band in the undistorted metal phase, defined by: with a = 1, 2. We find that the two bands at X 1 = X 2 = 0 show almost the same value Z 1/2 ∼ 0.67, not inconsistent with more realistic calculations [24,72,159,173]. Such agreement, a priori not guaranteed, gives further support to our simple modelling.

IV. PHASE TRANSITION AT FINITE TEMPERATURE
Our main scope here is however to describe the first-order phase transition upon heating from the low-temperature M1 monoclinic insulator to the hightemperature rutile metal. In general, we can envisage a phase transition primarily driven either by the electron entropy or by the lattice one.
Indeed, we note that the electron free energy of the metal solution, which is metastable at T = 0, must drop faster upon raising temperature than the insulator free energy since the metal carries more electron entropy than the insulator. This effect alone, that is ignoring lattice entropy, would be able to drive a first-order transition when insulator and metal free energies cross. On the other hand, since the distorted ground state breaks the Z 2 × Z 2 symmetry of the adiabatic lattice potential Φ eff (X 1 , X 2 ) in Fig. 3, we might expect such symmetry to be recovered by raising temperature only because of lattice entropy effects, i.e., ignoring the electronic contribution to entropy.
In reality, both effects should combine to drive the transition. However, dealing together with lattice and electron entropies within our computational scheme would imply to calculate the adiabatic potential Φ eff (X 1 , X 2 ) at any temperature, which is a rather heavy task. For this reason, in what follows we shall analyse separately electron and lattice entropy effects, and at the end argue what would happen should they act together.

A. Electron-driven transition
Let us first neglect the lattice entropy and study the temperature evolution of the free energies of the two inequivalent minima that we found at zero temperature. For that, we need to evaluate the electronic entropy, which can be obtained through: The last equality corresponds to a change of integration variable from the temperature T to the adiabatic potential Φ eff , which is also the internal energy. From the entropy S we can estimate the free energy: which, we emphasise once more, does not include the lattice contribution to entropy. We shall compare the free energy of the undistorted metal solution at X = 0, with that of the distorted insulator at X = 0. In principle, the equilibrium displacement in the insulator should change with temperature. In practice, since the entropy of the insulator is negligible for all temperatures under consideration, we shall fix X at the T = 0 value. The temperature evolution of the metal and insulator free energies so obtained are shown in Fig. 6. As expected, the larger entropy of the metal pushes its free energy below the insulator one at relatively low temperature, T el ∼ 0.021, Figure 6. (Color online) Temperature evolution of the free energy at the two inequivalent minima X1 = X2 = 0 (dots) and X1 = 1.5, X2 = 2.1 (squares) observed at zero temperature for U = 1.5. The first-order transition occurs at T el ∼ 0.021 ∼ 320 K, of the same order of magnitude as the experimental value 340 K.
substantially smaller than the insulating gap, and thus justifying our assumption of frozen X. T el identifies the insulator-metal transition, which is evidently first order since the two free energies cross with different slopes. Incidentally, T el ∼ 0.021 in half-bandwidth units, corresponds to ∼ 320 K for a realistic bandwidth of 2.6 eV, which has the right order of magnitude when compared with the true critical temperature of 340 K.

B. Lattice driven transition
We now move to study the properties of the latticedriven transition. For that, we first need to model the lattice dynamics. However, since the tetragonal R to monoclinic M1 transition is a complex structural transformation, with martensitic features, especially in films [174][175][176][177][178], our modelling ought to be oversimplified, and aimed just to get qualitatively reasonable results, with no pretension of quantitative accuracy.
As a first step, we must relax our previous assumption of a global antiferrodistortive mode, and instead introduce a displacement field, i.e., a site dependent displacement X i = X 1i , X 2i . We assume that X i feels the local adiabatic potential Φ eff X i of Fig. 3a, temperature independent since we are neglecting the electron entropy. In addition, we suppose that the displacements of nearest-neighbour sites are coupled to each other by an SO(2) ∼ = U (1) invariant term that tends to minimise the strain. With those assumptions the classical Hamiltonian 0.0 0.5 reads: where X denotes a configuration of all the displacement vectors. The model (13) is equivalent to a generalized XY -model, where X i plays the role of two-component spin of variable length, while J > 0 is the conventional spin stiffness. Both length and direction of X i are controlled by the anisotropic potential Φ eff X i , which is not invariant under U (1) but under separate X 1 → −X 1 and X 2 → −X 2 transformations, i.e., Z 2 × Z 2 . The phase diagram of an XY model in presence of an anisotropy term that lowers U (1) down to Z n is already known [179][180][181]. In particular, the anisotropy Z n for n ≥ 4 is a dangerously irrelevant perturbation that does not change the XY universality class of the transition [180,181]. Our specific case study, where U (1) → Z 2 × Z 2 , has not been considered yet, at least to our knowledge, but it should most likely change the XY universality class, which is what we are going to investigate in the following.
We study the classical model Eq. (13) at different temperatures by the standard Monte Carlo (MC) method [182] on a three-dimensional cubic lattice of side N x .
In Fig. 7(a) we plot the modulus of the average displacement, X , as function of the temperature. For small system size (e.g. N x = 10) X shows a smooth crossover in temperature. However, increasing N x unveils the existence of a continuous phase-transition at a critical value T c of the temperature, which is controlled by the value of J because Φ eff has been calculated earlier.
Since J is unknown, we have preferred to use T c as the unit of temperature in Fig. 7 and in those that follow. In order to better reveal the second order character of the transition, we also show in Fig. 7(a) the fit with a mean-field square-root behaviour. The fit is rather good, although we known that close to the transition the actual critical behaviour must deviate from mean-field.
A closer look to the temperature dependence of the order parameter uncovers a non-trivial two-step evolution, which is more evident in Fig. 7(b), where we show the specific heat C v = ∂ E /∂T vs. T . Indeed, C v clearly displays two peaks that are suggestive of two distinct transitions. The first transition at T = T c , below which X acquires a finite value, is followed by a second one at lower T = T d < T c .
In order to understand the nature of both transitions, in Fig. 8 we show at T > T c , left panels, T d < T < T c , middle panels, and T < T d , right panels, the endpoint distribution after N s = 4 × 10 5 MC sweeps of the lattice of the N 3 x displacement vectors superimposed to the potential landscape in the (X 1 , X 2 ) space (top panels), and a real space snapshot within a single layer of the cubic lattice (bottom panels). At high temperature, T > T c , the X i 's cover homogeneously the whole potential landscape, see top-left panel, without any appreciable spatial correlation, see the bottom-left panel. Lowering T slightly below T c , we observe a significant change in the displacement distribution, see middle panels. Specifically, the system seems to break ergodicity first along X 2 , in the simulation corresponding to the figure it localises in the X 2 > 0 half-plane, while it is still uniform along X 1 . Consequently, clusters of parallel displacement vectors form in real space. The alignment direction has X 2 > 0 for all clusters, while the X 1 component changes from cluster to cluster, see bottom-middle panel. Only below T d , full ergodicity breakdown occurs, with the system trapped around just one of the four equivalent minima, in the figure that with X 2 > 0 and X 1 > 0. In other words, the Z 2 × Z 2 symmetry of the model Eq. (13) gets broken in two steps upon cooling: first, the Z 2 symmetry X 2 → −X 2 spontaneously breaks, and next, the residual X 1 → −X 1 symmetry, leading to two consecutive Isinglike transitions. This is summarised in Fig. 9, where we see that at T c X 2 becomes finite, and thus also X , while X 1 is still zero. Only below T d also X 1 acquires a finite average value.
Translated in the language of VO 2 , these results suggest the existence of an intermediate monoclinic phase for T d < T < T c where the V atoms are displaced only within the basal plane, i.e., the chains are tilted but not yet dimerised. In our model Hamiltonian (1), such phase with X 1 = 0 describes a monoclinic metal, which, as discussed in Sec. I, has been reported in several experiments [81][82][83][125][126][127][128][129][130][131][132][133]. Only below T d , both components of the antiferrodistortive displacement are finite, leading to the M1 insulating phase.
In conclusion, without including the electron entropy x endpoints of the calculation. Bottom panels: displacement field configuration within a single plane of the cubic lattice, with the same parameters of the top panels. If Xi = Xi (cos θi, sin θi), the color code represents θi ∈ [0, 2π], and the arrow length Xi . At high temperature T > Tc (left panels) Xi have random length and orientation, thus covering homogeneously the entire potential landscape. For T d < T < Tc (center panels) the displacement orientation shows breaking of the Z2 symmetry X2 → −X2. At lower temperature T < T d < Tc, also the residual Z2 symmetry X1 → −X1 gets broken; most of the Xi's have length and direction corresponding to just one of the potential global minima.
we find two transitions that look continuous and in the Ising universality class: one at T d between a monoclinic insulator and a monoclinic metal, and another at T c > T d from the monoclinic metal to a rutile one. On the contrary, neglecting the lattice entropy and just including the electronic one, we found in Sec. IV A a single first-order transition at T el , directly from the monoclinic insulator to the rutile metal. We can try now to argue what we could have obtained keeping both entropy contributions still within the Born-Oppenheimer adiabatic approximation. Evidently, if T el T c the scenario should not change qualitatively with respect to the two-transition one uncovered in this section. We cannot exclude that the electron entropy and all the lattice effects we did not include in the simple model (13) could change the transitions into first-order ones, but we do expect still two distinct transitions. On the contrary, if T el < T d we would predict a single first-order transition like in Sec. IV A.
The experimental evidence supporting the existence of a monoclinic metal phase intruding between the M1 insulator and R metal [81][82][83][125][126][127][128][129][130][131][132][133] suggest that, should our modelling be indeed representative of VO 2 , then the Hamiltonian parameters should be such that T el T c . We also observe by comparing Fig. 6 with 7 that the loss of lattice entropy upon cooling across the transitions overwhelms that of electron entropy, suggesting lattice driven transitions in agreement with experimental [89] and theoretical [183] proposals.

V. CONCLUSIONS
We have constructed a minimal model that we believe contains all essential ingredients to correctly capture the physics of the metal-insulator transition in vanadium dioxide.
The model comprises two orbitals per site, one mimicking the a 1g and the other the e π g , thus neglecting the twofold nature of the latter, which broaden into two bands. The a 1g band has a double peak structure reflecting its bonding character along the rutile c-axis, while the e π g one is structureless. Both have the same bandwidth and centre of gravity. The density corresponds to one electron per site, i.e., the two bands are at quarter filling. The electrons feel an on-site Hubbard repulsion, and are coupled to two zone-boundary lattice modes, corresponding, respectively, to the basal plane component, i.e., the tilting of the Vanadium chains, and out-of-plane component, responsible of the chain dimerisation, of the antiferrodistortive displacement that acquires a finite expectation value below the transition from the high temperature rutile structure to the low temperature monoclinic one (M1). Using realistic Hamiltonian parameters and assuming the Born-Oppenheimer adiabatic approximation, we find at low temperatures phase coexistence between a stable distorted insulator, the monoclinic M1 insulator, and a metastable undistorted metal, the rutile metal. Upon rising temperature, we observe a two-step transition. First, the dimerisation component of the antiferrodistortive displacement melts, leading to a transition from the monoclinic insulator to a monoclinic metal. At higher temperature also the tilting component disappears, and the monoclinic metal turns into the rutile one. Such a two-transition scenario, not in disagreement with experiments, is mostly driven by the lattice entropy, also in accordance with experiments. One of the messages of our model calculation is that the electron-electron interaction has the role to effectively enhance the coupling to the lattice, stabilising a distorted phase otherwise metastable in the absence of interaction. This also implies that we could have obtained similar results with weaker electronic correlations but stronger electron-lattice coupling. This conclusion is actually supported by the phenomenology of Niobium dioxide NbO 2 , which, mutatis mutandis, is akin to that of VO 2 . NbO 2 also undergoes a metal-insulator transition, though at substantially higher temperature of T MIT ∼ 1080 K [184][185][186][187]. This transition occurs slightly below a structural one at T s ∼ 1123 K [188], from a high-temperature rutile structure to a low-temperature body centred tetragonal one that locally resembles the M1 phase of VO 2 [189][190][191][192]. However, the single 4d-electron in Nb 4+ is expected to be less correlated than the 3d-electron in V 4+ . This loss of correlations, testified by the VO 2 M2 phase having no counterpart in NbO 2 [193], and by the efficacy of ab initio methods to describe NbO 2 [74,[194][195][196], is actually overcompensated by the increase in covalency due to the broader spatial distribution of the 4d orbitals [197], which, in turn, yields a stronger coupling with the zone-boundary lattice modes, and thus a higher transition temperature. The intermediate poorly metallic phase for T MIT < T < T s [188,189,198,199] is thus the counterpart of the monoclinic metal in VO 2 , although the former is more clearly established than the latter.