Exact Solutions of Topological Superconductor Model with Hubbard Interactions

The strong correlated systems have been a difficult subject in condensed matter physics for decades Fradkin (2013). Due to the lack of small parameters, the traditional methods such as perturbation expansion or mean field approximation cannot provide reliable calculations for systems with strong enough couplings. In these situations, exact solvable models are very valuable for understanding the behaviors in strong coupling limit. Unfortunately, most exact solvable models are limited to spatial dimension one and they are usually solved by the Bethe ansatz which is quite complicated mathematically Bethe (1931); N. Andrei and Lowenstein (1983); Karabach et al. (1997); Wang et al. (2015). Exact solvable models in higher dimensions are certainly worth close studies.

The past decade witnessed the rising interests in the study of topological matters Hasan and Kane (2010); Qi and Zhang (2011). This new trend also inspired new ideas on constructing exact solvable models in dimension two (2D) Azaria et al. (1987); Sengupta et al. (2008). A famous example is the Kitaev model Kitaev (2006); Baskaran et al. (2007); Mandal and Jayannavar (2020); Miao et al. (2017a) on the honeycomb lattice which supports nontrivial topological order. Although its original form is a spin model, if one rewrite the Kitaev model in terms of Majorana fermions, then it is evident that there are infinite many conserved quantities, which makes this model exact solvable. More recently, this line of developments has been revived and a large class of exact solvable BCS-Hubbard models have been proposed in a series of works Chen et al. (2018a); Li et al. (2019); Miao et al. (2019, 2017b). These types of models are studied under the name of the Falicov-Kimball model many years ago. The common feature of these modes is the appearance of flat bands for certain appropriately chosen parameters. There exist many quadratic terms of flat band fermions which can serve as the infinite conserved quantities. This will transfer the four-fermion interacting Hubbard term to a quadratic fermion term. Then the interacting model can be easily solved as a simple quadratic fermion model. Other than the exact solvability, these models also provide a platform to study the interplay between the strong correlations and topology Ezawa (2018, 2017).

In the present paper, we propose a topological superconductor model with equal spin pairing, based on the Qi-Wu-Zhang model Qi et al. (2006); Asbóth et al. (2016) which is a Chern insulator. In section II, we present the a detailed construction of the model Hamiltonian. We also point out that, for certain choice of pairing gap, one half of its energy bands become flat bands. In section III, we transfer the model Hamiltonian to the real space and introduce the repulsive Hubbard interaction term. In section IV, it is shown that the exact solutions of this model can be easily revealed by rewriting its Hamiltonian in terms of Majorana fermions Shen (2012); Chen et al. (2018b). We also provide detailed studies of the band structure and edge modes in both non-interacting limit and also the fullly interacting case. The Chern number and system energy of the interacting model are computed numerically.

We will briefly review the Qi-Wu-Zhang model, which is prototype of two-dimensional Chern insulators on square lattice. In the momentum space, the Hamiltonian in second quantized form is given by

Here $a,b$ labels the two orbital in each site. $\sigma_{j}$ for $j=1,2,3$ are Pauli matrices. The three coefficients in front of the Pauli matrices is given by

Since the hopping and potential terms of the above model do not depend on the spin, we can treat this model as a spinless fermion model. The Qi-Wu-Zhang model has no internal symmetries and belongs to the class A of the tenfold way classification of topological insulators Schnyder et al. (2008); Chiu et al. (2016). In two-dimension, the topological property of class A is characterized by the $\mathbb{Z}$ classification, which is simply the Chern number of the occupied energy bands. It is well known that for a two-band model, the Chern number can be written as a winding number of the mapping from $k$-space to a unit sphere Qi et al. (2008) as

Here $R=\sqrt{R_{1}^{2}+R_{2}^{2}+R_{3}^{2}}$ and $\hat{R_{i}}=R_{i}/R$ is unit vector along $R_{i}$. Applying the above formula to the Qi-Wu-Zhang model, it is easy to see that the Chern number of lower band depends on the parameter $m$ as follows

The non-trivial Chern number appears when all $R_{i}$ can change signs in the $k$-space.

In order to construct an exact solvable model from the above Chern insulator, we first extend the above model into a topological superconductor by adding certain pairing terms. The pairing terms can be treated as equal spin pairing, which have similar structure as the hopping terms. These pairing terms are specially constructed such that the model can support flat bands for some appropriate choice of parameters. In experiments, this type of pairing terms can be physical realized in the superconductor-topological insulator heterostructure by superconducting proximity effect Fu and Kane (2008); Xu et al. (2014); Sau et al. (2010); Stanescu et al. (2010). To express the superconductor model, we introduce the Nambu spinor as follows

Then the Hamiltonian can be written as the standard Bogoliubov-de Gennes (BdG) form as

Here the kinetic term and pairing terms are given by

where $\sigma_{0}$ is 2 by 2 identity matrix. The topological superconductor has the following particle-hole symmetry

here $\tau_{1}$ is the Pauli matrix applying to the Numbu spinor space. Due to this particle-hole symmetry, this model belongs to the class D of tenfold way classification. In dimension two, the class D is also classified by integer group $\mathbb{Z}$, which is again labeled by the Chern number of the occupied energy bands.

The topological superconductor Hamiltonian in the first quantized form can be written more explicitly as a $4\times 4$ matrix

Here $R_{\pm}=R_{1}\pm iR_{2}$. In this form, it is easy to compute the four energy bands, which are given by

If one assumes that $\Delta=1$, there are two bands become exactly flat with $E=0$. The appearance of flat bands is crucial for construction of an exact solvable model with Hubbard interaction. Since the flat band fermions has zero energy, their fermion operators do not evolve with time and can be treated as constants. Therefore any four fermion interaction term involving two flat band fermions and two other fermions will reduce to a quadratic fermion term. This in turn makes the interacting model exact solvable. We will get back to this point in details in the next section.

Before introducing the Hubbard interaction term, we can first discuss the topological property of the above model. From now on, we will always assume that the flat band condition $\Delta=1$ is satisfied. This situation is also the focus of the next section. In this case, the eigenvector of the lowest band can be analytically expressed as

Plug the above result into the definition of Chern number

After some algebraic calculations, we find the same result of the Chern number as for a two-band model

Therefore, the Chern number of the model Eq.(11) is the same as the Chern number of the Qi-Wu-Zhang model in Eq.(7).

The result of Chern number can also be reflected through the number of edge modes by the bulk-edge correspondence. To investigate the edge modes of the topological superconductor model of Eq.(11), we realize this model on a cylinder-shape lattice. To be concrete, we assume an open boundary condition along $x$ axes and keep the periodic boundary condition along $y$ axes. The resulting Hamiltonian is still a function of $k_{y}$. In the left panel of Figure 1, we plot the eigenvalues as a function of $k_{y}$. We have assumed that $\Delta=1$ which is the flat band condition, and $m=0.8$ corresponding to $C=1$. One can see that beside the two valence and conduction bands, there are two perfect flat bands located at zero energy. Other than the 4 bulk bands, there are two chiral edge modes connecting the valence and conduction bands. In the right panel of Figure 1, the amplitude of one of the edge modes is shown as a function of the lattice sites. One can see the amplitude is localized at one end of the cylinder, which confirms that this is indeed an edge mode.

As a side note, we mention that it is also convenient to rewrite the topological superconductor model in a more compact form as

Here we introduce $\Gamma_{ij}=\tau_{i}\otimes\sigma_{j}$ with $\tau_{i}$ acting on the Nambu spinor space and $\sigma_{j}$ acting on the orbital space.

In order to introduce the four-fermion Hubbard interaction term, we extend the topological superconductor model Eq.(11) to include the spin degree of freedom. In the first quantized form, the Hamiltonian can be written as

Here $s_{0}$ is a $2\times 2$ identity matrix acting on the spinor space and $H_{TS}$ is from Eq.(18).

Since the Hubbard interaction is an on-site repulsion in real space, we can also transfer the hopping and pairing terms of the above model to the real space as follows

Here $n=(n_{x},n_{y})$ labels the lattice site on a square lattice and $s=1,2$ labels the spin up and down. The two-component fermion operator is defined as $c_{n,s}=(c_{n,1,s},c_{n,2,s})^{T}$. The second index labels the two orbital. We also introduce $\hat{x}$ and $\hat{y}$ representing the unit vector along the $x$ and $y$ direction respectively. The spin dependence of above Hamiltonian is trivial, because the spin up and down behave exactly the same. Nevertheless, we still keep the spin index explicit in the above equations, because the Hubbard interaction term we are about to introduce will depend on the spin.

To verify the correctness of Eq.(25) and Eq.(26), we can make the following fourier transformation back to momentum space,

Here $L_{x}$ and $L_{y}$ are the lattice site number along the $x$ and $y$ axes. For the hopping terms, it is straightforward to see that the transformation gives back the Qi-Wu-Zhang model in momentum space. For the pairing terms, we make use of the following identities

Then it is easy to see that Eq.(26) reproduce the desired pairing terms in momentum space.

Now we are in the position to define the Hubbard interaction term. It is an on-site repulsive interaction between spin up and spin down fermions. Explicitly, it can be written as

For $U>0$, the zero and double occupation of each orbital will be punished by raising the system energy with $U/4$. On the other hand, the single occupation is favored by lowering the system energy with $U/4$.

Collect all the above results, we have arrived at the full interacting topological superconductor Hamiltonian in real space

In the next section, we will show that this model can be exactly solved when the flat band condition $\Delta=1$ is satisfied.

The exact solvability of the interacting model Eq.(31) can be best understood by introducing the Majorana fermion operators as follows

Here $s$ labels the spin up and down. It is easy to see that they satisfy the following anti-commutation relations

Here $m,n$ labels the lattice sites, $j,j^{\prime}$ labels the two orbital, $s,s^{\prime}$ labels the spin degree of freedom.

Substituting the Majorana representation into Eq.(31), after some algebraic calculation, we find that the hopping and pairing terms can be written as

Here $\Delta_{1}=1+\Delta$ and $\Delta_{2}=1-\Delta$. In the meanwhile, he Hubbard interaction can be expressed as

When the flat band condition $\Delta=1$ is satisfied, all kinetic terms for $a_{n,j,s}$ in Eq.(35) are vanished. If we define $D_{n,j}=4ia_{n,j,\uparrow}a_{n,j,\downarrow}$, then it is easy to see the $[D_{n,j},H]=0$ for all sites $n$ and $j=1,2$, which makes $D_{n,j}$ a c-number. From the anti-commutation relations, we find $D_{n,j}^{2}=1$, therefore $D_{n,j}=\pm 1$. Suppose the total number of lattice sites is $N$, then there are $2^{2N}$ choices of different configurations of $D_{n,j}$. For a given configuration of $D_{n,j}$, the Hubbard interaction term reduces to a quadratic term of Majorana fermions $b_{n,j,s}$. Therefore, the full interacting model becomes quadratic in fermion operators, which can be easily solved.

Among all the configurations of $D_{n,j}$, there are two regular configurations obviously require more detailed examination. In contrast to the random configuration which can only be solved in real space, the regular configuration allows one to transfer the Hamiltonian back to momentum space, which greatly simplifies the computation of eigenvalues. One of the regular configurations is the uniform distribution with $D_{n,j}=1$ for all sites and orbital. If we transfer to the momentum space and rewrite the Majorana fermion back to ordinary fermion operators, we find that the Hubbard interaction term can be written as

Here $\Gamma_{ij}=\tau_{i}\otimes\sigma_{j}$ with $\tau_{i}$ acting on the Nambu spinor space and $\sigma_{j}$ acting on the orbital space. The matrices $s_{j}$ acts on the spin space. In this case, the 4 energy bands can be simple expressed as

The other 4 are flat bands with $E=0$.

The other regular configuration is the staggered distribution with $D_{n,1}=1$ and $D_{n,2}=-1$ for all sites. Transferring back to the original fermion operators, the Hubbard interaction can be written as

In this case, 4 out of the 8 energy bands can be written as

and the other 4 are zero energy flat bands. We will see later that the staggered configuration is the true ground state of the full interacting model.

Now we can investigate the band structure and edge modes of the interacting model. To this end, we can again realize the above model on a cylinder geometry with open boundary along $x$ axes. Then one can numerically diagonalize the resulting Hamiltonian which depends on the remaining momentum $k_{y}$. For the uniform configuration, we plot the eigenvalues as a function of $k_{y}$ in the left panel of Figure 2. One can see that the interaction term has lifted the degeneracy between the spin up and down fermions. Thus, there are two valence bands below $E=0$ and two conduction bands above zero energy. Each band can be mapped to itself by reflecting about $k_{y}=\pi$. There are clearly two pairs of edge modes signaling the non-trivial topology of these bulk bands. The wave-function of a typical edge mode is shown in the right panel of Figure 2. Its amplitude is localized at the right end, confirming it is truly an edge mode. Similarly, we show the eigenvalues as a function of $k_{y}$ for the staggered configuration in the left panel of Figure 3. Again, there are 4 different energy bands for the two orbital and two spins. In contract to the uniform case, for the staggered case, the two valence bands can be mapped to each other by reflecting about $k_{y}=\pi$. There also appears two pairs of chiral edge modes with its wave-function shown in the right panel of Figure 3.

The number of edge modes is also reflected by the Chern number of the bulk bands. For the interacting model, the eigenvectors are too complicated to support an analytical expression, thus it is impossible to find a simple result of Chern number as in Eq.(3). But it is still straightforward to compute the Chern number numerically for each separate bands. For the fixed interaction strength $U=2$, the numerical result shows that for the uniform case, the Chern numbers of the lower two bands are

Here the $C_{i}$ with $i=1,2$ labels the lowest and second lowest band. $\Theta(x)$ is the step function. For the staggered case, the Chern numbers are

The above results of Chern number indeed match the number of edge modes discussed above.

At last, we want to study the energetics of the model in order to determine which configuration is the ground state of the system. The model is realized on a $10\times 10$ lattice with open boundaries along both directions. Then we randomly choose $D_{n,j}$ and numerically find out all eigenvalues. The system energy is obtained by the summation of all negative eigenvalues. For such a large system, it quite expensive to enumerate all possible disordered configurations. Therefore, in Figure 4, we plot the system energies for 100 randomly chosen disordered configurations. The black and red thick lines corresponds the energies of the staggered and uniform configuration. One clear see that staggered configuration is the true ground state of the system. On the other hand, the uniform configuration gives rise to the highest system energy.

We have constructed a topological superconductor with equal spin pairing and Hubbard repulsive interaction, which can be exactly solved when the pairing gap equals to the hopping constant. The topology of this interacting model is examined in detail. The edge modes are displayed explicitly by putting the model on a cylinder geometry and numerically solve the band structure. It is found that the number of edge modes exactly matches the Chern number of occupied bands, which is consistent with the bulk-edge correspondence. Finally, the true ground state with staggered configuration is determined by comparison with many random configurations.