Generalization of Lieb-Wu wave function inspired by one-dimensional ionic Hubbard model

TL;DR

This paper generalizes the Lieb-Wu wave function inspired by the ionic Hubbard model, providing a new analytical approach that is valid in certain regimes and useful for cold atom experiments.

Contribution

It introduces a generalized Bethe ansatz wave function incorporating an ionic parameter, extending the Lieb-Wu solution to the ionic Hubbard model.

Findings

The generalized wave function reduces to Lieb-Wu in the limit of zero ionic parameter.

The two-particle scattering matrix satisfies the Yang-Baxter equation.

Numerical solutions of the generalized Bethe equations yield ground state energies in the thermodynamic limit.

Abstract

With the ionic Hubbard model (IHM) in mind, we construct a non-trivial generalization of the Bethe ansatz (BA) wave function which naturally generalizes the Lieb-Wu wave function with an ionic parameter $\Delta$, and reduces to Lieb-Wu solution in the limit $\Delta\to 0$. The resulting two-particle scattering matrix satisfies the Yang-Baxter equation. To the extent that the unit cells with more than two electrons (Choy-Haldane issue) are avoided on average, our wave function represents an effective soluiton for the one-dimensional IHM. The Choy-Haldane issue restricts the validity of our solution to low-filling and large $U\gtrsim 4$. This regime is attainable in cold atom realizations of the IHM. For this regime, we numerically solve the generalized Bethe equations and compute the ground state energy in the thermodynamic limit.

Equations

$$H_{\mathrm{IHM}}=H_{\mathrm{t}}+H_{\mathrm{U}}+H_{\Delta },$$

$$H_{\mathrm{t}}=$$

$$H_{\mathrm{U}}=$$

$$H_{\Delta }=$$

$$S^x=\frac{1}{2}\sum _{j=1}^L\sum _{\eta }\left(c_{j,\uparrow }^{{\eta }^{†}}{c}_{j,\downarrow }^{\eta }+c_{j,\downarrow }^{{\eta }^{†}}{c}_{j,\uparrow }^{\eta }\right),$$

$$S^y=\frac{i}{2}\sum _{j=1}^L\sum _{\eta }\left(c_{j,\downarrow }^{{\eta }^{†}}{c}_{j,\uparrow }^{\eta }-c_{j,\uparrow }^{{\eta }^{†}}{c}_{j,\downarrow }^{\eta }\right),$$

$$S^z=\frac{1}{2}\sum _{j=1}^L\sum _{\eta }\left(c_{j,\uparrow }^{{\eta }^{†}}{c}_{j,\uparrow }^{\eta }-c_{j,\downarrow }^{{\eta }^{†}}{c}_{j,\downarrow }^{\eta }\right).$$

$$|F⟩=$$

$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{{\boldsymbol{\Sigma}}}({\boldsymbol{X}},{\boldsymbol{\eta}})\equiv F_{\sigma_{1}\ldots\sigma_{N}}(x_{1},\eta_{1};\ldots;x_{N},\eta_{N})=\mathcal{A}\prod_{j=1}^{N}e^{ik_{j}x_{j}}\chi_{\eta_{j}}(\theta_{j})\sum_{Q}A_{\sigma_{1},...,\sigma_{N}}(Q)\tilde{\mathbb{H}}(x_{Q})}$$

$$h_N$$

$$\left[-2t\cos (\frac{k}{2})\left(\cos (\frac{k}{2}){\tau }^x+\sin (\frac{k}{2}){\tau }^y\right)+\Delta {\tau }^z\right]F_{\sigma }^{\eta }(x)=E{F}_{\sigma }^{\eta }(x),$$

$$\begin{array}{rl}& F_{\sigma }^{\eta }(x)=A_{\sigma }{e}^{ikx}{\chi }_{\eta }(\theta ),\\ & E(\theta )=\Delta \cosh (\theta ),\\ & k(\theta )=2{\cos }^{-1}(\frac{\Delta }{2t}\sinh (\theta )).\end{array}$$

$${\chi }_b(\theta )=-\varkappa (\theta )\frac{\cosh \left(\frac{\theta }{2}\right)}{\sqrt{\cosh (\theta )}},{\chi }_a(\theta )=\frac{\sinh \left(\frac{\theta }{2}\right)}{\sqrt{\cosh (\theta )}}.$$

$$F_{{\sigma }_1,{\sigma }_2}(x_1,{\eta }_1;x_2,{\eta }_2)=$$

$$\left(A_{{\sigma }_1,{\sigma }_2}\tilde{\mathbb{H}}(x_1-x_2)+B_{{\sigma }_1,{\sigma }_2}\tilde{\mathbb{H}}(x_2-x_1)\right)$$

$$\times e^{i(k_1{x}_1+k_2{x}_2)}{\chi }_{{\eta }_1}({\theta }_1){\chi }_{{\eta }_2}({\theta }_2)$$

$$-\left(A_{{\sigma }_2,{\sigma }_1}\tilde{\mathbb{H}}(x_2-x_1)+B_{{\sigma }_2,{\sigma }_1}\tilde{\mathbb{H}}(x_1-x_2)\right)$$

$$\times e^{i(k_1{x}_2+k_2{x}_1)}{\chi }_{{\eta }_1}({\theta }_2){\chi }_{{\eta }_2}({\theta }_1),$$

$$\text{S}=\frac{1}{2}\left(\text{I}+\text{P}\right)+\frac{1}{2}\left(\text{I}-\text{P}\right)\text{s},$$

$${\text{s}}^{\text{b}\text{b}}=\frac{{\Xi }^{\text{b}\text{b}}-\frac{U}{t}}{{\Xi }^{\text{b}\text{b}}+\frac{U}{t}},{\text{s}}^{\text{a}\text{a}}=\frac{{\Xi }^{\text{a}\text{a}}-\frac{U}{t}}{{\Xi }^{\text{a}\text{a}}+\frac{U}{t}},$$

$$\begin{array}{r}{\text{s}}^{\text{b}\text{b}}={\text{s}}^{\text{a}\text{a}}=\frac{\left(\sin \left(\frac{k_1}{2}\right)-\sin \left(\frac{k_2}{2}\right)\right)-i\frac{U}{2t}}{\left(\sin \left(\frac{k_1}{2}\right)-\sin \left(\frac{k_2}{2}\right)\right)+i\frac{U}{2t}},\end{array}$$

$${\text{s}}^{\text{b}\text{a}}={\text{s}}^{\text{a}\text{b}}=\frac{0-\frac{U}{t}}{0+\frac{U}{t}}.$$

$$\begin{array}{rl}\text{s}=& (|{\chi }_b({\theta }_1){|}^2|{\chi }_b({\theta }_2){|}^2{\Xi }^{\text{b}\text{b}}+|{\chi }_b({\theta }_1){|}^2|{\chi }_a({\theta }_2){|}^2\times 0\\ & +|{\chi }_a({\theta }_1){|}^2|{\chi }_b({\theta }_2){|}^2\times 0+|{\chi }_a({\theta }_1){|}^2|{\chi }_a({\theta }_2){|}^2{\Xi }^{\text{a}\text{a}}-U/t)\\ & /(|{\chi }_b({\theta }_1){|}^2|{\chi }_b({\theta }_2){|}^2{\Xi }^{\text{b}\text{b}}+|{\chi }_b({\theta }_1){|}^2|{\chi }_a({\theta }_2){|}^2\times 0\\ & +|{\chi }_a({\theta }_1){|}^2|{\chi }_b({\theta }_2){|}^2\times 0+|{\chi }_a({\theta }_1){|}^2|{\chi }_a({\theta }_2){|}^2{\Xi }^{\text{a}\text{a}}+U/t).\end{array}$$

$$\begin{array}{rl}{\text{S}}_{i,j}=& \frac{\left(\alpha ({\theta }_i)-\alpha ({\theta }_j)\right){\text{I}}_{i,j}+i\frac{2U}{\Delta }{\text{P}}_{i,j}}{\left(\alpha ({\theta }_i)-\alpha ({\theta }_j)\right)+i\frac{2U}{\Delta }}.\end{array}$$

$$\begin{array}{rl}{\text{S}}_{i,j}=& \frac{\left(\alpha (k_i)-\alpha (k_j)\right){\text{I}}_{i,j}+i\frac{U}{t}{\text{P}}_{i,j}}{\left(\alpha (k_i)-\alpha (k_j)\right)+i\frac{U}{t}}.\end{array}$$

$$F_{{\sigma }_1,...,{\sigma }_N}(x_1,{\eta }_1;...;x_N,{\eta }_N)$$

$$=\mathcal{A}\prod _{j=1}^N{e}^{i{k}_j{x}_j}{\chi }_{{\eta }_j}({\theta }_j)\sum _{Q}S(Q)A_{{\sigma }_1,...,{\sigma }_N}\tilde{\mathbb{H}}(x_Q).$$

$$\begin{array}{rl}{F}_{{\sigma }_1,...,{\sigma }_N}& (...;x_j=1,{\eta }_j;...)\\ & =F_{{\sigma }_1,...,{\sigma }_N}(...;x_j=L+1,{\eta }_j;...).\end{array}$$

$$\begin{array}{r}(Z_j{)}_{{\sigma }_1,...,{\sigma }_N}^{{\sigma }_1^{\prime },...,{\sigma }_N^{\prime }}{A}_{{\sigma }_1^{\prime },...,{\sigma }_N^{\prime }}(Q)=e^{-i{k}_{j}L}{A}_{{\sigma }_1,...,{\sigma }_N}(Q),\end{array}$$

$$\begin{array}{r}{e}^{i{k}_{j}L}=\prod _{\delta =1}^M\frac{{\Lambda }_{\delta }-{\alpha }_j-ic/2}{{\Lambda }_{\delta }-{\alpha }_j+ic/2},j=1,2,...,N\end{array}$$

$$\begin{array}{r}-\prod _{\delta =1,\delta \ne \gamma }^M\frac{{\Lambda }_{\delta }-{\Lambda }_{\gamma }+ic}{{\Lambda }_{\delta }-{\Lambda }_{\gamma }-ic}=\prod _{j=1}^N\frac{{\Lambda }_{\gamma }-{\alpha }_j-ic/2}{{\Lambda }_{\gamma }-{\alpha }_j+ic/2},\gamma =1,2,...,M\end{array}$$

$$\begin{array}{r}E=-\sum _{j=1}^N\sqrt{{\Delta }^2+4t^2{\cos }^2\left(\frac{k_j}{2}\right)},P=\sum _{j=1}^N{k}_j.\end{array}$$

$$L{k}_j=2\pi I_j+\sum _{\delta =1}^M\Theta (2{\alpha }_j-2{\Lambda }_{\delta }),j=1,2,...,N$$

$$\sum _{j=1}^N\Theta (2{\alpha }_j-2{\Lambda }_{\gamma })=2\pi I_{\gamma }+\sum _{\delta =1}^M\Theta ({\Lambda }_{\delta }-{\Lambda }_{\gamma }),\gamma =1,2,...,M$$

$$\begin{array}{r}\rho (k_j)=\frac{1}{L(k_{j+1}-k_j)},\sigma ({\Lambda }_{\gamma })=\frac{1}{L({\Lambda }_{\gamma +1}-{\Lambda }_{\gamma })}.\end{array}$$

$$\begin{array}{r}{\int }_{-Q}^Q\rho (k)dk=\frac{N}{L},{\int }_{-B}^B\sigma (\Lambda )d\Lambda =\frac{M}{L}.\end{array}$$

Full text

1 INTRODUCTION

Bethe Ansatz (BA) is a powerful method to construct exact wave functions for the ground and excited states of vast classes of one-dimensional Hamiltonians [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. It provides solutions of Heisenberg model [15], Lieb-Liniger Model [16], Hubbard model [17, 18], XXZ [19, 20], XYZ [21, 22], Temperley-Lieb spin chain [23], Thirring model [24, 25] and Chiral-Invariant Gross-Neveu Hamiltonian  [9].

Among all these, the celebrated Hubbard model addresses the competition between kinetic energy and the on-site Coulomb repulsion between electrons of opposite spins [26, 27]. An interesting deformation of this model is the so called ionic Hubbard model [28]. The IHM is the Hubbard model plus an additional alternating (ionic) scalar potential of strength $\Delta$ (see Fig. 1). The essential phyiscs of this model is the competition between the inoic potential $\Delta$ and the Hubbard term $U$. This model in arbitrary dimension and at half-filling has been the subject of extensive studies with various techniques [29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]: In infinite dimension the method of dynamical mean field theory (DMFT) has been applied to investigate the phase transitions of the IHM at half filling [40, 54, 53, 42, 45, 36, 51]. For large $U$ ($\Delta$) the Mott (band insulating) phase is stabilized. While with cellular DMFT one obtains a direct transition between Mott and band insulating phases [51], within a single-site DMFT one obtains an intermediate conducting phase, i.e. when $U$ and $\Delta$ become comparable, the parent conducting states is restored [40, 45, 53]. This agrees with coherent potential treatment of Hoang [49]. A similar picture can be obtained from continuous unitary transformations [43]. This method has been extensively applied by Hafez and co-workers to study the excitaiton sepctrum of this model [47, 48, 47, 50, 55, 57].

In two space dimensions the picture remains unclear. Determinantal quantum Monte Carlo study gives a metallic phase [41]. Building on an orthogonal metallic state [61], interesting possibility of semiconductor of spinons has been proposed [60]. On triangular lattice, crossovers between Mott and charge transfer and covalent insulators and magnetically ordered states is reported [44]. Aligia used charge and bond operator formalism to study the phase diagram of IHM in various dimensions [34]. Mancini has employed his composit operator method to study IHM [38].

In one dimension, Fabrizio and co-workers used the bosonization to address the competition between $U$ and $\Delta$ at half-filling. They find an spontaneously dimerized phase [35] between the Mott and band insulating phase [29]. This agrees with numerical calculations of Otsuka and Nakamura based on level-crossing and renormalization group method based on exact diagonalization [62]. IHM in 1D was subsequently studied with density matrix renormalization group (DMRG) [31, 33, 46] and Quantum Monte Carlo [30, 52] which support the bond ordered state. This model has also been studied using an interesting two-site entropy as a diagnosis tool of phase transitions [39]. The opitcal response of an extension of IHM has also been studied by Maeshima and Yonemitsu [37].

Given the recent cold atom realization of IHM [56] and the invention of superlattice modulation spectroscopy to probe the bond order [59], and a very appealing recent proposal to realized 1D IHM in twisted bilayer GeSe [63], an exact solution of this model in 1D would be desirable. The first question one faces in this direction is the quantum integrability of the model. In our previous work [hosseinzadeh2019integrability], using energy level spacing statistics as well as the eigenvalues of the one-particle density matrix, we have established the many-body localization for IHM. Therefore, there should exist an exact solution. However the Choy-Haldane issue [64] prevents the Bethe ansatz from solving the IHM: Beyond two electrons per unit cell (containing two sites or ”orbitals” with energies $\pm\Delta$) in Fig. 1, the BA wave function will not be eigen-state of the IHM [64, 65, 66]. The essence of Choy-Haldane issue is that, there are configurations in the Hilbert space of the actual IHM that are not accessible by Bethe ansatz wave functions. In the limit of large Hubbard $U$ which prevents the double occupancy of every orbital, the Choy-Haldane issue is ”effectively” avoided at low-filling. Furthermore, in the $\Delta=0$ limit, there is no Choy-Haldane issue. Therefore it is hoped that the Bethe ansatz wave function is a reasonable approximation of the exact eigen functions of the IHM in the following two regimes: (i) small $\Delta$ and arbitrary filling and (ii) large $U$ and low filling and the Choy-Haldane issue will causes parametrically small effects in the above regimes.

2 MODEL HAMILTONIAN

Hamiltonian describing the one-dimensional IHM, is an innocent-looking extension of the Hubbard model in which the lattice sites are assigned alternating electrostatic potentials $\pm\Delta$ as shown in Fig. 1. The Hamiltonian of IHM is given by

[TABLE]

where $H_{\rm t}+H_{\rm U}$ and $H_{\Delta}$ are the standard Hubbard Hamiltonian and the ionic part respectively which are described by,

[TABLE]

The parameters $t,U$ and $\Delta$ are hopping integral between two adjacent sublattices, on-site Hubbard interaction and ionic terms, respectively. In this paper, we set the energy unit by $t=1$. The superscripts $a,b$ are introduced to label the two sublattices with ionic potentials $\pm\Delta$, respectively. Therefore $c_{j,\sigma}^{\eta^{\dagger}}$ and $c_{j,\sigma}^{\eta}$ are the creation and annihilation operators for an electron of spin $\sigma$ in the Wannier state at the $j$th unit cell as in Fig. 1 and $\eta=a,b$ labels the internal orbital degree of freedom of every unit cell that corresponds to two sublattices. The Fermionic anticommutation relations is given by $\{c_{j,\sigma}^{\eta^{\dagger}},c_{j^{\prime},\sigma^{\prime}}^{\eta^{\prime}}\}=\delta_{j,j^{\prime}}\delta_{\eta,\eta^{\prime}}\delta_{\sigma,\sigma^{\prime}}$. They act in a Fock space with the pseudovacuum $|{0}\rangle$ defined by $c_{j,\sigma}^{\eta}|{0}\rangle=0$.

The model has translational symmetry in the $j$ and rotational invariance in spin space, with the corresponding SU(2) generators given by,

[TABLE]

The number operators $\hat{n}_{\uparrow}=\sum_{j=1}^{L}\sum_{\eta}c_{j,\uparrow}^{\eta^{\dagger}}c_{j,\uparrow}^{\eta}$ and $\hat{n}_{\downarrow}=\sum_{j=1}^{L}\sum_{\eta}c_{j,\downarrow}^{\eta^{\dagger}}c_{j,\downarrow}^{\eta}$ correspond to $\uparrow$ and $\downarrow$ spins, from which the total number of electrons will become, $\hat{n}=\sum_{j=1}^{L}\sum_{\eta,\sigma}c_{j,\sigma}^{\eta^{\dagger}}c_{j,\sigma}^{\eta}$. Obviously the commutation relations $\left[H,\hat{n}\right]=\left[H,\hat{n}_{\uparrow}\right]=\left[H,\hat{n}_{\downarrow}\right]=0$ hold that imply the total number of electrons $N$, down-spin electrons, $M$, and up-spin electrons, $(N-M)$, are constants. Therefore we can label the eigenstates in sectors specified by $N$ and $M$.

In at attempt to overcome the Choy-Haldane issue [64], a projection of kinetic energy term into a subspace with no more than two electrons per unit cell was proposed [65, 66]. In our case it will amount to define a constrained IHM (CIHM) by $H_{\rm CIHM}={\cal P}H_{\rm t}{\cal P}+H_{\rm U}+H_{\Delta}$ where ${\cal P}$ is the projection operator into the subspace with no more than two electrons per unit cell. This deformation of IHM restricts the Hilbert space in a way that it avoids the Choy-Haldane issue. But unfortunately it spoils the $U\to 0$ limit. In fact the $U\to 0$ limit of the Bethe ansatz wave function falls outside the portion of Hilbert space of the CIHM. However, within the IHM itself, when the Hubbard $U$ is large enough, and the filling is low-enough, the double occupancy of individual orbitals will be effectively avoided. Therefore our Bethe wave function obtained in this paper will be effective description of the spectrum of IHM in this limit.

In the following we will construct a generalization of Bethe wave function of Lieb and Wu for non-zero ionic potential $\Delta$. We show that its two-particle scattering matrix satisfies the Yang-Baxter equation. This wave functions allows us to explore strong $U$ limit of the IHM at low filling.

3 IHM MOTIVATED BETHE ANSATZ WAVE FUNCTION

In this section we discuss the coordinate Bethe ansatz wave function inspired by IHM. In this sense it will be a generalization of Lieb-Wu wave function for nonzero $\Delta$. From this equation we construct the two-particle scattering matrix and show that it satisfies the Yang-Baxter equation. Let the Hilbert space of $N$ electrons be $\mathcal{H}_{N}$ that contains a vacuum state $|0\rangle$. In this Hilbert space, we search for wave functions of the following Bethe ansatz form,

[TABLE]

where $j=1,\ldots,N$ labels the electrons, $x_{j}$ and $\sigma_{j}$ are position and spin of $j$’th electron, $\eta_{j}=a,b$ is the sublattice label of the $j$’th electron and

[TABLE]

is the BA wave function in coordinate basis the various components of which are further specified in the following. The Schrödinger equation $H|F\rangle=E|F\rangle$ in the coordinate basis becomes $h_{N}F_{\boldsymbol{\Sigma}}({\boldsymbol{X}},{\boldsymbol{\eta}})=EF_{\boldsymbol{\Sigma}}({\boldsymbol{X}},{\boldsymbol{\eta}})$ with the first quantized Hamiltonian $h_{N}$ given by,

[TABLE]

where $\tau_{j}^{i}$ is the $i$th ($i=0,x,y,z$) unit and Pauli matrix acting in sublattice space $a,b$ of the $j$’th electron.

In the one-particle sector ($N=1$) the Schrödinger equation becomes,

[TABLE]

whose solution are,

[TABLE]

where $\theta$ is the rapidity of the particle, $k$ is the wave vector, $A_{\sigma}$ is a spinor in spin space and $\chi_{\eta}(\theta)$ is a spinor in sublattice space whose components are of the following form,

[TABLE]

In this equation, $\varkappa(\theta)=\frac{\Delta}{2t}\sinh(\theta)+i\sqrt{\frac{4t^{2}}{\Delta^{2}}-\sinh^{2}(\theta)}$. There are two bands with positive and negative energies. Eq. (6) is for positive energy states. The negative energy states are simply obtained by $\theta\to i\pi-\theta$.

In construction of BA for models (such as Hubbard model) with only one orbital degree of freedom per unit cell, one breaks up the region of spatial variables $(x_{1},\ldots,x_{N})$ into $N!$ corners obtained by operating permutation $Q$ of indices $1,2\ldots,N$ on $x_{1}<x_{2}\ldots<x_{N}$ and looking for a wavefunction in each of these corners. For $N=2$ there are $N!=2$ corners, namely $x_{1}<x_{2}$ and $x_{2}<x_{1}$. In the case of IHM where there is additional orbital degree of freedom labelled by $\eta=a,b$, we need to slightly adjust the notation to take proper care of the spatial ordering of electron coordinates. In this case the spatial coordinates are $(x_{1},\eta_{1};...;x_{N},\eta_{N})$ and a modified step function will be introduced to handle the ordering of the actual location of electrons. Let us see how this is done for $N=2$ electrons. Higher $N$ values is similarly treated. The Hamiltonian for two electrons is $h_{2}=\sum_{j=1}^{2}\left(-2t\cos\left(\frac{k_{j}}{2}\right)\left(\cos\left(\frac{k_{j}}{2}\right)\tau_{j}^{x}+\sin\left(\frac{k_{j}}{2}\right)\tau_{j}^{y}\right)+\Delta\tau_{j}^{z}\right)+\frac{U}{2}\delta_{x_{1},x_{2}}\left(\tau_{1}^{0}\otimes\tau_{2}^{0}+\tau_{1}^{z}\otimes\tau_{2}^{z}\right)$. The Hubbard interaction $U$ operates when $x_{1}=x_{2}$ (i.e. the particles are in the same unit cell), and further $\eta_{1}=\eta_{2}$ meaning that they are on the same sublattice. Away from the boundary of the corners the Hamiltonian is free and the wavefunction can be written in the form of product of single particle solutions,

[TABLE]

Here $\tilde{\mathbb{H}}(x_{j}-x_{l})=\mathbb{H}(x_{j}-x_{l})+\frac{\delta_{x_{j},x_{l}}}{2}\left(\delta_{\eta_{j},b}\delta_{\eta_{l},a}-\delta_{\eta_{j},a}\delta_{\eta_{l},b}\right)$ is the extended Heaviside step function that separates the two corners by incorporating their orbital degree of freedom $\eta=a,b$ of the $j$’th and $l$’th electrons. This wave function is explicitly antisymmetric. Furthermore, $\mathbb{H}$ is the ordinary Heaviside step function and is regularized as by defining $\mathbb{H}(0)=\frac{1}{2}$. The total energy and total momentum for this wave function are given by $E=E_{1}+E_{2}=\Delta\left(\cosh(\theta_{1})+\cosh(\theta_{1})\right)$ and $k=k_{1}+k_{2}=2\left(\cos^{-1}(\frac{\Delta}{2t}\sinh(\theta_{1}))+\cos^{-1}(\frac{\Delta}{2t}\sinh(\theta_{2}))\right)$, respectively, that are conserved in all corners.

We introduce scattering matrix relating to the spin amplitudes in two corners of Eq. (8) as $B_{\sigma_{1},\sigma_{2}}=S_{\sigma_{1},\sigma_{2}}^{\sigma^{\prime}_{1},\sigma^{\prime}_{2}}A_{\sigma^{\prime}_{1},\sigma^{\prime}_{2}}$. This is determined by imposing two conditions [67]:

($i$) Continuity of the wavefunction in the corner boundary (Uniqueness). Considering the condition of the continuity of the wavefunction at $x_{1}=x_{2}$ and $\eta_{1}=\eta_{2}$ leads to,

[TABLE]

where $\textbf{I}_{\sigma_{1},\sigma_{2}}^{\sigma^{\prime}_{1},\sigma^{\prime}_{2}}=\delta^{\sigma^{\prime}_{1}}_{\sigma_{1}}\delta^{\sigma^{\prime}_{2}}_{\sigma_{2}}$ and $\textbf{P}_{\sigma_{1},\sigma_{2}}^{\sigma^{\prime}_{1},\sigma^{\prime}_{2}}=\delta^{\sigma^{\prime}_{2}}_{\sigma_{1}}\delta^{\sigma^{\prime}_{1}}_{\sigma_{2}}$ are the identity and exchange operators in the spin space, respectively. In the ordinary Hubbard model where there is one orbital degree of freedom per unit cell, s has to be a scalar (which is called spin scalar parameter in the BA literature). But in the present case where we have two orbitals (sublattices) $a,b$ in each unit cell, s will be a matrix in the orbital space that completely characterizes the scattering process.

($ii$) The Schrödinger equation on the corner boundary ($x_{1}=x_{2};\eta_{1}=\eta_{2}$). The Schrödinger equation for two electrons has four components labelled with sublattice indices $\textbf{bb},\textbf{ba},\textbf{ab}$ and aa. For two of the four cases, namely, ($\textbf{bb},\textbf{aa}$) Hubbard interactions occur between electrons. Therefore s for these components are obtained as follows (for details see appendix A),

[TABLE]

where rapidity related parts are $\mathfrak{\Xi}^{\textbf{b}\textbf{b}}=\left(\frac{\chi_{a}(\theta_{2})}{\chi_{b}(\theta_{2})}\left(1-e^{ik_{2}}\right)-\frac{\chi_{a}(\theta_{1})}{\chi_{b}(\theta_{1})}\left(1-e^{ik_{1}}\right)\right)$, and $\mathfrak{\Xi}^{\textbf{a}\textbf{a}}=\left(\frac{\chi_{b}(\theta_{1})}{\chi_{a}(\theta_{1})}\left(1-e^{-ik_{1}}\right)-\frac{\chi_{b}(\theta_{2})}{\chi_{a}(\theta_{2})}\left(1-e^{-ik_{2}}\right)\right)$. It can be seen that, by changing the sign of $\Delta$, these two amplitudes transform to each other, namely, $\textit{{s}}^{\textbf{b}\textbf{b}}(\Delta)=\textit{{s}}^{\textbf{a}\textbf{a}}(-\Delta)$ as expected. In the limit ($\Delta\rightarrow 0$) the scalar spin parameters in Eq. (10) reduce to,

[TABLE]

which has the expected form for the ordinary Hubbard model, except that the argument of the $\sin$ function is $k_{j}/2$ rather than $k_{j}$. This is because in the present formulation each unit cell contains two degrees of freedom corresponding to which the Brillouin zone will be folded into half of the case where only one degree of freedom in each unit cell is considered. For the other two cases ($\textbf{ba},\textbf{ab}$), the Hubbard interaction does not appear in their equations. In this case the parameter s can be written as $\textit{{s}}^{\textbf{b}\textbf{a}}=\textit{{s}}^{\textbf{a}\textbf{b}}=\zeta-\xi$, with the constraint that $\zeta+\xi=1$. This is because when two electrons do not interact with each other, they can either pass each other without change (scattering matrix is unit matrix I), or exchange (scattering matrix is solely P [67]). Any linear combination of the above two cases is also acceptable and this is the source of degree of arbitrariness in choosing $\zeta$ and $\xi$. The only constraint is $\zeta+\xi=1$ which stems from the solution of Schrödinger equation at the boundary. By choosing $\zeta=0$ and $\xi=1$, the spin scalar parameter will be $\textit{{s}}^{\textbf{b}\textbf{a}}=\textit{{s}}^{\textbf{a}\textbf{b}}=-1$ and we can write them as follows:

[TABLE]

Finally, by combining various components of the spin scalar in Eqs. (10) and (12), and by taking into account the probability of each $(a,b)$ components of the electrons, the spin scalar parameter becomes,

[TABLE]

By inserting Eq. (13) into Eq. (9) and after some algebra, and using the definition $\alpha(\theta)\equiv\tanh(\theta)\sqrt{4t^{2}/\Delta^{2}-\sinh^{2}(\theta)}$, the two particle scattering matrix will be,

[TABLE]

By changing the variable from the rapidity $\theta$ to momentum $k$ according to Eq. (6), for particles with positive energies we have $\alpha(k)\equiv\frac{2t\cos(k/2)\sin(k/2)}{\sqrt{\Delta^{2}+4t^{2}\cos^{2}(k/2)}}$ which then gives,

[TABLE]

From Eq. (15) it is clear that the limit of $\Delta\to 0$, this scattering matrix is exactly the same as the scattering matrix obtained for the standard Hubbard model. The only difference is that $U$ has been replaced with $2U$. This is because the present scattering matrix is for a unit cell containing an orbital degree of freedom $\eta=a,b$. Thus, each electron in interaction with another electron is experiencing a total potential of $2U$.

A straightforward calculation shows that the above $S$ matrix obtained for the IHM satisfies the Yang-Baxter equation (for details see appendix B). Therefore the following BA wavefunction for the $N$-particle system is consistent at three-particle scattering level and therefore it describes the exact wave function of some integrable model,

[TABLE]

Here $Q$ is a corner corresponding to one set of permutations, $S(Q)$ is the corresponding product of $S$-matrices, and $A_{\sigma_{1},...,\sigma_{N}}$ is the spin amplitude in some reference corner. The extended Heaviside function is defined below Eq. (8) which incorporates the sublattice structure into the definition of corners. The large $U$ and low-filling limit of the model described by the above wave function, approximately corresponds to the IHM.

Now we would like to study the system of $0<N\leqslant L$ electrons on a finite length $L$ with periodic boundary conditions. We can impose this condition in different regions in configuration space to the wavefunction as,

[TABLE]

For the configuration $Q$($x_{1}<x_{2}<...<x_{N}$) the above relation is translated to conditions on the spin amplitude of wavefunction in a single corner as,

[TABLE]

where $Z_{j}=S^{j,j-1}S^{j,j-2}\ldots S^{j,1}S^{j,N}S^{j,N-1}\ldots S^{j,j+2}S^{j,j+1}$, is the so called the transfer matrix. Again, we are encountered with a eigenvalue problem of diagonalization of the spin Hamiltonian ($Z_{j}$) by using of the quantum inverse scattering method for a state with $M$ down spins and $N-M$ up spins and obtain Bethe Ansatz equations as follows [17, 67],

[TABLE]

where $c=\frac{2U}{\Delta}$, $\alpha_{j}\equiv\alpha(\theta_{j})$ is defined in Eq. (14) and $M$ labels the number of down spins in a specific configuration. The spin variables $\Lambda_{\gamma}$ are coupled with momenta $k_{j}$ and must be determined self-consistently as detailed in the next section.

4 THERMODYNAMIC LIMIT AND GROUND STATE PROPERTIES

As already mentioned, we are going to solve these equations for the $0<N\leqslant L$ where $N$ is the total number of electrons and $L$ is the number of unit cells. We will consider the filling factors ($n=N/L\leqslant 1$) of less than a quarter filled bands. This amounts to avoiding the Choy-Haldane issue on average which makes our solution relevant to the IHM. From here on, we write equations in the wave vector space instead of rapidities. At low enough temperatures, only the lower band will be relevant and therefore, the $\alpha_{j}$ parameter in Eqs. (19) will be in the form of $\alpha_{j}=\frac{2t\cos(k_{j}/2)\sin(k_{j}/2)}{\sqrt{\Delta^{2}+4t^{2}\cos^{2}(k_{j}/2)}}$ where $c=U/t$. We will be interested in the limit $N,M,L\to\infty$ with $n=N/L$ and $M/L$ held constant.

For a finite system with $N$ electrons the total energy and momentum of the system given by,

[TABLE]

Eigenvalues and eigenfunctions of the system are characterized by set of momenta and spin variables $\left\{k_{j},\Lambda_{\gamma}\right\}$ as follows: Taking the logarithm of Eqs. (19) we obtain,

[TABLE]

where $\Theta(x)=-2\tan^{-1}(x/c)$ with condition $-\pi<\Theta(x)<\pi$. The $\left\{I_{j},I_{\gamma}\right\}$ are integer or half-integer quantum numbers which labels the states. Using Eqs. (21) the momentum of the system becomes, the $P=\frac{2\pi}{L}\left(\sum_{j}I_{j}+\sum_{\gamma}I_{\gamma}\right)$. Solutions of Eqs. (19) in the thermodynamic limit, based on string hypothesis [68, 67], can have different forms, one of which is the real $k_{j}$’s and real $\Lambda_{\gamma}$’s. This is the case for the ground state of the system.

According to the Lieb-Mattis theorem [69] the ground state of the system on a bipartite lattice should be a total singlet. Therefore in our one-dimensional case, the ground state configuration will be specified by the consecutive arrangement of $\left\{I_{j}\right\}$ and $\left\{I_{\gamma}\right\}$ around zero. In the thermodynamic limit and for any consecutive arrangements of $\left\{I_{j}\right\}$ and $\left\{I_{\gamma}\right\}$ around zero, we define the densities of solutions $\rho(k)$ and $\sigma(\Lambda)$ by [3],

[TABLE]

The values of $\rho(k_{j})$ and $\sigma(\Lambda\gamma)$ are continuously distributed in intervals $\left[-Q,Q\right]$ and $\left[-B,B\right]$, respectively. The integration limits $Q\leqslant\pi$ and $B\leqslant\infty$ are specified by the following two conditions:

[TABLE]

We use the subscript $g$ in the following to indicate that the quantity under consideration is associated with the ground state of the system.

Algebraic manipulation of Eqs. (21) for solution densities at the ground state, gives two coupled integral equations [67]

[TABLE]

where $f(k)=\frac{d\alpha(k)}{dk}$, $K_{n}(x)=\frac{1}{\pi}\frac{(nc/2)}{(nc/2)^{2}+x^{2}}$ and its Fourier transform is $K_{n}(\omega)=\frac{1}{\sqrt{2\pi}}e^{-n\frac{c}{2}|\omega|}$. For the ground state, the $z$-component and total spin is equal to zero ($S_{z}=S=\frac{1}{2}(N-2M)$) such that one has $B_{g}=\infty$. Then by using the Fourier transform of $\sigma_{g}(\Lambda)=1/\sqrt{2\pi}\int_{-\infty}^{\infty}d\omega\sigma_{g}(\omega)e^{i\Lambda\omega}$ for $\rho_{g}(k)$ we have,

[TABLE]

In this equation, $R(x)\equiv(1/4\pi)\int_{-\infty}^{\infty}dy\frac{e^{i\frac{xy}{2}}}{1+e^{|y|}}=(1/\pi)\sum_{n=1}^{\infty}(-1)^{n+1}\frac{2n}{x^{2}+(2n)^{2}}$ is a real valued function. Thus the ground state energy can be written as,

[TABLE]

For quarter filling ($N/L=1$) one has $Q_{g}=\pi$, and the density of electron states ($\rho_{g}(k)$) and ground state energy can be analytically calculated as,

[TABLE]

where $E(\nu)=\int_{0}^{\pi/2}dx\sqrt{1-\nu\sin^{2}(x)}$ , is the complete elliptic integral of the second kind. The functions $\mathcal{J}_{0}$ and $\mathcal{J}_{1}$ are defined by $\mathcal{J}_{0}(\omega)\equiv\pi^{-1}\int_{0}^{\pi}dk\cos(\omega\alpha(k))$, and $\mathcal{J}_{1}(\omega)\equiv\omega\pi^{-1}\int_{0}^{\pi}dk\cos(\omega\alpha(k))f(k)\sqrt{\Delta^{2}+4t^{2}\cos^{2}\left(\frac{k}{2}\right)}$. In the limit $\Delta\to 0$, these functions reduce to the standard Bessel functions $J_{0}$ and $J_{1}$ of the Lieb-Wu solution [70]. Therefore our expression for the ground state energy in the limit of $\Delta\to 0$, correctly reduces to the Lieb-Wu result.

To obtain the ground state energy for different conditions, we numerically solve the equations (23), (25) and (26). Fig. 2 shows our results for ground-state energy. In panels (a) and (b) we plot the ground state energy as a function of $U$ and $n$ for $\Delta=0$ that corresponds to the standard Hubbard model [71]. In panel (c) we increase the ionic potential to $\Delta=2.0$. In the second column we plot the ground state energy as a function of $\Delta$ and $n$ for a fixed $U$. The values of $U$ in panels (d) and (e) are $U=4.0$ while in panel (f) we have a quite strong Hubbard $U=8.0$. First of all, the $\Delta=0$ results in (a) and (b) agree with the numerical results for the Hubbard model [71]. As can be seen in panel (b), at low densities the chemical potential (slope of the curve) is not much dependent on $U$. But by increasing $n$ towards quarter filling, the $U$ dependence of the slope (chemical potential) becomes important. This is in contrast to panel (c) with $\Delta=2.0$ where even near the quarter filling the slope is not much sensitive to the values of $U$. Near the quarter-filling ($n=1$) there is in average one electrons in a unit cell which occupies the negative $\Delta$ site. Therefore for large enough $\Delta$, one barely faces a situation with two electrons coming across each other at the same site and that is why the $U$ dependence upto quarter filling is not strong.

Comparing panels (e) and (f) corresponding to $U=4.0$ and $U=8.0$, respectively, it can be seen that the main factor controlling the shape of the curves is $\Delta$, and the $U$-dependence is already very weak. It appears that the large $U$ regime of validity of our wave function already extends down to $U\sim 4$.

In Fig. 3 we have plotted the ground state energy as a function of $\Delta$ and $U$ for a fixed density $n=1$ corresponding to quarter filling. As can be inferred from Eq. (1), the quantity $\partial E/\partial\Delta$ corresponds to the average sublattice-imbalance, namely the difference in the average population of the two sublattices. As expected in the $\Delta=0$ limit, irrespective of the value of $U$, due to translational invariance the average population of two sublattices must be the same. Therefore the slope of $E(\Delta)$ curve has to be zero. This can be clearly seen in both panels, and in particular in panel (b). For quite large values of $\Delta$, all curves with Hubbard $U$ values upto $8.0$ converge to the same slope $\approx 1$. This simply means that when $\Delta$ is very large, one sublattice (with negative $\Delta$) is fully occupied, while the other sublattice is completely empty. Another interesting feature that can be seen in both (a) and (b) panels is that, for a constant finite $\Delta$ the grounds state energy reduces by increasing the Hubbard $U$. This is better seen for smaller $\Delta$. This can be explained as follows: In small $U$ limit and at quarter (or even less) filling the electrons dominantly occupy negative $\Delta$ sites and positive $\Delta$ sites are empty. Due to kinetic term, electrons may come to positive $\Delta$ sites. Strong $U$ can take advantage of this and lower the ground state energy by superexchange mechanism. In this way the positive $\Delta$ sites will act like effective ligand sites for the superexchange mechanism. As such, superexchange is suppressed by large $\Delta$.

5 SUMMARY AND OUTLOOK

In this work, in an attempt to find an exact solution, we generalized the Lieb-Wu solution of the standard Hubbard model to include the effect of non-zero ionic potential $\Delta$. Every unit cell in IHM will have an additional sublattice label, $\eta=a,b$. This acts like an additional orbital degree of freedom for every unit cell, albeit with different energies $\pm\Delta$. We constructed the two-particle scattering matrix in the space of $\eta$ and spin. For this purpose we introduced extensions of the step function to properly define various ”corners” in the configuration space. The scattering matrix of the IHM satisfies the Yang-Baxter equation which indicates that our Bethe ansatz wave function, actually describes some integrable model (related to IHM). Even without knowing this model, our wave function is a highly non-trivial wave function the properties of which is worth investigation.

Using our Bethe Ansatz wave function, we constructed the ground state energy as a function of various parameters $n,\Delta,U$. In the $\Delta=0$ limit, our results agree with the standard Hubbard model. This wave function can be used to study various correlations in the ground state, as well as the excitation spectrum with particular attention to the spin and charge excitations and the effect of $\Delta$ in these excitations.

Since in every unit cell of the IHM there are two atomic orbitals with on-site energies $\pm\Delta$ (see Fig. 1), by Choy-Haldane issue there will be configurations in the Hilbert space which are not reachable by Bethe wave functions. In the limit of large $U$, and at low fillings, this issue becomes less important. As such, our generalization of the Lieb-Wu wave function represents an effective solution of the IHM in the large $U$ and low-filling limit. Since the IHM has already been realized in a new solid state system based on twisted bilayer GeSe [63], as well as in cold atom setting [56] where $U$ is arbitrarily tunable, our effective exact solution can be a useful tool to extract information about the properties of these systems.

In the absence of an analytic solution for IHM, our generalization of the Lieb-Wu wave function is a promising playground to learn about the competition between the Hubbard $U$ and the ionic potential $\Delta$ in two directions: (i) An outstanding problem would be to search for solutions outside the family of Bethe wave functions (i.e. some sort of generalizations of our $\Delta$-dependent Bethe ansatz wave function) that are able to overcome the Choy-Haldane issue. (ii) Although it would be nice to know the Hamiltonian to which our wave function corresponds, even without knowing what is the Hamiltonian, the wave-function itself contains a non-trivial deal of information about the competition between two important insulating phases, namely the Mott and band insulators.

ACKNOWLEDGMENTS

This work was supported by grant number G960214 of the research deputy, Sharif University of Technology. SAJ was supported by Iran Science Elites Federation (ISEF). We are grateful to Natan Andrei for many insightful comments. AH thanks R. Ghadimi for help in numerical computations. SAJ thanks Perimeter Institute of theoretical physics for hosting a visit during which this work was completed.

APPENDICES

Appendix A: Scattering matrix

For two electrons, the wavefunction for the sublattice components is,

[TABLE]

The Schrödinger equation for this wavefunction is given by the following matrix relation,

[TABLE]

where

[TABLE]

Then in the case $x_{1}=x_{2}\equiv x$, for component $F_{\sigma_{1},\sigma_{2}}(x,b;x,b)$ we have,

[TABLE]

By doing the calculations and using $\textbf{S}=\frac{1}{2}\left(\textbf{I}+\textbf{P}\right)+\frac{1}{2}\left(\textbf{I}-\textbf{P}\right)\textbf{{s}}^{bb}$, the $\textbf{{s}}^{bb}$ component of spin scalar becomes,

[TABLE]

Repeating these calculations for the $F_{\sigma_{1},\sigma_{2}}(x,a;x,a)$ component, the $\textbf{{s}}^{aa}$ component of spin scalar obtain as,

[TABLE]

Appendix B: Yang-Baxter Equation

The Yang-Baxter equation or the star-triangle relation is a consistency equation which if it is established for a scattering matrix, the system (Hamiltonian) governing the dispersion will be integrable. In the case of our problem, this equation is written in the form,

[TABLE]

where the scattering matrix is,

[TABLE]

Assuming the two-particle scattering matrix at the brief form of $\textbf{S}_{i,j}\equiv\Upsilon_{ij}\textbf{I}+\Omega_{ij}\textbf{P}$, and by performing matrix multiplication the left-hand side(LHS) of Eq. (34) is,

[TABLE]

Performing these calculations in a similar manner for the right-hand side(RHS) of Eq. (34) indicates that RHS is exactly equal to LHS.