Mott Transition in the Hubbard Model on Anisotropic Honeycomb Lattice with Implications for Strained Graphene: Gutzwiller Variational Study

The modification of interatomic distances due to high pressure leads to exotic phenomena, including metallicity, superconductivity and magnetism, observed in materials not showing such properties in normal conditions. In two-dimensional crystals, such as graphene, atomic bond lengths can be modified by more than 10 percent by applying in-plane strain, i.e., without generating high pressure in the bulk. In this work, we study the strain-induced Mott transition on a honeycomb lattice by using computationally inexpensive techniques, including the Gutzwiller Wave Function (GWF) and different variants of Gutzwiller Approximation (GA), obtaining the lower and upper bounds for the critical Hubbard repulsion (U) of electrons. For uniaxial strain in the armchair direction, the band gap is absent, and electron correlations play a dominant role. A significant reduction in the critical Hubbard U is predicted. Model considerations are mapped onto the tight-binding Hamiltonian for monolayer graphene by the auxiliary Su–Schrieffer–Heeger model for acoustic phonons, assuming zero stress in the direction perpendicular to the strain applied. Our results suggest that graphene, although staying in the semimetallic phase even for extremely high uniaxial strains, may show measurable signatures of electron correlations, such as the band narrowing and the reduction in double occupancies.

A notable exception, however, is a honeycomb lattice, for which relatively simple techniques, including GWF [15] or CPA [16,17], provide reasonable approximations for the critical Hubbard interaction, differing from the QMC value, U c = 3.86(1) t 0 [18] (with t 0 being the nearest-neighbor hopping integral and the number in parenthesis denoting uncertainty for the last digit), by less than 10%.For a comparison, the HF method gives U (HF) c = 2.23 t 0 [19] for the same lattice.
Since the advent of graphene [20,21] a half-filled, fermionic honeycomb-lattice systems have attracted renewed attention, as they emulate several field-theoretical phenomena in condensed matter [22].The effective Hubbard model for monolayer graphene with on-site interaction U eff ≈ 1.6 t 0 was proposed [23], suggesting that large isotropic strain may drive this system from semimetallic towards Mott-insulating phase [24,25] in analogy with high pressure changing properties of various bulk materials [26][27][28][29][30]. Effects of electron correlations are usually more pronounced in graphene nanosystems, where quantum fluctuations are reduced and magnetic moments may form near free edges [31][32][33] (although defining metallic and insulating states for a nanosystem is more cumbersome than for a bulk system [34,35]).We further notice that artificial graphene-like systems allow one to tune the interaction in a wider range than actual graphene [36][37][38][39].Yet another possibility to study electron correlations has open with the fabrication of twisted bilayer graphene [40,41].
A separate issue concerns the bandgap opening due to spatial rearrangement of atoms in strained graphene (a so-called two-dimensional Peierls instability), which may turn the system into insulator before the Mott transition occurs [42][43][44][45][46].To the contrary, weak electron-phonon interaction of the Holstein type, which may appear in graphene on some substrates, is predicted to favor the semimetallic phase [47].
In this paper, the discussion is limited to a honeycomb lattice strained along main crystallographic axes (see Fig. 1) supposing that the bipartite structure of the lattice is preserved under strain.In turn, there are two different values of the nearestneighbor hopping integral in a single-particle Hamiltonian, t x and t y , corresponding to electron hopping along the zigzag direction (t x ) or along the armchair direction (t y ).To obtain a direct mapping between the strain applied and the hopping integrals, a version of the Su-Schrieffer-Heeger (SSH) model [48] is developed, with microscopic parameters adjusted to match elastic properties of graphene [49].Once fixed strain is applied in a selected direction, the lattice is allowed to relax along the perpendicular direction to reach a conditional energy minimum (a zero perpendicular stress case).We further focus our attention on the strain applied along armchair direction, for which the system evolves towards a collection of weakly-coupled one dimensional chains [50][51][52] allowing one to expect that, once the effective Hubbard model is considered, the Mott transition may appear for smaller value of U eff than for isotropic strain.
The remaining part of the paper is organized as follows.In Sec.II, we briefly present approximate approaches to the effective Hubbard Hamiltonian (including HF, GWF, and GA).Also in Sec.II, we show some original data, illustrating how these approaches work for anisotropic honeycomb lattice.In Sec.III, we discuss our numerical results concerning the phase < l a t e x i t s h a 1 _ b a s e 6 4 = " p G z q z T p O R k 8 e A E a 2 9 R q X 8 l k 3 A Q E = " > A A A B 7 n i c b V D L S g M x F L 1 T X 7 U + O u r S T b A I r s q M K L q S o h u X F e w D 2 l I y 6 Z 0 2 N D M T k k y x l P 6 G G x E 3 C n 6 K v + D f m L a z a e u B w M k 5 J + S e G 0 j B t f G 8 X y e 3 s b m 1 v Z P f L e z t H x w W 3 a P j u k 5 S x b D G E p G o Z k A 1 C h 5 j z X A j s C k V 0 i g Q 2 A i G D z O / M U K l e R I / m 7 H E T k T 7 M Q 8 5 o 8 Z K X b f Y H l G F U n N h b y 9 3 X t c t e W V v D r J O / I y U I E O 1 6 / 6 0 e w l L I 4 w N E 1 T r l u 9 J 0 5 l Q Z T g T O C 2 0 U 4 2 S s i H t 4 2 Q + 7 p S c W 6 l H w k T Z E x s y V 5 d y N N J 6 H A U 2 G V E z 0 K v e T P z P a 6 U m v O 1 M e C x T g z F b f B S m g p i E z L q T H l f I j B h b Q p n i d k L C B l R R Z u y G C r a 6 v 1 p 0 n d Q v y / 5 1 2 X u 6 K l X u s y X k 4 R T O 4 A J 8 u I E K P E I V a s A g h T f 4 h C 9 H O q / O u / O x i O a c 7 M 0 J L M H 5 / g O F W 4 7 U < / l a t e x i t > " x > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " s r 7 T h 5 K 2 E X D N s X S T i u f N o E 4 a V 6 E = " > A A A B 7 n i c b V D L S g M x F L 3 j s R s S N g p / i L / g 3 p u 1 s 2 n o g c H L O C b n n B l J w b T z v 1 9 n Y 3 N r e 2 S 3 s F f c P D o 9 K 7 v F J Q y e p Y l h n i U h U K 6 A a B Y + x b r g R 2 J I K a R Q I b A a j x 5 n f H K P S P I l f T C a x G 9 F B z E P O q L F S z y 1 1 x l S h 1 F z Y W 3 b v 9 d y y V / H m I O v E z 0 k Z c t R 6 7 k + n n 7 A 0 w t g w Q b V u + 5 4 0 3 Q l V h j O B 0 2 I n 1 S g p G 9 E B T u b j T s m F l f o k T J Q 9 s S F z d S l H I 6 2 z K L D J i J q h X v V m 4 n 9 e O z X h X X f C Y 5 k a j N n i o z A V x C R k 1 p 3 0 u U J m R G Y J Z Y r b C Q k b U k W Z s R s q 2 u r + a t F 1 0 r i q + D c V 7 / m 6 X H 3 I l 1 C A M z i H S / D h F q r w B D W o A 4 M U 3 u A T v h z p v D r v z s c i u u H k b 0 5 h C c 7 3 H 4 b X j t U = < / l a t e x i t > " y > 0 i < l a t e x i t s h a 1 _ b a s e 6 4 = " u a F 9 y A n v w v d G D C E O c m e U i b e x T K 8 = " > A A A B 3 n i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h d 0 o 6 D H o x W M C 5 g F J i L O T 3 m T I 7 I O Z X i G E X L 2 I e F H w k / w F / 8 Z J s p c k F g w U V T V 0 V / u J k o Z c 9 9 f Z 2 N z a 3 t n N 7 e X 3 D w 6 P j g s n p w 0 T p 1 p g X c Q q 1 i 2 f G 1 Q y w j p J U t h K N P L Q V 9 j 0 R w 8 z v / m C 2 s g 4 e q J x g t 2 Q D y I Z S M H J S j X Z K x T d k j s H W y d e R o q Q o d o r / H T 6 s U h D j E g o b k z b c x P q T r g m K R R O 8 5 3 U Y M L F i A 9 w M l 9 v y i 6 t 1 G d B r O 2 L i M 3 V p R w P j R m H v k 2 G n I Z m 1 Z u J / 3 n t l I K 7 7 4 Y e a 3 X l B p H s t n M 0 m w F 9 G B 5 C F n 1 F i p G f Q z P p r 2 S 2 W 3 4 s 5 B 1 o m X k z L k q P d L P 9 0 g Z m m E 0 j B B t e 5 4 b m J 6 G V W G M 4 H T Y j f V m F A 2 p g P M 5 j t O y a W V A h L G y j 5 p y F x d y t F I 6 0 n k 2 2 R E z V C v e j P x P 6 + T m v C u l 3 G Z p A Y l W w w K U 0 F M T G a F S c A V M i M m l l C m u N 2 Q s C F V l B l 7 l q K t 7 q 0 W X S f N a s W 7 r l S f b s q 1 + / w I B T i H C 7 g C D 2 6 h B o 9 Q h w Y w G M E b f M K X g 8 6 r 8 + 5 8 L K I b T v 7 n D J b g f P 8 B r 7 2 K 8 w = = < / l a t e x i t > ] ( j ) < l a t e x i t s h a 1 _ b a s e 6 : Top: Honeycomb lattice subjected to uniaxial strain in selected direction (see the coordinate system).Zoom-in visualizes the distance (bond length) dij between atoms i and j and in-plane angle with the vertex at site j ( (j)).Bottom: Hexagonal first Brillouin zone (FBZ) of the reciprocal lattice, with (dimensionless) basis vectors b1 = 2π/ √ 3 ( √ 3, −1) and b2 = 2π/ √ 3 (0, 2), and the symmetry points, K = (4π/3, 0) and K = (2π/3, 2π/ √ 3) coinciding with Dirac points in the absence of strain.Magnified area shows discretized FBZ for a finite system of N = 2NxNy atoms with periodic boundary conditions [see Eq. diagram of the effective Hubbard model with arbitrary parameters (t x t y and U eff ), the evolution of the model parameters in graphene subjected to uniaxial strain, and approximate formula relating the reduction of U c to strain-induced anisotropy of the Fermi velocity.The effects of electron correlations on selected measurable quantities are also presented in Sec.III.The concluding remarks are given in Sec.IV.
Next to the main text, in Appendix A, the Coherent Potential Approximation (CPA) is briefly described.In Appendix B, we present the auxiliary SSH model, proposed to relate the physical strain onto the microscopic parameters of the effective Hubbard model.1) with U = 0 displayed as a function of energy.Top: strain applied in the armchair direction (ty tx), bottom: strain applied in the zigzag direction (ty tx).The ratio t</t> [with t< = min (tx, ty) and t> = max (tx, ty)] is varied between the lines with the the steps of 0.2.A vertical offset is applied to each dataset except from the isotropic case (tx = ty).Inset shows the band gap, appearing for tx < 0.5 ty due to the Peierls transition.

A. The anisotropic Hubbard model
Our analysis of electron correlations on anisotropic honeycomb lattice starts from the Hamiltonian with the first sum running over pairs of nearest-neighbors ij and spin up/down orientations (s =↑, ↓), and the hoppingmatrix elements are given by (Without loss of generality, we suppose the coordinate system is oriented as depicted in Fig. 1.) Remaining symbols in Eq.
(1) are a creation (annihilation) operator for electron with spin s on the lattice site i, c † i,s (c i,s ), n is = c † i,s c i,s , and the on-site Hubbard repulsion U .We further limit our considerations to the ground state and suppose the half-filling, i.e., one electron per lattice site, n = n i↑ + n i↓ = 1.
In principle, ground-state properties of the model defined by Eqs. ( 1) and ( 2) can be discussed as functions of two dimensionless parameters, e.g., t y /t x and U/t x .The relation between parameters t x and t y and strain applied to graphene is discussed later in this section.But first, we briefly present approximate approaches capable to distinguish whether ground state of the Hamiltonian (1) is semimetallic or insulating.

B. Hartree-Fock approximation
Although a honeycomb lattice is bipartite and the antiferromagnetic order is possible, its peculiar band structure suppresses antiferromagnetism at small U [15].Since a single particle density of states (i.e, density of states at U = 0) is linear for low energies, see Fig. 2, there is no Fermi surface that could produce magnetic instability also for small U > 0.
Within the Hartree-Fock approximation, interaction part in the Hamiltonian (1) is replaced by where we have introduced the operator D = j n j↑ n j↓ measuring the number of double occupancies.We further impose the antiferromagnetic order, where λ i = 1 if i belongs to one sublattice (A), or λ i = −1 if i belongs to the other sublattice (B), and m is the magnetization (|m| n), and the half filling (n = 1).The above yields the HF ground-state energy per site where the factor 2 accounts for s =↑, ↓, and the summation runs over quasimomenta k ≡ (k x , k y ) in the first Brillouin zone, namely with N x,y being the number of unit cells in x, y direction, N = 2N x N y (the periodic boundary conditions are imposed).
For sufficiently large number of points in the momentum space, say N x , N y 10 3 , one can usually works with a square inverse lattice (omitting the term ∝ n x in the expression for k y ) [53]; nevertheless, the discretization of (k x , k y ) as given in Eq. ( 6) becomes crucial when discussing the finite-size effects for small N .The single-particle energies for anisotropic honeycomb lattice are given by with Next, the density of states is defined as with two parts corresponding to the conduction (E > 0) and valence (E < 0) band, and satisfying the normalization conditions: given by Eq. ( 5) brought us to In case the solution of Eq. ( 10) does not exist, the minimum of E [15].This can be easily understood for the case of unstrained (or uniformly strained) lattice, for which t x = t y = t 0 and the density of states can be approximated by with a cut-off energy of Λ = √ 3π 1/2 t 0 2.33268 t 0 .The above is equivalent, for |E| Λ, to ρ Λ (E) = 2A|E|/ N π( v F ) 2 , with A being the system area, v F = Generalized Gutzwiller wavefunction, allowing antiferromagnetic order, was applied in Ref. [15] to find out that correlated, but paramagnetic solution remains stable up to the region of the Mott semimetal-insulator transition.Although several features of the solution are altered when employing more advanced techniques [18], the values of U c following from < l a t e x i t s h a 1 _ b a s e 6 4 = " y S 1 D 5 Q R t x u i P j  < l a t e x i t s h a 1 _ b a s e 6 4 = " z 7 6 o 2 Z F w a 7 4 Y q x 5 J + q i 6 A W 9 o U / 0 p T 1 r r 9 q 7 9 j G K z m n j P z t o A t r 3 H x l 3 k 3 o = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " z 7 6 o 2 Z F w a 7 4 Y q x 5 J + x < l a t e x i t s h a 1 _ b a s e 6 4 = " v + h P w y r h M Y 5 I 5 s v 5 9 given in Table I. (Statistical errorbars are to small to be shown on the plots.)GWF are surprisingly close to those obtained within largescale computer simulations for a honeycomb lattice.Investigating the variational wavefunction where |ψ 0 (m) denotes a Slater determinant corresponding to a given magnetization in Eq. ( 5) and η is another variational parameter (quantifying the role of electron correlations), one needs to minimize the ground-state energy with respect to η and m.In many cases, the system may prefer to reduce m (even to m = 0) and increase η, allowing to expect that, in general, U . Several approximated techniques for calculating the averages in Eq. ( 13) were developed [7][8][9][10]15].Here, we apply Variational Monte Carlo (VMC), described in details in Ref. [10].To determine the value of U (GWF) c , we have directly followed the procedure proposed by Martelo et al. [15].For a fixed value of the gap (∆ ≡ U m), the energy difference , where the parameter η is optimized independently for m = 0 and m = 0, changes sign at some U = U 0 (∆).Numerical extrapolation of U 0 (∆) with ∆ → 0 allows one to determine the critical value of U (GWF) c . Selected examples, for t y t x (i.e., strain applied in the armchair direction) and the system size of N = 200 sites (N x = N y = 10), are presented in Fig. 3.
For more details of the simulation, see Ref. [55].

D. Gutzwiller Approximation and its variants
To efficiently study the effects of electron correlations present in |Ψ GWF , see Eq. ( 12), one can also adopt the Gutzwiller Approximation (GA) and find out how the number of double occupancies is reduced comparing to the HF solution |ψ 0 (m) .Within GA, which is exact in the infinite dimension limit, the correlation functions where the band-narrowing factor q({ρ (0) ll s }; {d l }) depends only on the single-particle density matrix elements ρ (0) ll s = c † ls c l s 0 with l, l = i, j and s =↑, ↓ (here, . . .0 is the expectation value over the uncorrelated state; i.e., a single Slater determinant such as |ψ 0 (m) ), and d l = n l↑ n l↓ GWF (l = i, j) being the average double occupancies.The {d i } variables are further regarded as variational parameters to be determined by minimizing the Gutzwiller energy functional, Several forms of the band-narrowing factor q({ρ (0) ll s }; {d l }), being equivalent in the infinite dimension limit but producing slightly different results when applied to the system of a finite dimensionality, are used among the literature [8,[56][57][58][59][60][61][62].For the diagonal elements ρ (0) iis = n is 0 parametrized as in Eq. ( 4) with n = 1, one can impose d i ≡ d for all sites and rewrite the expression given in  ) bound following from the Gutzwiller Approximation (GA) and the Néel-state Gutzwiller Approximation (NGA).The results obtained from the Statistically-consistent Gutzwiller Approximation (SGA) are also given.The system size is defined by Nx = Ny = 10 for VMC simulations; remaining results correspond to the limit of Nx = Ny → ∞.Ref. [61] as The variable d is bounded as 0 d .This brought us to Numerical minimization of , with respect to (m, d), truncates the optimization of both the density matrix {ρ (0) i,j } and the parameters {d i }.For the linear density of states ρ Λ (E), see Eq. ( 11), one can easily find closed-from expression for E (GA) G ; the minima corresponding to m = 0 appear for U > U (Λ,GA) c = 1.270Λ = 2.963 t 0 , with the critical value lying between HF [19] and QMC [18] results for isotropic honeycomb lattice.
A slightly more accurate (but also more computationally expensive) approach can be constituted by parametrizing the uncorrelated state |ψ 0 not only via the magnetization m, as in the above, but via all independent parameters of the density matrix {ρ (0) i,j }.In particular, the auxiliary single-particle Hamiltonian determining {ρ (0) i,j } contains the renormalized hopping integrals ( tx and ty ) which may differ from t x and t y in the multiparticle Hamiltonian (1).The resulting method, called the Statistically-consistent Gutzwiller Approximation (SGA), is presented in details in Ref. [58].
Both (S)GA and GWF methods can be regarded as improvements to mean-field (HF) solution, including some classes of quantum fluctuations.Since not all fluctuations are included, the AF order is artificially favored when searching for the energy minimum, and therefore these methods usually underestimate the value of U c .In order to bound U c from the top, we employ the scheme proposed by Martelo et al. [15], in which two solutions are compared: The paramagnetic GA solution, corresponding m = 0 in Eq. ( 17), with a complementary variational wavefunction, where κ is a variational parameter, T = ij ,s t ij (c † is c js + H.c.) is the kinetic-energy part of the Hamiltonian (1), and |Ψ U →∞ is the ground state for U → ∞.The critical value U c is than estimated by finding a crossing point of E (GA) G , Eq. < l a t e x i t s h a 1 _ b a s e 6 4 = " T F F i + N g S T B S z G y I y x A o T Y 8 + T t 9 W 9 5 a K r p H 5 V 9 q 7 L l a d K s X q R H S E H p 3 A G J f D g B q r w C D X w g Y C C N / i E L 0 c 4 r 8 6 7 8 z G P r j n Z n x N Y g P P 9 B / r H j T 8 = < / l a t e x i t > U (GA) c < l a t e x i t s h a 1 _ b a s e 6 4 = " P q m w v m P u r f Z t G U s w y I w y S I M 0 J r c = " , see Table I.] (17), with a fixed m = 0 and optimized d, and the variational energy E B corresponding |Ψ B , Eq. ( 18), with optimized κ.
In the limit of infinite dimensions, the variational energy E B associated with the state |Ψ B can be evaluated exactly, since the ground state |Ψ U →∞ is the Néel antiferromagnet [65].The variational energy reads where are the kinetic energy per site and the sublattice magnetization (respectively).The so-called Néel-Gutzwiller Approximation (NGA) is constituted by substituting the density of states given by Eq. ( 9) into Eqs.( 21), (22), and the subsequent minimization of with respect to κ. Selected numerical results, for t y t x , are presented in Fig. 4.
It is worth mentioning that (S)GA can be systematically improved, approaching the GWF solution, by including consecutive corrections following from the relevant diagrammatic expansion [8,60,62] (we further notice that the detailed scheme for a honeycomb lattice is missing so far).Similar approach for |Ψ B is difficult due to necessity of determining the ground state of the Heisenberg model (|Ψ U →∞ ) as a first.
Selected numerical values of U c , following from the methods described in this Section, are compared in Table I.
A substantially different approach, the Coherent Potential Approximation (CPA), in which one considers random scattering of electrons with a given spin on motionless electrons with the opposite spin (instead of imposing some spin order), is described in Appendix A.

A. Phase diagram
Our central results are presented in Fig. 5, where we display the phase diagram for the Hamiltonian (1) with the strain applied in armchair direction (t y t x ).A single-particle spectrum is gapless in such a case (see also Fig. 2) since the positions of Dirac cones do not merge [66]; therefore, metalinsulator (if occurs) must be driven by electron-electron interaction.Most of the methods which we have presented in Sec.II, i.e., HF, GA, and NGA, share a common feature that they allows one to take the limit of N → ∞ numerically, and the results are free of finite-size (and statistical) errors.Same applies to SGA (see Ref. [58]) and CPA described in Appendix A. The case of GWF is different, since VMC simulations produced considerable statistical errorbars (a triple standard deviation is marked for each datapoint) and may be biased due to possible systematic errors following from a limited system size of N = 200 (N x = N y = 10).
Despite the limited accuracy of VMC simulation results, they typically lie between the GA and CPA values (up to the errorbars), allowing to regard the last to methods as providing approximate lower (GA) and upper (CPA) bounds to the value of U c .However, it must be noticed that the 'exact' numerical value of U (QMC) c = 3.86 t 0 of Ref. [18] (available only for the isotropic case, t x = t y = t 0 ) significantly exceeds U (CPA) c = 3.49 t 0 , and therefore the CPA results cannot be considered as upper bound to U c in a rigorous manner.
L 0 t a r 9 a 7 9 T G P Z q z 0 z z E s w P r + A 7 B V j H 4 = < / l a t e x i t > (t

SM (CSM) MI
< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 y 5 Q 0  1) and ( 2), with ty tx corresponding a gapless single-particle spectrum, see Fig.  12); thin dash-dotted line represents a power-law fit given by Eq. ( 23) with thin solid lines bounding the statistical uncertainty (yellow area).Quantum Monte Carlo value for the isotropic case, Uc/t0 = 3.86 [18], is also mark (full circle).Bottom: A zoom in, with trajectories following from the SSH model for strained graphene When searching for a computationally-inexpensive technique providing a safe upper bound to U c , one should rather refer to Neél-state Gutzwiller Approximation (NGA).
The relation between SGA and the above-mentioned methods is more complex, since more variational parameters defining the single-particle state |ψ 0 are optimized.In brief, the SGA ground-state energy lower or equal than the obtained from GA, leading to U . However, the mutual relation between SGA and VMC results cannot be determined a priori, as the former provides better optimization of |ψ 0 , whereas the latter put more emphasize on accurate calculation of averages in Eq. ( 13).Looking at the results presented in Fig. 5, we may conclude that U , finding SGA as slightly less accurate, but promising (due to much lower computational costs) counterpart to VMC.
The VMC results concerning U c can be rationalized within a power law, with least-square fitted parameters, as follows (Here, a single standard deviation is given for each parameter.)The line given by Eq. ( 23) [dashed-dotted], surrounded by the area [yellow] marking the statistical uncertainty, is further regarded as a border between semimetallic (SM) and Mott-insulating (MI) regions in the phase diagram.The former is further divided by marking the correlated-semimetal range (CSM), an appearance of which can be attributed to the fact that the HF approximation no longer produces a correct paramagnetic solution (m = 0).Such a computation-oriented notion cannot be regarded as a thermodynamic phase per se; however, prominent effects of electron correlations, i.e., the band narrowing and the reduction of double occupancies are gradually amplified when the interaction is increased.These effects are further discussed in next subsection, where we describe the behavior of measurable quantities when passing the CSM range and approaching the metal-insulator boundary.Three of the methods (HF, GA, and GWF) indicate U c → 0 for t y /t x → 0, coinciding with the exact solution for the Hubbard chain [3,4], giving an insulating phase at arbitrarily small U > 0. In contrast, CPA and NGA produces U c > 0 in such a limit, showing that these are inapplicable in the limit of weakly-coupled chains, despite producing a reasonable results in the isotropic case.(In particular, when comparing to the value of U c /t 0 = 10 given by DMFT [67].) Two striking features of the data shown in Fig. 5 are that most of the VMC datapoints do not match the GA line within the errorbars, but -on the other hand -the points for t y /t x < 0.8 match the CPA results surprisingly close.The above may indicate a role of finite-size effects in VMC simulations (notice that both GA and CPA solutions correspond to the N → ∞ limit).By manipulating the system sizes used for HF and GA calculations we found that shrinking to N x = N y = 10 usually produces U c /t x enlarged by 0.1 (HF) or 0.2 (GA) comparing to the large-system limit.Therefore, one could roughly estimate U (GWF) c /t x to be reduces by 0.2−0.3 when enlarging the system for t y ≈ t x .This quantity is comparable but smaller than the deviation from the 'exact' QMC result of Ref. [18], namely U suggesting that, in search for more accurate VMC results, one should first include additional variational parameters (such as Jastrow factors [68,69]), while the role of system size is rather secondary.
Also in Fig. 5 (bottom panel) we depict the trajectories followed by a sheet of graphene subjected to a strain in armchair direction (ε y > 0) and allowed to relax in the perpendicular (i.e., zigzag) direction.The hopping matrix elements in the < l a t e x i t s h a 1 _ b a s e 6 4 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 p 9 9 Figure 6: The evolution of critical Hubbard interaction with armchair strain strain (ty tx), approximated by Eq. ( 26) [solid lines] for two values of U (0) c (specified on the plot) adjusted to match the zero-strain results obtained from HF and GWF methods.Remaining lines are same as in Fig. 5.
Hamiltonian (1) are parametrized according to [48,70,71] where β = 2 (red solid line; open symbols) or β = 3 (blue solid line; closed symbols) is the dimensionless electronphonon coupling parameter.The bond-length variations (δd ij ) are adjusted to minimize the ground-state energy for an auxiliary Su-Schrieffer-Heeger model.(For more details, see Appendix B.) The effective Hubbard repulsion is approximated as with the coefficients U = 3.63 t 0 , and V 01 = 2.03 t 0 taken from Ref. [23], and . . .j(i) denoting the average over three nearest neighbors j of the site i. (Due to our suppositions on the symmetry, Eq. ( 25) produces same value for all sites.)Depending in the electron-phonon coupling β, we find the strain of ε y = 0.25 (for β = 2) or ε y = 0.20 (for β = 3) is necessary to approach U (HF) c , being a conventional border of the CSM range.This values are comparable with to the maximal strain of ≈ 0.20 reported in experiments.The phase of Mott Insulator seems inaccessible by applying mechanical strains to graphene, although modification of the equilibrium U/t 0 ratio due to substrate effect may possibly enhance the interaction effects.Below, we discuss the effects of electron correlations which should be visible also in CSM (or even SM) phase.

B. Effects of strain on measurable quantities
Earlier in this paper, we point out that a model assuming linear density of states, Eq. ( 11), parametrized by the Fermi velocity a zero energy, gives the critical Hubbard interaction U (X) c that differs by only 5% from the values obtained using < l a t e x i t s h a 1 _ b a s e 6 4 = " l q H z r t q T o D V i 7 .(For the numerical values, see Table I.) the actual density of states in the absence of strain, for the two methods, i.e., X = HF and X = GA.It is reasonable to expect, that for strains introducing anisotropy of the Fermi velocity [72,73] the value of U c will be affected predominantly via a change of the cut-off energy Λ, related to the Fermi velocity.For t y t x and strains limited to experimentally-accessible values, one can set Λ ∝ √ t x t y , leading to where U for t y = t x , we find (in Fig. 6) that the evolution of U c /t x with increasing strain is approximated by Eq. ( 26) quite well for both HF and GWF methods (similar agreement is observed for GA results, omitted in Fig. 6), but not for CPA, which predicts much weaker effects of strain.
Our results (in particular, a systematic shrinking of the SM phase, as well as the CSM range, with increasing strain) suggest that some measurable signatures of electron correlations should be visible in strained system also for U < U c .These expectation is further supported with the data presented in Fig. 7, where we display the average kinetic energy per site (quantifying the band narrowing) and the average double occupancy as functions of U .This time, the GWF results obtained fro VMC simulations do not differer significantly from GA results (see datapoints and black solid lines, respectively), while HF (dashed lines) predicts qualitatively different behavior, particularly for T displayed versus U in Figs.7(a), 7(b), and 7(c), but the differences between HF and Gutzwillerbased techniques are also apparent for n i↑ n i↓ (see remaining panels in Figs. 7).
For δ t = 0.5, see Figs. 7(c) and 7(f), corresponding to the strain of ε y = 0.22 for β = 3, and in the interval of U/t x = 1.5 ÷ 2 being relevant for graphene, the values of T obtained from GWF or GA are reduced by more than 20% comparing the δ t = 0 situation; see Fig. 7(a).(Notice that the above-mention reduction includes the change of t x , used as an energy unit in Fig. 7; the details are given in Appendix B.) Also in the intermediate case, δ t = 0.25, a 10% reduction is noticed, see Fig. 7(b).For n i↑ n i↓ , the effect of strain is less pronounced, but we still have an approximately 10% reduction for the δ t = 0.5 case [see Fig. 7(f)] compared to the δ t = 0 case [Fig.7(d)], following from both GWF or GA methods for the interval of U/t x = 1.5 ÷ 2.

IV. CONCLUDING REMARKS
We have investigated the mutual effect of electron-electron interaction, modeled by a Hubbard term in the secondquantized Hamiltonian, and geometric strains applied to a half-filled honeycomb lattice, quantified (at a first step) via two arbitrary values of the nearest-neighbor hopping integrals: one for bonds inside zigzag lines parallel to a selected direction, and the other for remaining bonds.Related problems were widely studied in the existing literature [24,66,72,73]; therefore, our attention has focussed on the case when strain is applied in the armchair direction (i.e., the hopping integrals connecting different zigzag lines are suppressed).In such a case, energy spectrum for noninteracting system remain gapless for arbitrary high strains, since the Dirac cones do not merge.In turn, the semimetal-insulator transition may occur only due to interactions.Also, the system gradually evolves, with increasing strain, towards a collection of weakly-coupled Hubbard chains, allowing to expect a considerable reduction of the critical Hubbard repulsion.
Several computational methods are compared, finding that the Hartree-Fock (HF) approximation, Gutzwiller Approximation (GA), and Gutzwiller Wave Function (GWF) treated within Variational Monte Carlo simulations, all predict qualitatively-similar shrinking of the semimetallic phase with increasing strain.Two remaining methods, the Coherent Potential Approximation (CPA) and so-called Neél-state GA, produce slightly different shapes of the semimetal-insulator boundary, but in the CPA case the results are numerically close to these obtained from GWF (provided that the strain is weak or moderate).
Phase diagram for the parametrized model is supplemented with calculations of trajectories, followed by monolayer graphene strained in armchair direction and allowed to relax along the perpendicular (i.e., zigzag) direction.These calculations were performed employing modified Su-Schrieffer-Heeger Hamiltonian, including the harmonic terms for bonds and angles (with the parameters fixed to reproduce elastics properties for in-plane small deformations), and the term describing coupling between electrons and the lattice, with dimensionless parameter varied between the possible values.
Albeit the critical and the actual Hubbard repulsion approach each other with increasing strain, we find that the semimetal-Mott insulator transition point cannot be achieved in monolayer graphene subjected to non-destructive deformations.Instead, one can observe the effects of electron correlations, such as the bandwidth renormalization or the reduction of double occupancies, which are well-pronounced (and affected by an applied strain) much before the transition.Probably, the above-mentioned effects will also be relevant in novel two-dimensional materials, predicted to sustain geometric deformations up to about 30% without structural demages [74,75].When looking for experimental realization of the Mott insulator on a honeycomb lattice, one should rather focus on artificial graphene-like systems [36][37][38][39].
What is more, we have shown that a basic version of the Gutzwiller Approximation already captures crucial correlation effects in (strained) graphene, allowing one to expect that proper generalizations of this method, such as the diagrammatic expansion [8,60,62], may lead to the development of versatile computational tools for studying graphene and related Dirac systems, providing low-(or moderate-) costs counterparts for the Quantum Monte Carlo methods.
L 0 t a r 9 a 7 9 T G P Z q z 0 z z E s w P r + A 7 B V j H 4 = < / l a t e x i t > (t L 0 t a r 9 a 7 9 T G P Z q z 0 z z E s w P r + A 7 B V j H 4 = < / l a t e x i t > (t < l a t e x i t s h a _ b a s e = " / R J l o V h S u S B P + z D l j u e w D = " ty/t0, = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " N P j w < l a t e x i t s h a 1 _ b a s e 6 4 = " I 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 + t 6 J B j a x 4 q D Here, t 0 = 2.7 eV is the equilibrium nearest-neighbor hopping integral for π electrons in monolayer graphene, is the dimensionless parameter quantifying the electron-phonon coupling (later, we perform main calculations for β = 2 and β = 3), and δd ij = d ij − d 0 is the bond-length change, calculated with respect to the equilibrium length of d 0 = a/ √ 3 = 0.142 nm.The second and third term in H SSH represent potential energy describing the covalent bonds [49], with (j) denoting the angles having a common vertex at a given lattice site j (see Fig. 1).
The parameters K d = 40.67eV/Å 2 , K θ = 5.46 eV/rad 2 , and θ 0 = π/3, are adjusted to restore the in-plane elastic coefficients of bulk graphene in the β = 0 case [49].For β = 0, a correction to the effective potential energy per bond (2/3) ∂ 2 0 /∂d 2 ij {dij =d0} = 0, provided that the equilibrium bond length is unaffected.To guarantee the last condition, we introduce the effective (β-dependent) equilibrium length where we have substituted the equilibrium kinetic energy per site 0 = −1.57460t 0 .(Notice that a standard constrain to the SSH model, ij d ij = const., is irrelevant when studying graphene with a global strain.) Next step is the optimization of the ground-state energy, with respect to in-plane atomic positions {R j }.In this paper, we limit the discussion to atomic arrangements preserving the bipartite structure structure of the lattice and the two mirror symmetries.For a fixed strain in the selected direction, ε > = max(ε x , ε y ), there are two parameters left to be optimized: the elongation in the perpendicular direction, ε < = min(ε x , ε y ), and the length d y of the bonds parallel to y-axis.The length of the remaining bonds, belonging to zigzag line, is given by Numerical values of the hopping matrix elements, following from the optimization procedure for β = 2, 3, and the two directions of strain are displayed in Fig. 8.In principle, the optimization scheme similar to the above can be constructed also for the Hamiltonian containing both lattice degrees of freedom and the Hubbard repulsion, H SSH + U D. For instance, if U < U c and the paramagnetic phase can be assumed, one can refer directly to Eq. ( 19), and setup a self-consistent procedure sharing main features of the ED-ABI method [51], but not limited to small systems.However, the term ∝ U 2 in Eq. (19), introducing the coupling between lattice and electron correlations, is preceded by a small prefactor and thus the corrections to atomic arrangements due to U > 0 are insignificant (notice that the Hubbard repulsion for a monolayer in equilibrium is U eff ≈ 1.6 t 0 [23]).For these reasons, we decided not to pursue this direction here, limiting our presentation to the data obtained directly via minimizing E (SSH) G ({R j }) in Eq. (B3), leaving the geometric parameters and the correlation effects decoupled.
7 k R G S U o Y i c W g I F W M Y j b r y v p S o y A 1 t o Q L L e 2 G T A y 5 5 o L s R f K 2 u r d a d J 0 0 y i X v u l S u 3 R Q r 9 9 k R c n A O F 3 A F H t x C B R 6 h C n U Q g P A G n / D l P D u v z r v z s Y h u O N m f M 1 i C 8 / 0 H x m y I n A = = < / l a t e x i t > j < l a t e x i t s h a 1 _ b a s e 6 4 = " / D x U O I Y O 2 j w 7 d 6 R n w w 5 p m u o u k 1 M = " > A A A B 3 n i c b V B N S 0 J B F L 3 P v s y + r J Z t h i R o J e 9 Z U E u p T U u F / A A V m z d e d X L e B z P 3 B S J u 2 0 S 0 K e g n 9 R f 6 N 4 3 6 N m o H B g 7 n n O H e c / 1 Y S U O u + + t k N j a 3 t n e y u 7 m 9 / Y P D o / z x S d 1 E i R Z Y E 5 G K d N P n B p U M s U a S F D Z j j T z w F T b 8 0 f 3 M b 7 y g N j I K H 2 k c Y y f g g 1 D 2 p e B k p e p z N 1 9 w i + 4 c b J 1 4 K S l A i k o 3 / 9 P u R S I J M C S h u D E t z 4 2 p M + G a p F A 4 z b U T g z E X I z 7 A y X y 9 K b u wU o / 1 I 2 1 f S G y u L u V 4 Y M w 4 8 G 0 y 4 D Q 0 q 9 5 M / M 9 r J d S / 7 U x k G C e E o V g M 6 i e K U c R m X V l P a h S k x p Z w o a X d k I k h 1 1 y Q v U j O V v d W i 6 6 T e q n o X R V L 1 e t C + S 4 9 Q h b O 4 B w u w Y M b K M M D V K A G A h D e 4 B O + n C f n 1 X l 3 P h b R j J P + O Y U l O N 9 / x + a I n Q = = < / l a t e x i t >d ij < l a t e x i t s h a 1 _ b a s e 6 4 = " J 6 m j y 8 X K M j V t s Q N + 0 k n r 9 8 3 C N o k = " > A A A B 4 3 i c b V D L S g M x F L 3 j s 9 Z X 1 a W b Y B F c l Z k q 6 L L o x m U F + 4 C 2 l E z m T p s 2 k x m S j F C G f o E b E T c K / o + / 4 N + Y t r N p 6 4 H A 4 Z w T 7 j 3 X T w T X x n V / n Y 3 N r e 2 d 3 c J e c f / g 8 O i 4 d H L a 1 H G q G D Z Y L G L V 9 q l G w S U 2 D D c C 2 4 l C G v k C W / 4 = " P t I 0 A K / + O z 5 e T 8 w 7 I u p + b l 7 d 5 BA = " > A A A B 8 n i c b V D L S s N A F L 2 p r 1 p f 1 S 7 d B I t Q N y V R Q Z d F N y 4 r 2 A e 0 p U w m k 3 b s Z B J m b o Q Q + i N u R N w o + C H + g n / j t M 2 m r Q c G D ue c m b n n e r H g G h 3 n 1 y p s b G 5 t 7 x R 3 S 3 v 7 B 4 d H 5 e O T t o 4 S R V m L R i J S X Y 9 o J r h k L e Q o W D d W j I S e Y B 1 v c j / z O y 9 M a R 7 J J 0 x j N g j J S P K A U 4 J G G p Y r W T 9 k R C e K + U S O B J v W n i + G 5 a p T d + a w 1 4 m b k y r k a A 7 L P 3 0 / o k n I J F J B t O 6 5 T o y D j C j k 1 D x Z 6 i e a x Y R O y I h l 8 5 G n 9 r m R f D u I l D k S 7 b m 6 l C O h 1 m n o Figure 1: Top: Honeycomb lattice subjected to uniaxial strain in selected direction (see the coordinate system).Zoom-in visualizes the distance (bond length) dij between atoms i and j and in-plane angle with the vertex at site j ( (j)).Bottom: Hexagonal first Brillouin zone (FBZ) of the reciprocal lattice, with (dimensionless) basis vectors b1 = 2π/ √ 3 ( √ 3, −1) and b2 = 2π/ √ 3 (0, 2), and the symmetry points, K = (4π/3, 0) and K = (2π/3, 2π/ √ 3) coinciding with Dirac points in the absence of strain.Magnified area shows discretized FBZ for a finite system of N = 2NxNy atoms with periodic boundary conditions [see Eq. (6)].(The values of Nx = 4 and Ny = 5 are used for illustration only.)

Figure 2 :
Figure 2: Density of states for the Hamiltonian (1) with U = 0 displayed as a function of energy.Top: strain applied in the armchair direction (ty tx), bottom: strain applied in the zigzag direction (ty tx).The ratio t</t> [with t< = min (tx, ty) and t> = max (tx, ty)] is varied between the lines with the the steps of 0.2.A vertical offset is applied to each dataset except from the isotropic case (tx = ty).Inset shows the band gap, appearing for tx < 0.5 ty due to the Peierls transition.
Unlike for square lattice, for which one gets m = 0 for any U > 0 [5], on a honeycomb lattice the minimization gives m = 0 for U U (HF) c and m = 0 for U > U (HF) c o 6 n U l 7 x 0 g I 3 + D G G D e a + D f + g n 9 j g d k A n q T J y T m n u f d c P 5 b C o O v + O p m 1 9 Y 3 N r e x 2 b m d 3 b / 8 g f 3 h U N 1 G i G a + x S E a 6 6 V P D p V C 8 h g I l b 8 a a 0 9 C X v O E P b 6 d + 4 5 l r I y L 1 D B g r e 4 B O + H O G 8 O u / O x z y a c d I / x 7 A A 5 / s P C r S M N g = = < / l a t e x i t > D B g r e 4 B O + H O G 8 O u / O x z y a c d I / x 7 A A 5 / s P C r S M N g = = < / l a t e x i t > /t x < l a t e x i t s h a 1 _ b a s e 6 4 = " y S 1 D 5 D B g r e 4 B O + H O G 8 O u / O x z y a c d I / x 7 A A 5 / s P C r S M N g = = < / l a t e x i t > /t x < l a t e x i t s h a 1 _ b a s e 6 4 = " v + h P w y r h M Y 5 I 5 s v 5 9 M +

Figure 3 :
Figure 3: (a)-(d) Main: Energy difference between the antiferromagnetic Gutzwiller variational energy E (GWF) G (m) [see Eq. (13)] and the paramagnetic solution E (GWF) G (0) obtained from VMC simulations as a function the on-site Hubbard repulsion (U ).The parameter η is optimized for a fixed m = ∆/U (or m = 0); ∆ is varied between the lines from ∆/tx = 0.25 to ∆/tx = 1, with the steps of 0.25.The hopping anisotropy ty/tx is varied between the panels.Inset shows the value of U = U0(∆) at which ∆E (GWF) G changes sign for a given ∆.The extrapolation to ∆ → 0 yields the critical values of U (GWF) c

1 4 ( 1 − 4 U ( 1 − m 2 )
m 2 ), with the upper limit corresponding to the average double occupancy in the uncorrelated state |ψ 0 (m) ).The kinetic energy term can be estimated by referring to the Hatree-Fock energy E even for m being away from the minimum of E (HF) G 9 1 b L r p K m l c V 7 7 p S f a y W a x f 5 E Q p w C m d w C R 7 c Q A 0 e o A 4 N Y D C C N / i E L y d 0 X p 1 3 5 2 M e X X P y P y e w A O f 7 D z e Y i y w = < / l a t e x i t > t y = t x (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " n K e L L p P G p e 1 e 2 5 X H S r l 6 N j t C A Y 7 h B M 7 B h R u o w g P U o A 4 M U n i D T / i y B t a r 9 W 5 9 T K M r 1 u z P E c z B + v 4 D n z S M e Q = = < / l a t e x i t > t y = 0.5 t x t e x i t s h a 1 _ b a s e 6 4 = " 8 b 3 0 r 8 W 5 8 z K M F I / 9 z h B Z g f P 8 B k f a N d w = = < / l a t e x i t >t y = 0.25 t x < l a t e x i t s h a 1 _ b a s e 6 4 = " p Q o / o R H c R / X M t V d W / O U b T Y 5 i r T I = " > A A A B 5 H i c b V D L S g M x F L 3 x W e u r 6 t J N s A j i o s 5 I U Z c F E V 1 W s A 9 o y 5 B J M 2 1 o 5 k F y R y y l f + B G x I 2 C 3 + M v + D e m 7 Wz a e i B w O O e E e 8 / 1 E y U N O s 4 v W V l d W 9 / Y z G 3 l t 3 d 2 9 / Y L B 4 d 1 E 6 e a i x q P V a y b P j N

4 9 3 Figure 4 :
Figure 4: (a)-(d).The lower and the upper bounds to the critical Hubbard repulsion Uc for anisotropic honeycomb lattice estimated by comparing different versions of the Gutzwiller Approximation described in the text.The lower bound (U (GA) c ) coincides with the splitting of the Gutzwiller energy for paramagnetic state, E (GA) G (m = 0), given by Eq. (19) [red dashed line] and the variational energy E (GA) G , see Eq. (17), with the parameters (m, d) optimized numerically [blue solid line], both displayed as functions of U .The value of U (GA) c is obtained via the extrapolation with m → 0, similarly as for the VMC results in Fig. 3.The intersection of E (GA) G (m = 0) with E NGA G , see Eq. (20) [green dashed-dotted line] yields the upper bound (U (NGA) c ).The value of ty/tx ratio is varied between the panels.[For the numerical values of U (GA) c and U (NGA) c t e x i t s h a 1 _ b a s e 6 4 = " 5 3 Z z 4 C E z k r 8 Y 0 g T b H / 7 n + F 6 U G p E = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q D 8 b d E N S D S M C L x w j m A c m y z E 4 m y

H g 1 3 o 2 P"Figure 5 :
Figure 5: Top: Phase diagram for the Hubbard model on anisotropic honeycomb lattice, see Eqs. (1) and (2), with ty tx corresponding a gapless single-particle spectrum, see Fig. 2. [Here, t> = tx and t< = ty.]Lines depict the critical Hubbard repulsion estimated within the Hartree-Fock method [short dashed], Gutzwiller Approximation [thick solid], Statistically-consistent GA [long dashed-double dotted], Coherent Potential Approximation [dotted], and the Néelstate GA [long dashed-dotted].Datapoints with errorbars are obtained from VMC simulations for the Gutzwiller Wave Function, see Eq. (12); thin dash-dotted line represents a power-law fit given by Eq. (23) with thin solid lines bounding the statistical uncertainty (yellow area).Quantum Monte Carlo value for the isotropic case, Uc/t0 = 3.86[18], is also mark (full circle).Bottom: A zoom in, with trajectories following from the SSH model for strained graphene (see Appendix B) for β = 2 (red line/open symbols) and β = 3 (blue line/closed symbols).Different datapoints for each value of β correspond to the applied strain εy varied from εy = 0.05 to 0.25 with the steps of 0.05.(GA and VMC results are omitted for clarity.)Remaining Labels/colored areas: the semimetallic phase (SM) [blue] with the correlated-semimetal range (CSM) [light blue] and the Mott insulator (MI) [magenta].
Figure 5: Top: Phase diagram for the Hubbard model on anisotropic honeycomb lattice, see Eqs. (1) and (2), with ty tx corresponding a gapless single-particle spectrum, see Fig. 2. [Here, t> = tx and t< = ty.]Lines depict the critical Hubbard repulsion estimated within the Hartree-Fock method [short dashed], Gutzwiller Approximation [thick solid], Statistically-consistent GA [long dashed-double dotted], Coherent Potential Approximation [dotted], and the Néelstate GA [long dashed-dotted].Datapoints with errorbars are obtained from VMC simulations for the Gutzwiller Wave Function, see Eq. (12); thin dash-dotted line represents a power-law fit given by Eq. (23) with thin solid lines bounding the statistical uncertainty (yellow area).Quantum Monte Carlo value for the isotropic case, Uc/t0 = 3.86[18], is also mark (full circle).Bottom: A zoom in, with trajectories following from the SSH model for strained graphene (see Appendix B) for β = 2 (red line/open symbols) and β = 3 (blue line/closed symbols).Different datapoints for each value of β correspond to the applied strain εy varied from εy = 0.05 to 0.25 with the steps of 0.05.(GA and VMC results are omitted for clarity.)Remaining Labels/colored areas: the semimetallic phase (SM) [blue] with the correlated-semimetal range (CSM) [light blue] and the Mott insulator (MI) [magenta].
Figure 5: Top: Phase diagram for the Hubbard model on anisotropic honeycomb lattice, see Eqs. (1) and (2), with ty tx corresponding a gapless single-particle spectrum, see Fig. 2. [Here, t> = tx and t< = ty.]Lines depict the critical Hubbard repulsion estimated within the Hartree-Fock method [short dashed], Gutzwiller Approximation [thick solid], Statistically-consistent GA [long dashed-double dotted], Coherent Potential Approximation [dotted], and the Néelstate GA [long dashed-dotted].Datapoints with errorbars are obtained from VMC simulations for the Gutzwiller Wave Function, see Eq. (12); thin dash-dotted line represents a power-law fit given by Eq. (23) with thin solid lines bounding the statistical uncertainty (yellow area).Quantum Monte Carlo value for the isotropic case, Uc/t0 = 3.86[18], is also mark (full circle).Bottom: A zoom in, with trajectories following from the SSH model for strained graphene (see Appendix B) for β = 2 (red line/open symbols) and β = 3 (blue line/closed symbols).Different datapoints for each value of β correspond to the applied strain εy varied from εy = 0.05 to 0.25 with the steps of 0.05.(GA and VMC results are omitted for clarity.)Remaining Labels/colored areas: the semimetallic phase (SM) [blue] with the correlated-semimetal range (CSM) [light blue] and the Mott insulator (MI) [magenta].

4 5 w 7 3
n h q k U B l 3 3 1 9 n Y 3 N r e 2 S 3 s F f c P D o + O S y e n T a M y z b j P l F S 6 H V L D p U i 4 j w I l b 6 e a 0 z i

Figure 7 :
Figure 7: (a)-(c) Average kinetic energy per site and (d)-(f) average double occupancy displayed as functions of the Hubbard repulsion U for ty/tx = 1 (top), ty/tx = 0.75 (middle), and ty/tx = 0.5 (bottom).Thick dashes line marks the Hartree-Fock results, thick solid line represents the Gutzwiller Approximation.Datapoints depict the VMC results for GWF with m = 0 (red open symbols) and optimized m (blue closed symbols); thin lines are guide for the eye only.Shaded area marks the correlated semimetallic phase, bounded by U (HF) c and U (GWF) c

6 <
l a t e x i t s h a 1 _ b a s e 6 4 = " J U R e + W r 1 u u f O 5 h t j u g f I 7 6 E p p q I

h 7 7 3 < l a t e x i t s h a 1 _ 3 < l a t e x i t s h a 1 _ 2 Figure 8 :
Figure 8: The hopping-matrix elements [see Eqs.(2) and(24) in the main text] as functions of the strain εx > 0 (a), (b) or εy > 0 (c), (d).The electron-phonon coupling is fixed at β = 2 (dashed lines in all plots) or β = (solid lines).

Table I :
Critical values of the Hubbard repulsion U