Disorder effects on the Topological Superconductor with Hubbard Interactions

In the past decade, there are tremendous progress that have been achieved in the area of the topological properties of materials Hasan and Kane (2010); Qi and Zhang (2011). However, most of these studies are restricted to the free fermion systems or weakly interacting systems. It is well known that the understanding of strongly correlated systems such as Hubbard model Arovas et al. (2022); Hirsch (1985); Schulz (1990); Fradkin (2013) is very difficult since the perturbation theory is not reliable for these systems due to the lack of small expansion parameters. Recently, there emerges a series of topological models with Hubbard interactions which are exact solvable due to the existence of flat bands. This type of models are partly inspired by the Kitaev model Kitaev (2006) on the honeycomb lattice and are extensively studied by several authors Chen et al. (2018); Li et al. (2019); Miao et al. (2017, 2019); Ezawa (2018). They provide some interesting examples for understanding how the interactions affect the topological properties. Along this line, we also proposed an exactly solvable topological superconductor with Hubbard interaction Deng and He (2021), which enriches this class of models.

In the real materials, disorders or impurities are un-avoidable. In view of this, the studies of the disordered systems are also an important topic in condensed matter physics. A classical example of disorder effects is the Anderson localization Abrahams et al. (1979) which induces metal-insulator transitions in free fermion systems. An intriguing question is that how the disorders will affect the exact solvable interacting topological models we mentioned above. It is well-known that the topological phase can persist until the disorder strength reaches some critical values. After that, the system will return to the trivial phase. However, it is also found that the disorders can induce nontrivial topology in the system that is trivial in the clean limit. Ever since the pioneer work of the so-called topological Anderson insulator Li et al. (2009); Shen (2012) has been carried out, many other theoretical models have been proposed to understand these disorder induced topological phases Chen et al. (2017a, b); Klier et al. (2019); Guo et al. (2010); Guo (2010); Shapourian and Hughes (2016); Chen et al. (2015). The disorder induced topological phase transitions were also discussed in works Louvet et al. (2016); Ryu and Nomura (2012). In this paper, we will investigate how the disorders influence the topology of our interacting topological superconductor. The exact solvability allows us to perform a detailed study of the interplay between the interactions and disorders. Furthermore, it is also possible to make reliable calculation when the disorders are imposed on the coupling constants.

The rest of this paper is organized as follows. In section II, we present the TS model with Hubbard interactions and show that this model is exactly solvable when flat band condition is satisfied. In section III, we introduce the real space Chern number formula which is suitable for investigations of the topology of disordered systems. The self-consistent Born approximation is also discussed, which can provide more insights about the influence of disorders. The numerical results of the disordered interacting TS model based on the above mentioned two methods are presented in section IV. Finally, we briefly conclude in section V.

In this section, we introduce the model of topological superconductor which can support flat bands with a suitable choice of parameters. Owing to the presence of flat bands, this type of model is exactly solvable when a Hubbard interaction term is included Miao et al. (2017); Ezawa (2018). The Hamiltonian of the topological superconductor (TS) in momentum space can be expressed in the conventional Bogoliubov-de Gennes form as Deng and He (2021)

Here we have defined $\psi_{\mathbf{k}}=(c_{1,\mathbf{k}},c_{2,\mathbf{k}},c_{1,-\mathbf{k}}^{\dagger},c_{2,-\mathbf{k}}^{\dagger})^{T}$. The first sub-indices denote the two different orbitals on each site. The lattice momentum $\mathbf{k}$ takes values from the Brillouin zone of a square lattice. The hopping and pairing terms are given by

Where $\sigma_{j}$ with $j=1,2,3$ are Pauli matrices and $\sigma_{0}$ is 2 by 2 identity matrix. The 3 coefficients $R_{j}$ are chosen as follows

One can see that the hopping term alone is equivalent to the famous Qi-Wu-Zhang model of Chern insulator Qi et al. (2006).

It is also convenient to rewrite the topological superconductor model as

Here we defined $\Gamma_{ij}=\tau_{i}\otimes\sigma_{j}$ with Pauli matrices $\tau_{i}$ and $\sigma_{j}$ acting on the Nambu spinor space and the orbital space respectively. One can see that the 4 energy bands of $H_{TS}$ is given by

Therefore, when $\Delta=\pm 1$, there will appear two flat bands $E=0$. Note that the pairing term is chosen to share the same structure as the hopping term. This is the reason why the $H_{TS}$ can support two flat bands with special choice of parameters.

The topological properties of the model with Eq.(3) is the same as the Qi-Wu-Zhang model. One find that the Chern number $C=\pm 1$ when $|m|<2$ and $C=0$ otherwise. When $C=\pm 1$, the above model is topologically nontrivial, which also indicates that chiral edge modes will appear if we impose an open boundary condition to the system. One the other hand, for $C=0$, the system is topologically trivial with no edge modes.

Up to now, the above topological superconductor is a model of spinless fermions. It is straightforward to include the spin degree of freedom which gives rise to the Hamiltonian of the first quantization form as

Here $s_{0}$ is a $2\times 2$ identity matrix acting on the spinor space and $H_{TS}$ is from Eq. (3). Now we can define the Hubbard interaction as the on-site repulsion between spin up and spin down fermions as follows

Since the Hubbard interaction term is expressed in the real space, it is more convenient to transfer the hopping and pairing terms of Eq.(10) to the real space. In the second quantization form, The resulting $H_{\textrm{CI}}$ and $H_{\textrm{pair}}$ terms can be written as

Here $\mathbf{n}=(n_{x},n_{y})$ labels the lattice site on a square lattice and $s=\uparrow,\downarrow$ labels the spin up and down. The two-component fermion operator is defined as $c_{\mathbf{n},s}=(c_{\mathbf{n},1,s},c_{\mathbf{n},2,s})^{T}$, where the second index labels the two orbitals. We also introduce $\hat{x}$ and $\hat{y}$ representing the unit vector along the $x$ and $y$ direction respectively. Now the full interacting topological superconductor can be summarized as

In order to understand the reason why the above interacting model is exactly solvable, we introduce the Majorana fermion operators as follows

If the Majorana fermion operators is substituted into the hopping and pairing terms of Eq.(12) and (13), we find that the terms involving $a_{\mathbf{n},j,s}$ are all vanished when the flat band condition $\Delta=1$ is satisfied (Detailed calculations can be found in Deng and He (2021)). In terms of the Majorana fermions, the Hubbard interaction can be expressed as

If we define $D_{\mathbf{n},j}=4ia_{\mathbf{n},j,\uparrow}a_{\mathbf{n},j,\downarrow}$, it is easy to see that $[D_{\mathbf{n},j},H]=0$ in all sites and orbitals. Therefore, all $D_{\mathbf{n},j}$ are actually constant $c$-numbers. Due to the anti-commutation relations of Majorana fermions, one can see that $D_{\mathbf{n},j}^{2}=1$ and thus $D_{\mathbf{n},j}=\pm 1$. When the distribution of $D_{\mathbf{n},j}$ is given, the Hubbard interaction term becomes quadratic in terms of fermion operators, which makes the interacting model solvable.

Among all possible the $D_{\mathbf{n},j}$ configurations, there are two special ones with the potential to give rise the ground state of the interacting topological superconductor. The first one is the uniform distribution with $D_{\mathbf{n},j}=1$ for all sites and orbitals, which is called the ferromagnetic (FM) configuration. The other one is the staggered distribution with $D_{\mathbf{n},1}=1$ and $D_{\mathbf{n},2}=-1$ for all sites, which is called anti-ferromagnetic (AFM) configuration.

If we transfer the Majorana fermions back to ordinary fermions, the Hubbard interaction term in the first quantized form for the FM and AFM configurations can be written as

With this result, the interacting topological superconductor of Eq.(14) can be converted into a simple free fermion model, for which the Chern number can be easily computed. Accordingly, we map out the phase diagram of the interacting topological superconductor in Fig.1. The topological non-trivial phase is roughly located inside the area of $m\in(-2,2)$ and $U\in(-4,4)$. The yellow and blue regions indicate the phase with Chern number $C=\pm 2$, while the green region indicate the topological trivial phase with $C=0$. We would like to mention that the true ground state of the model is the AFM configuration.

In this paper, we will focus on the disorder effects of the interacting topological superconductor of Eq.(14). We always assume that the flat band condition $\Delta=1$ is satisfied, such that the interacting model can be exactly solved. In order to keep the solvability of this interacting model, we impose the disorders either on the coefficient $m$ or the coupling constant $U$. More explicitly, we will make the following replacement $m\to m_{\mathbf{n}}=m+\omega_{\mathbf{n}}$ or $U\to U_{\mathbf{n}}=U+\omega_{\mathbf{n}}$. In other words, we consider two types of disorder models as follows

Here $\omega_{\mathbf{n}}$ located at each lattice site are random variables with uniform distribution inside the interval $[-W/2,W/2]$. $W$ is the disorder strength. We have numerically checked that the AFM configuration is also the true ground state even when the disorders are introduced. Therefore, the disorder potential term in Eq.(21) takes the form of Eq.(19).

In this section, we discuss two different methods to investigate the topology of the interacting TS model with disorders. Usually, the Chern number is defined as an integral in the momentum space, but this form is not suitable for disordered system. Here we introduce a Chern number formula in real space, which can compute this topological index to high precision. The other method is the so-called self-consistent Born approximation, which can convert the disorder system into a uniformed system with renormalized parameters. Then, it is convenient to understand how the disorder affects topological phase transitions.

The Chern number formula defined in $\mathbf{k}$-space can be written as

Here $P_{\mathbf{k}}$ is the projection operator in the $\mathbf{k}$-space, which is define as

where $|\psi_{n}\rangle$ represents the $n$-th eigenvector with the eigen-energy $E_{n}$, and $E_{f}$ represents the Fermi energy. This operator $P_{\mathbf{k}}$ projects the wave-function into the subspace spanned by the filled energy bands.

For the disordered systems, it is not possible to directly make $\mathbf{k}$-space calculations. Therefore, one must seek a method to compute the Chern number in real space. If we diagonalize the Hamiltonian in real space, we can also construct a projection operator $P$ in real space. Next, we can approximate the derivative by finite differences. Then we will arrived at a Chern number formula in real space. This approach was suggested a few years ago by Prodan Prodan et al. (2010); Prodan (2011). We will provide a detailed derivation of this formula in the appendix A. In this section, we simply give out the real space Chern number formula as follows

Here the dot means element-wise product of the two matrices. “Im” means taking the imaginary part of the following quantities. The matrices $K_{x}$ and $K_{y}$ are given in the appendix A.

Now we turn to the self-consistent Born approximation which is widely used in the many works such as Chen et al. (2017a, b); Klier et al. (2019). This method usually provide a useful way to understand topological phase transitions driven by disorders. The essence of this method is to incorporate the effects of the disorder into a self-energy. Then the self-energy will in turn modify the parameters of the original model. Now the disordered model is actually converted into a uniform model with renormalized parameters. Then we can determine the boundary of topological phase transition by the phase diagram of Fig. 1 with renormalized $m$ or $U$. In other words, the way that renormalized $m$ or $U$ depends on the disorder strength $W$ directly reflect the influence of disorders on topological phase transitions.

In the self-consistent Born approximation, the self-energy $\Sigma$ is defined as Shen (2012)

Here $H_{\textrm{full}}(\mathbf{k})$ is the interacting TS model of Eq.(14) in the momentum space. $V_{\mathbf{n}}$ is the disorder potential in real space. The $\langle\cdots\rangle_{\textrm{dis}}$ denotes the disorder average. For the two disorder models of Eq.(20) and (21), the disorder potential can be written as

where $\omega_{\mathbf{n}}$ is the random number we defined before.

If we expand Eq.(25) to the first order, we find that the disorder-induced self-energy is determined by the equation Groth et al. (2009); Hua et al. (2019) as follows

Here $E_{F}$ is the Fermi energy which is close to zero in our model. The integration is over the whole Brillouin zone.

For the disorder model of Eq.(20), by iterating Eq.(27), we can solve the self-energy and find that $\Sigma=\delta m\,\Lambda_{1}$ plus some constant matrices which can be ignored. Now the effective Hamiltonian is $H=H_{\textrm{full}}+\Sigma$. Therefore, one can see that the self-energy only renormalize the parameter $m$ to $\widetilde{m}$, which is given by

Similarly, if we consider the disorder model of Eq.(21), we find that the self-energy $\Sigma=\delta U\,\Lambda_{2}$ and the coupling $U$ is renormalized to $\widetilde{U}=U+\delta U$.

In this section, we present the numerical results generated by the two methods that we discussed in the previous section. We start from with the real space Chern number for the disordered interacting topological superconductor. First, as a simple example, we consider the disorder model of Eq.(20) with $U=0$. In other words, we are working with Eq.(10) without the Hubbard interaction and the disorders are imposed on the potential term as in Eq.(20). Here we study the disorder effects on two different cases with $m=1.2$ or $m=2.2$ corresponding to topological or trivial phases of the clean TS model respectively. The Chern number as a function of disorder strength $W$ is shown in Fig. 2.

For the case of $m=1.2$ (blue curve), we have the Chern number $C=2$ in the clean limit. As the disorder strength $W$ increases, one can see that $C$ stays at value $2$ for a while and then start to drop at about $W=4.8$. It gradually decreases to zero at about $W=9$ and after that $C$ stays at zero. Therefore, there exists a critical disorder region where the topological phase transition takes place. This indicates that when the disorders are strong enough, they can kill the topological nontrivial phases.

For the other case of $m=2.2$ (red curve), the model is topologically trivial when $W=0$. However, when the disorder strength increased, one can see that the Chern number rises from $0$ to $2$ at around $W=1.7$ and stays at the value 2 for a while. It then starts to drop at $W=5.8$ and back to zero at around $W=8$. Within a certain range of disorder strength, a plateau of non-zero Chern number is maintained, which indicates that disorder can drive the system into a topologically nontrivial phase. This is very similar to the so-called topological Anderson insulator. But here it is realized in a topological superconductor instead of insulator.

Next we turn on the Hubbard interaction in the disorder model of Eq.(20). Since we assume $\Delta=1$, this mode is still exactly solvable. The result of its real space Chern number is presented in the left panel of Fig.3. Here we still consider two different cases with $m=1.2$ (blue curve) and $m=2.2$ (red curve). For both cases, we assume $U=1.2$. The behavior is very similar to the non-interacting case of Fig.2. For $m=1.2$, the model is already topological for $W=0$. When $W$ increase to around $W=5$, the Chern number $C$ gradually drop to zero. One the other hand, the model is trivial for $m=2.2$ in the clean limit. But the disorders will drive $C$ to increase from zero to $2$ and then return to zero again. The plateau of $C=2$ within a certain range of $W$ suggests that the disorder can induce a topologically nontrivial phase even when the Hubbard interaction is presented. This result generalizes the phenomena of topological Anderson transition to interacting models.

Since the interacting TS model is exact solvable, we can also impose the disorders on the coupling constant, as the disorder model of Eq.(21). The Chern number of this model is shown in the right panel of Fig. 3. The blue and red curve corresponds to the $U=1.2$ and $U=4.1$ respectively. In both cases, we assume $m=1.2$. One can see that the $U=1.2$ system has $C=2$ in the clean limit and $C$ gradually decreasing to zero inside the region of $20<W<30$. According to the phase diagram of Fig.1, the system with $U=4.1$ is topologically trivial in the clean limit. Again, one can see that the disorders can generate a stable plateau of $C=2$ inside the region of $10<W<20$. This is again the signature of topology driven by disorders. In the model of Eq.(21), the disorders come from the random coupling constant, which is usually difficult to analyze in other types of models.

We also would like to mention that the error bars in both Fig. 2 and 3 represent the standard deviations of 30 random samples. The error bars in the transition regions are larger than the ones in the stable regions.

At last, we present the numerical results of effective medium theory based on the self-consistent Born approximation. We first consider the disorder model of Eq.(20) which renormalizes the parameter $m$. The renormalized $\widetilde{m}$ as a function of $W$ is displayed in the top left panel of Fig.4. One can see that for both $m=2.2,\,U=0$ (blue triangle) and $m=2.2,\,U=1.2$ (yellow dot) cases, the two renormalized $\widetilde{m}$ are almost identical. The dashed line represents the topological phase boundary of $\widetilde{m}=2$. The corresponding critical disorder strength is $W\approx 2$, which is close the region where $C=0$ rises to $C=2$ in the left panel of Fig.3. For larger $W$, we find that $-2<\widetilde{m}<2$, which suggest the system stays in the topological phase even when $W$ increases up to 10. On the other hand the Chern number curve in the left panel of Fig.3 shows that $C$ drop to zero at $W\approx 6$. This contradiction shows that the self-consistent Born approximation only works for the relatively weak disorders. It cannot capture the correct topological features when $W$ is too strong.

Similarly, we can also study the disorder model of Eq.(21) which renormalizes the coupling $U$. The resulting $\widetilde{U}$ is shown in the right panel of Fig.4. The dashed line represents the phase boundary of $\widetilde{U}=3.9$. One can see the corresponding critical $W\approx 10$ which is close to the first transition in the right panel of Fig.3. Again, the renormalized $\widetilde{U}$ fail to capture the second transition at strong $W$.

Another advantage of the self-consistent Born approximation is that it can eliminate the disorders such that the model becomes a uniform one again. Because of this, we can put the model on a cylinder geometry with open boundary along $x$-axes and closed boundary along $y$-axes. For example, the disorder model of Eq.(21) is converted to Eq.(14) with renormalized $\widetilde{U}$. Now the Hamiltonian as a function of $k_{y}$ is ready to be diagonalized. The resulting band structure is shown in the bottom panel of Fig.4. One can see the appears of the edge modes due to the non-trivial topology. It is also evident that the system preserves its flat bands even with the introduction of disorders, thereby maintaining its exact solvability. Similar conclusions can be drawn for the disorder model of Eq.(20).

We have studied the topological phase transitions of the disordered topological superconductor with Hubbard interaction where the disorders are either imposed on the potential term or the coupling constants. Throughout the whole paper, we stick to the flat-band condition which guarantees the exact solvability of this model even when disorders are included. It is found that the strong enough disorders are able to suppress the topological phases. For topological trivial system in the clean limit, the disorders also induce a phase with non-zero Chern number. This can be thought as a generalization of the topological Anderson transitions to interacting models. The results of real space Chern number is also confirmed by the self-consistent Born approximation. In the later method, we also demonstrated the influence of disorders on the system parameters and band structures.