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

Abstract

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.

1. Introduction

The Hubbard model, initially proposed to describe interaction-driven transition between conducting and insulating systems [1,2], needs to be carefully applied in low dimensions, where exact solutions (when available) [3,4] show substantially different ground-state properties than approximate solutions obtained using methods such as Hartree–Fock (HF) [2,5], GWF or Gutzwiller Approximation (GA) [6,7,8,9,10]. For this reason, computationally-expensive numerical techniques, such as Quantum Monte Carlo (QMC) [11] or a more recent tensor-network method [12,13,14], are usually employed for the Hubbard model in two dimensions, for which an exact solution is missing.

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\left(1\right)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_c^{\left(\mathrm{HF}\right)}=2.23t_0$ [19] for the same lattice.

Since the advent of graphene [20,21], 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_{\mathrm{eff}}\approx 1.6t_0$ was proposed [23], suggesting that large isotropic strain may drive this system from semimetallic towards the Mott-insulating phase [24,25] in analogy with high pressure changing properties of various bulk materials [26,27,28,29,30]. The 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 opened with the fabrication of twisted bilayer graphene [40,41].

A separate issue concerns the bandgap opening due to the spatial rearrangement of atoms in strained graphene (a so-called two-dimensional Peierls instability), which may turn the system into an insulator before the Mott transition occurs [42,43,44,45,46]. On the contrary, a 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 Figure 1) supposing that the bipartite structure of the lattice is preserved under strain. In turn, there are two different values of the nearest-neighbor 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 the 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 values of $U_{\mathrm{eff}}$ than for isotropic strain.

The remaining part of the paper is organized as follows. In Section 2, we briefly present approximate approaches to the effective Hubbard Hamiltonian (including HF, GWF, and GA). Moreover, in Section 2, we show some original data, illustrating how these approaches work for anisotropic honeycomb lattice. In Section 3, we discuss our numerical results concerning the phase diagram of the effective Hubbard model with arbitrary parameters ($t_x⩾t_y$ and $U_{\mathrm{eff}}$), the evolution of the model parameters in graphene subjected to uniaxial strain, and approximate formula relating the reduction in $U_c$ to strain-induced anisotropy of the Fermi velocity. The effects of electron correlations on selected measurable quantities are also presented in Section 3. The concluding remarks are given in Section 4.

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.

2. Model and Methods

2.1. The Anisotropic Hubbard Model

Our analysis of electron correlations on the anisotropic honeycomb lattice starts from the Hamiltonian

$$H=\sum _{\langle ij\rangle ,s}{t}_{ij}\left(c_{i,s}^{†}{c}_{j,s}+\mathrm{H}.\mathrm{c}.\right)+U\sum _j{n}_{j\uparrow }{n}_{j\downarrow },$$

(1)

with the first sum running over pairs of nearest-neighbors $\langle ij\rangle$ and spin up/down orientations ($s=\uparrow ,\downarrow$), and the hopping-matrix elements are given by

$$t_{ij}=\left\{\begin{array}{cc}-t_x\hfill & \mathrm{if} i,j \mathrm{belongs\; to\; same\; zigzag\; line},\hfill \\ -t_y\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(2)

(Without loss of generality, we suppose the coordinate system is oriented as depicted in Figure 1) The remaining symbols in Equation (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, $\overline{n}=\langle n_{i\uparrow }+n_{i\downarrow }\rangle =1$.

In principle, the ground-state properties of the model defined by Equations (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. However, first, we briefly present approximate approaches capable to distinguish whether ground state of the Hamiltonian (1) is semimetallic or insulating.

2.2. 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 Figure 2, there is no Fermi surface that could produce magnetic instability also for small $U>0$.

Within the Hartree–Fock approximation, the interaction part in the Hamiltonian (1) is replaced by

$$U\widehat{D}\stackrel{\mathrm{HF}}{=}U\sum _i\left(\langle n_{i\uparrow }\rangle n_{i\downarrow }+n_{i\uparrow }\langle n_{i\downarrow }\rangle -\langle n_{i\uparrow }\rangle \langle n_{i\downarrow }\rangle \right),$$

(3)

where we have introduced the operator $\widehat{D}={\sum }_j{n}_{j\uparrow }{n}_{j\downarrow }$ measuring the number of double occupancies. We further impose the antiferromagnetic order,

$$\langle n_{i\uparrow }\rangle =\frac{\overline{n}+{\lambda }_{i}m}{2},\langle n_{i\downarrow }\rangle =\frac{\overline{n}-{\lambda }_{i}m}{2},$$

(4)

where ${\lambda }_i=1$ if i belongs to one sublattice (A), or ${\lambda }_i=-1$ if i belongs to the other sublattice (B), and m is the magnetization ($\left|m\right|⩽\overline{n}$), and the half filling ($\overline{n}=1$). The above yields the HF ground-state energy per site

$$\begin{array}{cc}\hfill \frac{E_G^{\left(\mathrm{HF}\right)}}{N}& =-\frac{2}{N}\sum _{\mathbf{k}}\sqrt{E_{\mathbf{k}}^2+{\left(\frac{Um}{2}\right)}^2}+\frac{U(1+m^2)}{4},\hfill \end{array}$$

(5)

where the factor 2 accounts for $s=\uparrow ,\downarrow$, and the summation runs over quasimomenta $\mathbf{k}\equiv (k_x,k_y)$ in the first Brillouin zone, namely

$$\begin{array}{cc}\hfill k_x& =\frac{2\pi }{N_x}{n}_x,k_y=\frac{4\pi }{\sqrt{3}}\left(\frac{n_y}{N_y}-\frac{n_x}{2N_x}\right),\hfill \\ \hfill n_x& =0,1,\cdots ,N_x-1,n_y=0,1,\cdots ,N_y-1,\hfill \end{array}$$

(6)

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 a sufficiently large number of points in the momentum space, say $N_x,N_y\gtrsim {10}^3$, one can usually works with a square inverse lattice (omitting the term $\propto n_x$ in the expression for $k_y$) [53]; nevertheless, the discretization of $(k_x,k_y)$ as given in Equation (6) becomes crucial when discussing the finite-size effects for small N. The single-particle energies for anisotropic honeycomb lattice are given by

$$E_{\mathbf{k}}=t_x\sqrt{a_{\mathbf{k}}^2+b_{\mathbf{k}}^2},$$

(7)

with

$$\begin{array}{cc}\hfill a_{\mathbf{k}}=& -cos\left(\frac{k_x}{2}+\frac{\sqrt{3}{k}_y}{2}\right)-cos\left(\frac{k_x}{2}-\frac{\sqrt{3}{k}_y}{2}\right)-\frac{t_y}{t_x},\hfill \\ \hfill b_{\mathbf{k}}=& sin\left(\frac{k_x}{2}+\frac{\sqrt{3}{k}_y}{2}\right)-sin\left(\frac{k_x}{2}-\frac{\sqrt{3}{k}_y}{2}\right).\hfill \end{array}$$

(8)

Next, the density of states is defined as

$$\rho \left(E\right)=2N^{-1}\sum _{\mathbf{k}}\left[\delta \left(E-E_{\mathbf{k}}\right)+\delta \left(E+E_{\mathbf{k}}\right)\right],$$

(9)

with two parts corresponding to the conduction ($E>0$) and valence ($E<0$) band, and satisfying the normalization conditions: ${\int }_{-\infty }^{0}dE\rho \left(E\right)={\int }_0^{\infty }dE\rho \left(E\right)=1$. If $m\ne 0$, the minimization of $E_G^{\left(\mathrm{HF}\right)}$ given by Equation (5) brought us to

$$1={\int }_{E<0}dE\rho \left(E\right)\frac{U/2}{\sqrt{E^2+{(Um/2)}^2}}.$$

(10)

In case the solution of Equation (10) does not exist, the minimum of $E_G^{\left(\mathrm{HF}\right)}$ corresponds to $m=0$.

Unlike for square lattice, for which one obtains $m\ne 0$ for any $U>0$ [5], on a honeycomb lattice the minimization gives $m=0$ for $U⩽U_c^{\left(\mathrm{HF}\right)}$ and $m\ne 0$ for $U>U_c^{\left(\mathrm{HF}\right)}$ [15]. This can be easily understood in 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

$$\rho \left(E\right)\approx {\rho }_{\mathsf{\Lambda }}\left(E\right)=\left\{\begin{array}{cc}\frac{2}{\mathsf{\Lambda }{}^2}\left|E\right|\hfill & \mathrm{for}\left|E\right|⩽\mathsf{\Lambda },\hfill \\ 0\hfill & \mathrm{for}\left|E\right|>\mathsf{\Lambda },\hfill \end{array}\right.$$

(11)

with a cut-off energy of $\mathsf{\Lambda }={\left(\sqrt{3}\pi \right)}^{1/2}{t}_0\simeq 2.33268t_0$. The above is equivalent, for $\left|E\right|⩽\mathsf{\Lambda }$, to ${\rho }_{\mathsf{\Lambda }}\left(E\right)=2\mathcal{A}\left|E\right|/\left[N\pi {(\hslash v_F)}^2\right]$, with $\mathcal{A}$ being the system area, $v_F=\frac{1}{2}\sqrt{3}a{t}_0/\hslash$ the Fermi velocity, and a the lattice parameter [54]. It is straightforward to show that $m\ne 0$ appears above $U_c^{\left(\mathsf{\Lambda }\right)}=\mathsf{\Lambda }$, being not far from the value reported in Ref. [19].

The values of $U_c^{\left(\mathrm{HF}\right)}$ following from numerical minimization of $E_G^{\left(\mathrm{HF}\right)}$ given by Equation (5) for the actual density of states are presented in Section 3.

2.3. Gutzwiller Wavefunction

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 GWF are surprisingly close to those obtained within large-scale computer simulations for a honeycomb lattice. Investigating the variational wavefunction

$$\left|{\mathsf{\Psi }}_{\mathrm{GWF}}\right.〉=e^{-\eta \widehat{D}}\left|{\psi }_0\left(m\right)\right.〉,$$

(12)

where $\left|{\psi }_0\left(m\right)\right.〉$ denotes a Slater determinant corresponding to a given magnetization in Equation (5) and $\eta$ is another variational parameter (quantifying the role of electron correlations), one needs to minimize the ground-state energy

$$E_G^{\left(\mathrm{GWF}\right)}=\frac{\left.〈{\psi }_0\left(m\right)\right|e^{-\eta \widehat{D}}H{e}^{-\eta \widehat{D}}\left|{\psi }_0\left(m\right)\right.〉}{\left.〈{\psi }_0\left(m\right)\right|e^{-2\eta \widehat{D}}\left|{\psi }_0\left(m\right)\right.〉},$$

(13)

with respect to $\eta$ and m. In many cases, the system may prefer to reduce m (even to $m=0$) and increase $\eta$, allowing to expect that, in general, $U_c^{\left(\mathrm{GWF}\right)}⩾U_c^{\left(\mathrm{HF}\right)}$.

Several approximated techniques for calculating the averages in Equation (13) were developed [7,8,9,10,15]. Here, we apply Variational Monte Carlo (VMC), described in detail in Ref. [10]. To determine the value of $U_c^{\left(\mathrm{GWF}\right)}$, we have directly followed the procedure proposed by Martelo et al. [15]. For a fixed value of the gap ($\Delta \equiv Um$), the energy difference $\Delta E_G^{\left(\mathrm{GWF}\right)}=E_G^{\left(\mathrm{GWF}\right)}\left(m\right)-E_G^{\left(\mathrm{GWF}\right)}\left(0\right)$, where the parameter $\eta$ is optimized independently for $m=0$ and $m\ne 0$, changes sign at some $U=U_0(\Delta )$. Numerical extrapolation of $U_0(\Delta )$ with $\Delta \to 0$ allows one to determine the critical value of $U_c^{\left(\mathrm{GWF}\right)}$. 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 Figure 3. For more details of the simulation, see Ref. [55].

2.4. Gutzwiller Approximation and Its Variants

To efficiently study the effects of electron correlations present in $|{\mathsf{\Psi }}_{\mathrm{GWF}}\rangle$, see Equation (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 $|{\psi }_0\left(m\right)\rangle$. Within GA, which is exact in the infinite dimension limit, the correlation functions ${\langle c_{is}^{†}{c}_{js}\rangle }_{\mathrm{GWF}}=\left.〈{\mathsf{\Psi }}_{\mathrm{GWF}}\right|c_{is}^{†}{c}_{js}\left|{\mathsf{\Psi }}_{\mathrm{GWF}}\right.〉/\left.〈{\mathsf{\Psi }}_{\mathrm{GWF}}\right|{\mathsf{\Psi }}_{\mathrm{GWF}}〉$ are approximated by

$${\langle c_{is}^{†}{c}_{js}\rangle }_{\mathrm{GWF}}\stackrel{\mathrm{GA}}{=}q(\left\{{\rho }_{l{l}^{\prime }{s}^{\prime }}^{\left(0\right)}\right\};\left\{d_l\right\}){\langle c_{is}^{†}{c}_{js}\rangle }_0$$

(14)

where the band-narrowing factor $q(\left\{{\rho }_{l{l}^{\prime }{s}^{\prime }}^{\left(0\right)}\right\};\left\{d_l\right\})$ depends only on the single-particle density matrix elements ${\rho }_{l{l}^{\prime }{s}^{\prime }}^{\left(0\right)}={\langle c_{l{s}^{\prime }}^{†}{c}_{l^{\prime }{s}^{\prime }}\rangle }_0$ with $l,l^{\prime }=i,j$ and $s^{\prime }=\uparrow ,\downarrow$ (here, ${\langle \cdots \rangle }_0$ is the expectation value over the uncorrelated state; i.e., a single Slater determinant such as $|{\psi }_0\left(m\right)\rangle$), and $d_l={\langle n_{l\uparrow }{n}_{l\downarrow }\rangle }_{\mathrm{GWF}}$ ($l=i,j$) being the average double occupancies. The $\left\{d_i\right\}$ variables are further regarded as variational parameters to be determined by minimizing the Gutzwiller energy functional,

$$E_G^{\left(\mathrm{GA}\right)}=2\sum _{\langle ij\rangle ,s}q(\left\{{\rho }_{l{l}^{\prime }{s}^{\prime }}^{\left(0\right)}\right\};\left\{d_l\right\})t_{ij}{\rho }_{ijs}^{\left(0\right)}+U\sum _j{d}_j.$$

(15)

Several forms of the band-narrowing factor $q(\left\{{\rho }_{l{l}^{\prime }s}^{\left(0\right)}\right\};\left\{d_l\right\})$, being equivalent in the infinite dimension limit but producing slightly different results when applied to the system of a finite dimensionality, are used in the literature [8,56,57,58,59,60,61,62]. For the diagonal elements ${\rho }_{iis}^{\left(0\right)}={\langle n_{is}\rangle }_0$ parametrized as in Equation (4) with $\overline{n}=1$, one can impose $d_i\equiv d$ for all sites and rewrite the expression given in Ref. [61] as

$$\begin{array}{cc}\hfill q(\left\{{\rho }_{l{l}^{\prime }s}^{\left(0\right)}\right\};\left\{d_l\right\})& \equiv q(m,d)=\frac{4d}{1-m^2}\left[1-2d+\sqrt{{(1-2d)}^2-m^2}\right].\hfill \end{array}$$

(16)

The variable d is bounded as $0⩽d⩽\frac{1}{4}(1-m^2)$, with the upper limit corresponding to the average double occupancy in the uncorrelated state $|{\psi }_0\left(m\right)\rangle$). The kinetic energy term can be estimated by referring to the Hartree–Fock energy $E_G^{\left(\mathrm{HF}\right)}$, see Equation (5), as $2{\sum }_{\langle ij\rangle ,s}{t}_{ij}{\rho }_{ij}^{\left(0\right)}=E_G^{\left(\mathrm{HF}\right)}-\frac{N}{4}U(1-m^2)$ even for m being away from the minimum of $E_G^{\left(\mathrm{HF}\right)}$. This brought us to

$$\begin{array}{cc}\hfill \frac{E_G^{\left(\mathrm{GA}\right)}}{N}& =q(m,d)\times \left[-\frac{2}{N}\sum _{\mathbf{k}}\sqrt{E_{\mathbf{k}}^2+{\left(\frac{Um}{2}\right)}^2}+\frac{U{m}^2}{2}\right]+Ud.\hfill \end{array}$$

(17)

Numerical minimization of $E_G^{\left(\mathrm{GA}\right)}$, with respect to $(m,d)$, truncates the optimization of both the density matrix $\left\{{\rho }_{i,j}^{\left(0\right)}\right\}$ and the parameters $\left\{d_i\right\}$. For the linear density of states ${\rho }_{\mathsf{\Lambda }}\left(E\right)$, see Equation (11), one can easily find closed-from expression for $E_G^{\left(\mathrm{GA}\right)}$; the minima corresponding to $m\ne 0$ appear for $U>U_c^{(\mathsf{\Lambda },\mathrm{GA})}=1.270\mathsf{\Lambda }=2.963t_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 $|{\psi }_0\rangle$ not only via the magnetization m, as in the above, but via all independent parameters of the density matrix $\left\{{\rho }_{i,j}^{\left(0\right)}\right\}$. In particular, the auxiliary single-particle Hamiltonian determining $\left\{{\rho }_{i,j}^{\left(0\right)}\right\}$ contains the renormalized hopping integrals (${\tilde{t}}_x$ and ${\tilde{t}}_y$) 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 detail in Ref. [58].

Both (S)GA and GWF methods can be regarded as improvements to the 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 Equation (17), with a complementary variational wavefunction,

$$|{\mathsf{\Psi }}_B\rangle =e^{-\kappa \widehat{T}}|{\mathsf{\Psi }}_{U\to \infty }\rangle ,$$

(18)

where $\kappa$ is a variational parameter, $\widehat{T}={\sum }_{\langle ij\rangle ,s}{t}_{ij}(c_{is}^{†}{c}_{js}+\mathrm{H.c.})$ is the kinetic-energy part of the Hamiltonian (1), and $|{\mathsf{\Psi }}_{U\to \infty }\rangle$ is the ground state for $U\to \infty$. The critical value $U_c$ is than that estimated by finding a crossing point of $E_G^{\left(\mathrm{GA}\right)}$, Equation (17), with a fixed $m=0$ and optimized d, and the variational energy $E_B$ corresponding to $|{\mathsf{\Psi }}_B\rangle$, Equation (18), with optimized $\kappa$.

For $m=0$, the factor $q(m,d)$, Equation (16), reduces to a quadratic function of d and the functional $E_G^{\left(\mathrm{GA}\right)}$, Equation (17), reaches the minimum at $d=\frac{1}{4}\max \left[0,1-U/\left(8\right|{ϵ}_0\left|\right)\right]$, leading to a form originally derived by Gutzwiller [63,64]

$$E_G^{\left(\mathrm{GA}\right)}(m=0)=\left\{\begin{array}{cc}{\displaystyle {ϵ}_0+\frac{U}{4}-\frac{U^2}{64|{ϵ}_0|}}\hfill & \mathrm{for}U⩽8|{ϵ}_0|,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(19)

The symbol ${ϵ}_0$ is the kinetic energy per site for $U=0$, namely ${ϵ}_0={\int }_{E<0}dE\rho \left(E\right)E$ [for the definition of $\rho \left(E\right)$, see Equation (9)], taking the numerical value of ${ϵ}_0/t_x=-1.57460$ for $t_y/t_x=1$, ${ϵ}_0/t_x=-1.45540$ for $t_y/t_x=0.75$, ${ϵ}_0/t_x=-1.36218$ for $t_y/t_x=0.5$, or ${ϵ}_0/t_x=-1.29891$ for $t_y/t_x=0.25$.

In the limit of infinite dimensions, the variational energy $E_B$ associated with the state $|{\mathsf{\Psi }}_B\rangle$ can be evaluated exactly, since the ground state $|{\mathsf{\Psi }}_{U\to \infty }\rangle$ is the Néel antiferromagnet [65]. The variational energy reads

$$\frac{E_B}{N}={ϵ}_{\mathrm{kin}}+\frac{U(1-m^2)}{4}\equiv \frac{E_G^{\left(\mathrm{NGA}\right)}}{N}$$

(20)

where

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{kin}}& ={\int }_{E<0}dE\rho \left(E\right)E\left[-tanh\left(2\kappa E\right)\right],\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill m& ={\int }_{E<0}dE\frac{\rho \left(E\right)}{cosh\left(2\kappa E\right)},\hfill \end{array}$$

(22)

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 Equation (9) into Equations (21) and (22), and the subsequent minimization of $E_B\equiv E_G^{\left(\mathrm{NGA}\right)}$ with respect to $\kappa$. Selected numerical results, for $t_y⩽t_x$, are presented in Figure 4.

It is worth mentioning that (S)GA can be systematically improved by approaching the GWF solution, 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). A similar approach for $|{\mathsf{\Psi }}_B\rangle$ is difficult due to the necessity of determining the ground state of the Heisenberg model ($|{\mathsf{\Psi }}_{U\to \infty }\rangle$) as a first.

The selected numerical values of $U_c$, following on from the methods described in this section, are compared in

Table 1. Critical values of the Hubbard repulsion $U_c^{\left(\mathrm{GWF}\right)}$ obtained from VMC simulations (with standard deviations for the last digit specified in parentheses) compared with the upper ($U_c^{\left(\mathrm{GA}\right)}$) and the upper ($U_c^{\left(\mathrm{NGA}\right)}$) 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 $N_x=N_y=10$ for VMC simulations; the remaining results correspond to the limit of $N_x=N_y\to \infty$.

${\mathit{t}}_{\mathit{y}}/{\mathit{t}}_{\mathit{x}}$	${\mathit{U}}_{\mathit{c}}^{\left(\mathbf{GWF}\right)}/{\mathit{t}}_{\mathit{x}}$	${\mathit{U}}_{\mathit{c}}^{\left(\mathbf{GA}\right)}/{\mathit{t}}_{\mathit{x}}$	${\mathit{U}}_{\mathit{c}}^{\left(\mathbf{SGA}\right)}/{\mathit{t}}_{\mathit{x}}$	${\mathit{U}}_{\mathit{c}}^{\left(\mathbf{NGA}\right)}/{\mathit{t}}_{\mathit{x}}$
1.00	3.48(1)	2.804	3.12	5.281
0.75	2.91(1)	2.550	2.83	4.871
0.50	2.69(3)	2.241	2.47	4.508
0.25	2.24(1)	1.830	1.98	4.199

.

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.

3. Results and Discussion

3.1. Phase Diagram

Our central results are presented in Figure 5, where we display the phase diagram for the Hamiltonian (1) with the strain applied in the armchair direction ($t_y⩽t_x$). A single-particle spectrum is gapless in such a case (see also Figure 2) since the positions of Dirac cones do not merge [66]; therefore, the metal-insulator (if it occurs) must be driven by electron–electron interaction. Most of the methods which we have presented in Section 2, i.e., HF, GA, and NGA, share a common feature that they allow one to take the limit of $N\to \infty$ numerically, and the results are free of finite-size (and statistical) errors. The 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 us to regard the last two 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_c^{\left(\mathrm{QMC}\right)}=3.86t_0$ of Ref. [18] (available only for the isotropic case, $t_x=t_y=t_0$) significantly exceeds $U_c^{\left(\mathrm{CPA}\right)}=3.49t_0$, and therefore the CPA results cannot be considered as upper bound to $U_c$ in a rigorous manner. 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 $|{\psi }_0\rangle$ are optimized. In brief, the SGA ground-state energy is lower or equal to that obtained from GA, leading to $U_c^{\left(\mathrm{SGA}\right)}⩾U_c^{\left(\mathrm{GA}\right)}$. However, the mutual relation between SGA and VMC results cannot be determined a priori, as the former provides better optimization of $|{\psi }_0\rangle$, whereas the latter puts more emphasis on the accurate calculation of averages in Equation (13). Looking at the results presented in Figure 5, we may conclude that $U_c^{\left(\mathrm{GWF}\right)}\gtrsim U_c^{\left(\mathrm{SGA}\right)}$, 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

$$\frac{t_y}{t_x}=(0.0269\pm 0.0014)\times {\left(\frac{U_c}{t_x}\right)}^{2.93\pm 0.05}.$$

(23)

(Here, a single standard deviation is given for each parameter.) The line given by Equation (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 in double occupancies are gradually amplified when the interaction is increased. These effects are further discussed in the 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\to 0$ for $t_y/t_x\to 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 produce $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 Figure 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\to \infty$ 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) compared to the large-system limit. Therefore, one could roughly estimate $U_c^{\left(\mathrm{GWF}\right)}/t_x$ to be reduced by 0.2–0.3 when enlarging the system for $t_y\approx t_x$. This quantity is comparable but smaller than the deviation from the ‘exact’ QMC result of Ref. [18], namely $U_c^{\left(\mathrm{QMC}\right)}-U_c^{\left(\mathrm{GWF}\right)}\approx 0.4t_x$, 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.

Moreover, in Figure 5 (bottom panel), we depict the trajectories followed by a sheet of graphene subjected to a strain in the armchair direction (${\epsilon }_y>0$) and allowed to relax in the perpendicular (i.e., zigzag) direction. The hopping matrix elements in the Hamiltonian (1) are parametrized according to [48,70,71]

$$t_{ij}=-t_0\left(1-\beta \frac{\delta d_{ij}}{d_0}\right),$$

(24)

where $\beta =2$ (red solid line; open symbols) or $\beta =3$ (blue solid line; closed symbols) is the dimensionless electron–phonon coupling parameter. The bond-length variations ($\delta 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

$$U_{\mathrm{eff}}=U-V_{01}{〈\frac{d_0}{d_{ij}}〉}_{j\left(i\right)},$$

(25)

with the coefficients $U=3.63t_0$, and $V_{01}=2.03t_0$ taken from Ref. [23], and ${\langle \cdots \rangle }_{j\left(i\right)}$ denoting the average over three nearest neighbors j of the site i. (Due to our suppositions on the symmetry, Equation (25) produces same value for all sites.)

Depending on the electron–phonon coupling $\beta$, we find the strain of ${\epsilon }_y=0.25$ (for $\beta =2$) or ${\epsilon }_y=0.20$ (for $\beta =3$) is necessary to approach $U_c^{\left(\mathrm{HF}\right)}$, being a conventional border of the CSM range. These values are comparable with the maximal strain of ≈0.20 reported in experiments. The phase of the Mott Insulator seems inaccessible by applying mechanical strains to graphene, although modification of the equilibrium $U/t_0$ ratio due to the substrate effect may possibly enhance the interaction effects. Below, we discuss the effects of electron correlations which should also be visible in CSM (or even SM) phase.

3.2. Effects of Strain on Measurable Quantities

Earlier in this paper, we point out that a model assuming a linear density of states, Equation (11), parametrized by the Fermi velocity a zero energy, gives the critical Hubbard interaction $U_c^{\left(X\right)}$ that differs by only $5\%$ from the values obtained using 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 $\mathsf{\Lambda }$, related to the Fermi velocity. For $t_y⩽t_x$ and strains limited to experimentally-accessible values, one can set $\mathsf{\Lambda }\propto \sqrt{t_x{t}_y}$, leading to

$$U_c\approx U_c^{\left(0\right)}\left(1-\frac{1}{2}{\delta }_t\right),{\delta }_t=\frac{t_x-t_y}{t_x},$$

(26)

where $U_c^{\left(0\right)}$ a zero-strain value. Substituting the values of $U_c^{\left(\mathrm{HF}\right)}$ and $U_c^{\left(\mathrm{GWF}\right)}$ for $t_y=t_x$, we find (in Figure 6) that the evolution of $U_c/t_x$ with increasing strain is approximated by Equation (26) quite well for both HF and GWF methods (similar agreement is observed for GA results, omitted in Figure 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 also be visible in the strained system for $U<U_c$. This expectation is further supported by the data presented in Figure 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 from VMC simulations do not differ significantly from GA results (see datapoints and black solid lines, respectively), while HF (dashed lines) predict qualitatively different behavior, particularly for $\langle \widehat{T}\rangle$ displayed versus U in Figure 7a–c, but the differences between HF and Gutzwiller-based techniques are also apparent for $\langle n_{i\uparrow }{n}_{i\downarrow }\rangle$ (see remaining panels in Figure 7).

For ${\delta }_t=0.5$, see Figure 7c,f, corresponding to the strain of ${\epsilon }_y=0.22$ for $\beta =3$, and in the interval of $U/t_x=1.5÷2$ being relevant for graphene, the values of $\langle \widehat{T}\rangle$ obtained from GWF or GA are reduced by more than $20\%$ comparing the ${\delta }_t=0$ situation; see Figure 7a. (Notice that the above-mention reduction includes the change of $t_x$, used as an energy unit in Figure 7; the details are given in Appendix B.) Moreover, in the intermediate case, ${\delta }_t=0.25$, a $10\%$ reduction is noticed, see Figure 7b. For $\langle n_{i\uparrow }{n}_{i\downarrow }\rangle$, the effect of strain is less pronounced, but we still have an approximately $10\%$ reduction for the ${\delta }_t=0.5$ case [see Figure 7f] compared to the ${\delta }_t=0$ case [Figure 7d], following from both GWF or GA methods for the interval of $U/t_x=1.5÷2$.

4. Concluding Remarks

We have investigated the mutual effect of electron–electron interaction, modeled by a Hubbard term in the second-quantized 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 focused 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, the energy spectrum for the noninteracting system remains gapless for arbitrary high strains, since the Dirac cones do not merge. In turn, the semimetal-insulator transition may occur only due to interactions. Moreover, the system gradually evolves, with increasing strain, towards a collection of weakly-coupled Hubbard chains, allowing us to expect a considerable reduction in 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).

The phase diagram for the parametrized model is supplemented with calculations of trajectories, followed by monolayer graphene strained in the 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 the 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 in double occupancies, which are well-pronounced (and affected by an applied strain) much before the transition. The above-mentioned effects will probably 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-) cost counterparts for the Quantum Monte Carlo methods.