Topological Phases of an Interacting Majorana Benalcazar–Bernevig–Hughes Model

Abstract

We study the effects of Coulomb repulsive interactions on a Majorana Benalcazar–Bernevig–Huges (MBBH) model. The MBBH model belongs to the class of second-order topological superconductors ($HOTS{C}_2$), featuring robust Majorana corner modes. We consider an interacting strip of four chains of length L and perform a density matrix renormalization group (DMRG) numerical simulation based on a tensor-network approach. Study of the non-local fermionic correlations and the degenerate entanglement spectrum indicates that the topological phases are robust in the presence of interactions, even in the strongly interacting regime.

1. Introduction

In recent years, higher-order topological superconductors (HOTSCs) have attracted significant interest and have been investigated in depth as novel platforms to realize topological superconductivity [1,2,3,4,5]. Aside from their theoretical interest, such systems are potentially relevant for applications; in particular, devices based on two-dimensional HOTSCs hosting Majorana corner states have been proposed to implement braiding dynamics [6,7,8] and quantum gates for quantum computation [9,10,11,12]. HOTSCs have surface states that propagate along one-dimensional lines (hinges) or are localized at some points (corners) on the surface. In particular, for such systems, m-dimensional Majorana corner states can be realized in d-dimensional superconductors, with $m\le d-2$.

In Ref. [13], we introduced a model of second-order topological superconductor ($HOTS{C}_2$) based on Majorana fermions (Mfs) operators. The model is the equivalent of the Benalcazar–Bernevig–Hughes (BBH) model [14] for Dirac fermions. In the Majorana BBH model (MBBH), $C4$ symmetry and reflection symmetries ensure robustness of the corner states. In fact, the model belongs to the trivial two-dimensional $BDI$ class satisfying time-reversal, particle–hole and chiral symmetries [15]. Moreover, the crystalline symmetries ensure a quantized two-dimensional Zak phase [16,17,18]. When written in terms of Dirac fermions, the MBBH model is also equivalent to a model of Kitaev chains coupled by a staggered pairing coupling.

On the other hand, one of the main challenges in the study of quantum matter concerns the robustness of topological phases and topological superconductivity in the presence of interactions. The question was addressed by considering the effect of Coulomb repulsive interactions [19,20,21] and/or by developing a number-conserving theory [22] in order to go beyond the BCS approximation. In particular, some recent studies investigated the physical properties of HOTSCs when Coulomb interactions are taken into account [23,24].

The mean-field approximation offers a “cartoon” description of the topological phases within the single-particle approximation. On the other hand, adding interactions removes the single-particle approximation and opens the way to the investigation of more realistic settings. Thus motivated, in the present work, by means of a density matrix renormalization group (DMRG) numerical analysis, we consider the effects of repulsive Coulomb interactions on a specific version of the MBBH model, namely an interacting strip with four chains.

The strip setting is particularly interesting because it realizes the minimal model able to characterize the crossover from a 1D to a 2D system. It can thus be reduced to an effective one-dimensional lattice that is a suitable lattice for simulations with DMRG technique implemented by matrix product states (MPS) and, on the other hand, it can shed light on the topological robustness expected when one considers fully two-dimensional geometries.

The results of our careful analysis suggest that the quasi one-dimensional system supports two robust zero energy fermionic modes localized at the two edges of the strip and robust, even in the limit of strong interactions. Although the topology of the system is intrinsically two-dimensional, the quasi one-dimensional case provides signatures of the robustness of the topological regime against repulsive interactions. In particular, these robust fermionic modes are expected to split into corner Majorana modes for a fully two-dimensional square lattice system with length L equal to the width N. The twofold degenerate entanglement spectrum, reported in Section 3.2, is a further possible signature of topological order. A complete proof of such a robustness requires to test the system in the fully 2D geometry using three tensor networks, and we expect to accomplish this task in the near future.

Despite being characterized by a specific set of parameters, our system could be realized within the framework of network models [25,26,27,28], instead of traditional condensed matter setups. The experimental implementation of the former was successfully realized using meta-materials platforms, as optical fibers and coupled ring resonators [29,30,31]. A theoretical proposal of network models has already been introduced to describe some condensed matter systems made of chiral Majorana modes [25] or to describe physical phenomena, such as the quantum Hall effect [32]. Network models enable a large degree of control on the parameters subject to various constraints, including particle–hole symmetry. In fact, network models with particle–hole symmetric spectra were already probed in recent experiments [33,34].

The paper is organized as follows. In Section 2, we review the main properties of the mean-field MBBH model introduced in Ref. [13]. In Section 3, we study the particular case of an interacting four-chain strip; in Section 3.1, we discuss the analytical operations needed to transform the strip into an effective one-dimensional system and to perform the DMRG simulations with MPSs. In Section 3.2, we discuss the main numerical results. We track the fermionic correlations along the strip geometry, and we study the trend of the entanglement spectrum from the perturbative to the strongly interacting regime. Conclusions are drawn in Section 4. Appendix A contains details on the matrix product operator (MPO) tensors built for our model.

2. Majorana BBH Model

In this section, we briefly discuss the main topological properties of the BBH model. This provides the appropriate starting point before treating the case of the interacting four-chain strip. Following Ref. [13], we consider a system of Majorana fermions with staggered couplings (w, v) confined in a two-dimensional lattice and described by the Hamiltonian:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_0=& \frac{i}{2}[w\sum _{m,l=1}^{L,N}{a}_{m,l}{b}_{m,l}+v\sum _{m,l=1}^{L-1,N}{b}_{m,l}{a}_{m+1,l}+w\sum _{l=1}^{N-1,2}\sum _{m=1}^L\left(b_{m,l}{b}_{m,l+1}-a_{m,l}{a}_{m,l+1}\right)\hfill \\ & +v\sum _{l=2}^{N-1,2}\sum _{m=1}^L\left(b_{m,l}{b}_{m,l+1}-a_{m,l}{a}_{m,l+1}\right)].\hfill \end{array}\end{array}$$

(1)

$a_{m,l}$ and $b_{m,l}$ are the Majorana operators belonging to a complex fermion operator $c_{m,l}=(a_{m,l}+i{b}_{m,l})/2$. L and N are, respectively, the length and the width of the system, and m and l are the lattice sites. The synthetic $\pi$ flux per plaquette gives rise to the staggered couplings via the Peierls substitution [18].

A scheme of the lattice is reported in Figure 1a, where the translational invariant unit cell is indicated by the blue square plaquette surrounding the Majorana operators a and b of Equation (1), which, on the other hand, are represented by two circles of different colors. When Mfs operators are expressed in terms of fermion operators, the model reduces to Kitaev chains coupled by a staggered pairing coupling:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_0^{\prime }=& \sum _{m,l=1}^{L,N}\mu c_{m,l}^{†}{c}_{m,l}+\sum _{m,l=1}^{L-1,N}(t{c}_{m,l}^{†}{c}_{m+1,l}+\Delta c_{m,l}{c}_{m+1,l}+h.c.)\hfill \\ & +\sum _{m,l=1}^{L,N-1}({\Delta }_1{c}_{m,l}{c}_{m,l+1}+h.c),\hfill \end{array}\end{array}$$

(2)

with $\mu =w$, $t=\Delta =-v/2$ and ${\Delta }_1$:

$$\begin{array}{c}\hfill {\Delta }_1=\left\{\begin{array}{c}-iw,l=odd\hfill \\ -iv,l=even.\hfill \end{array}\right.\end{array}$$

(3)

The system belongs to the second-order topological superconductors (HOTSC2) class, featuring corner Majorana states [13]. The simultaneous presence of crystalline symmetries and standard symmetries ($𝒞$, $𝒫$, $𝒯$) is essential to ensure the realization of second-order topological superconductivity. Indeed, the system satisfies chiral ($𝒞$), particle–hole ($𝒫$) and time-reversal ($𝒯$) symmetries, belonging to the trivial two-dimensional $BDI$ class of the ten-fold classification [15]. Nonetheless, robust corner Majorana modes still appear because of the $C_4$ symmetry, i.e., a rotation of $\theta =\pi /2$ exchanging $x\to y$, $y\to -x$ and reflection symmetries $m_x$, $m_y$, corresponding, respectively, to $x\to -x$, $y\to -y$.

The phase diagram can be obtained by projection of the 2D Zak phase [17,18], which is quantized as $\mathbf{P}=(1/2,1/2)$ in the topological phase and vanishes in the trivial phase ($\mathbf{P}=(0,0)$). On the other hand, the Berry connection vanishes. The phase diagram is reported in Figure 1b, where the green and red colors correspond respectively to the topological and the trivial regimes. The band structure along the crystallographic path $M\Gamma XM$, in the first Brillouin zone, features topological phase transition points. Indeed, the spectrum closes when $\left|w\right|=\left|v\right|$, while it remains gapped both in the trivial and nontrivial phases when $\left|w\right|\ne \left|v\right|$ (see Figure 1c,d). In the topological regime, the zero-energy Majorana modes are localized at the four corners of the strip geometry and decay exponentially inside the bulk, as shown in Figure 1e.

3. Strip of Four Interacting Chains

In order to go beyond the mean field treatment, we add repulsive Coulomb interactions to the MBBH model with the aim of gaining insight on how the interactions affect the stability of the topological phases, which is still an open problem despite intense investigations [23,35,36,37].

Here, we restrict our analysis to a quasi one-dimensional limit consisting of $N=4$ chains of length L. From a computational point of view, the DMRG method [38,39] works at its best when applied to one-dimensional many-body systems associated to gapped Hamiltonians with short-range interactions. For this reason, the simplified limit of four interacting chains analyzed by DMRG techniques based on a tensor-network approach [40,41,42] reduces the computational efforts and increases the power of the numerical algorithms. The strip limit is expected to shed light on the topological properties of the fully two-dimensional Majorana BBH model, which realizes the thermodynamic limit of our strip ($N=2$, L).

3.1. Analytical Results

We include a local repulsive interaction of strength $U>0$ to the Hamiltonian in Equation (2), $H=H_0^{\prime }+H_I$, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_0^{\prime }=& \sum _{m,l=1}^{L,4}\mu c_{m,l}^{†}{c}_{m,l}+\sum _{m,l=1}^{L-1,4}(t{c}_{m,l}^{†}{c}_{m+1,l}+\Delta c_{m,l}{c}_{m+1,l}+h.c.)+\sum _{m,l=1}^{L,3}({\Delta }_1{c}_{m,l}{c}_{m,l+1}+h.c)\hfill \\ \hfill H_I=& \sum _{m,l=1}^{L-1,4}U\left(n_{m,l}{n}_{m+1,l}\right),\hfill \end{array}\end{array}$$

(4)

with $n_{m,l}=c_{m,l}^{†}{c}_{m,l}$. The efficiency of the DMRG algorithm with MPSs for one-dimensional systems [41] stems from the internal structure of the MPSs state [42]. Therefore, we transform the strip geometry of four chains with L sites into a single chain made of $4L$ sites, as shown pictorially in Figure 2a. Introducing a single fermionic species $c_j$ ($c_j^{†}$) and defining a “curvilinear abscissa” which maps the inter-chain interaction terms into nearest-neighbor interaction terms and the intra-chain interaction terms into fourth-neighbour interaction terms, see Figure 2b, the Hamiltonian in Equation (4) reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H=& \sum _{j=1}^{4L}\mu c_j^{†}{c}_j+\sum _{j=1}^{4L-4}(t{c}_j^{†}{c}_{j+4}+\Delta c_j{c}_{j+4}+h.c.)+\hfill \\ & \sum _{j=1}^{4L-1}({\Delta }_1{c}_j{c}_{j+1}+h.c)+\sum _{j=1}^{4L-4}U\left(n_j{n}_{j+4}\right).\hfill \end{array}\end{array}$$

(5)

As fermionic operators of different sites anti-commute, lattice fermionic fields are nonlocal, and it is therefore impossible for such fields to determine a local matrix representation, i.e., a representation in which matrices representing the fermions commute when belonging to different lattice sites. The correct on-site mapping is provided by the Jordan–Wigner transformation, a highly non-local mapping between fermionic operators and spin $1/2$ operators that is a particular case of the general Klein transformation in quantum field theory [43]. On each site, an empty fermionic occupation is mapped into an up spin and an occupied one into a down spin. The nonlocal part of this mapping is the so-called Jordan–Wigner string and fixes the (anti)commutation relation between sites by counting the parity of overturned sites to the left of the spin on which it is applied. This transformation explicitly breaks the translational invariance of the model by singling out a particular site as the origin for each lattice string:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{c}_j=e^{-i\pi {\sum }_{l=1}^{j-1}{c}_l^{†}{c}_l}{\sigma }_j^{+}\hfill \\ c_j^{†}={\sigma }_j^{-}{e}^{i\pi {\sum }_{l=1}^{j-1}{c}_l^{†}{c}_l}\hfill \\ n_j=\frac{1-{\sigma }_j^z}{2},\hfill \end{array}\right.\end{array}$$

(6)

where the aforementioned string parity of the overturned sites is $e^{-i\pi {\sum }_{l=1}^{j-1}{c}_l^{†}{c}_l}$. The operators ${\sigma }_j^{(+,-)}=({\sigma }_j^x\pm i{\sigma }_j^y)/2$ are the well-known combination of Pauli matrices, and the last relation in Equation (6) allows to express the parity operator of the fermionic site j as $e^{-i\pi c_j^{†}{c}_j}={\sigma }_j^z$. Using the algebra of the spin-$1/2$ operators and imposing commutation of the Pauli matrices defined on different lattice sites, it is straightforward to show that

$$\begin{array}{c}\hfill \begin{array}{cc}& c_j{c}_{j+1}=-{\sigma }_j^{+}{\sigma }_{j+1}^{+}\hfill \\ & c_j^{†}{c}_{j+4}={\sigma }_j^{-}\left(\prod _{l=j+1}^{j+3}{\sigma }_l^z\right){\sigma }_{j+4}^{+}\hfill \\ & c_j{c}_{j+4}=-{\sigma }_j^{+}\left(\prod _{l=j+1}^{j+3}{\sigma }_l^z\right){\sigma }_{j+4}^{+}.\hfill \end{array}\end{array}$$

(7)

The Hamiltonian Equation (5) in terms of spin-$1/2$ operators becomes

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H=& \sum _{j=1}^{4L}\mu n_j+\sum _{j=1}^{4L-4}(t{\sigma }_j^{-}{\sigma }_{j+1}^z{\sigma }_{j+2}^z{\sigma }_{j+3}^z{\sigma }_{j+4}^{+}-\Delta {\sigma }_j^{+}{\sigma }_{j+1}^z{\sigma }_{j+2}^z{\sigma }_{j+3}^z{\sigma }_{j+4}^{+}+h.c.)+\hfill \\ & \sum _{j=1}^{4L-1}({\Delta }_1{\sigma }_j^{+}{\sigma }_{j+1}^{+}+h.c)+\sum _{j=1}^{4L-4}U{n}_j{n}_{j+4}.\hfill \end{array}\end{array}$$

(8)

We see that in Equation (8), there appear explicitly the parity strings between sites that are not nearest neighbors; these strings encode the non-local character of the mapping.

Equation (8) can be represented as a matrix product operator (MPO), i.e., as a site-by-site decomposition of the Hamiltonian into the product of matrices containing operators that act only on one site:

$$\begin{array}{c}\hfill H=\sum _{\left\{w\right\},\left\{\sigma \right\},\left\{{\sigma }^{\prime }\right\}}{M}_{w_1}^{{\sigma }_1^{\prime },{\sigma }_1}{M}_{w_1,w_2}^{{\sigma }_1^{\prime },{\sigma }_1}\cdots M_{w_{L-1}}^{{\sigma }_L^{\prime },{\sigma }_L}|{\sigma }_1^{{}^{\prime }}{\sigma }_2^{{}^{\prime }}\cdots {\sigma }_L^{{}^{\prime }}\rangle \langle {\sigma }_L^{}{\sigma }_{L-1}^{}\cdots {\sigma }_1^{}|,\end{array}$$

(9)

with $M_{w_{k-1},w_k}^{{\sigma }_k^{{}^{\prime }},{\sigma }_k}$ a fourth rank tensor. The physical indices (${\sigma }_k^{{}^{\prime }},{\sigma }_k=1,2$) depend on the dimension of the single-site Hilbert spaces, while the tensor indices ($w_{k-1},w_k=1,\cdots ,14$) depend on the structure of Hamiltonian in Equation (8).

Not all sites have square MPOs. Indeed, the first and last tensors ($M_{w_1}^{{\sigma }_1^{{}^{\prime }},{\sigma }_1}$, $M_{w_{L-1}}^{{\sigma }_L^{{}^{\prime }},{\sigma }_L}$) are of rank three and are crucial to reproduce the full Hamiltonian, Equation (8). The analytic expressions of local tensors, the corresponding graph representations of the tensor structure, and other technical details are reported and discussed in Appendix A.

Once the MPO structure of the Hamiltonian is obtained, simulations with the DMRG code are performed. The results of the numerical simulations and the emerging physical properties are discussed in the following subsection.

3.2. Numerical Results

Topological phases are usually identified according to the following criteria: (i) evidence of degenerate ground states with different parities, (ii) evidence of nonlocal fermionic correlations between the edges, (iii) a degenerate entanglement spectrum [20,22]. Here, we focus on properties (ii) and (iii).

We consider a strip of length $L=48$, width $N=4$, and set the parameters of the model at the point of exact topological degeneracy, $w=0.1$ and $v=1$, in the phase diagram reported in Figure 1b.

In Figure 3, we report the non-local fermion correlations as computed on the many-body ground state $|GS$〉. We use the shorthand notation <$c_{m,l}^{†}{c}_{m^{\prime },l^{\prime }}$> in place of 〈$GS\left|c_{m,l}^{†}{c}_{m^{\prime },l^{\prime }}\right|GS$〉. In Figure 3a,b, we report, respectively, <$c_{1,1}^{†}{c}_{1,l}$> and <$c_{L,4}^{†}{c}_{L,l}$> which are the fermionic correlations along the short dimension of the strip, while Figure 3c–f show, respectively, the correlations of the four corners along the length L of the strip: <$c_{1,1}^{†}{c}_{m,1}$>, <$c_{L,1}^{†}{c}_{m,1}$>, <$c_{1,4}^{†}{c}_{m,4}$>, and <$c_{L,4}^{†}{c}_{m,4}$>. Here, $l=1,2,3,4$ is the chain index while $m=1$, ⋯, L denotes the site index. We also draw a comparison between three cases corresponding to three different values of the interaction strength ($U=0.25$, $0.5$, 1) and the non-interacting case $U=0$.

The presence of robust fermionic modes localized at the edge of the strip is clearly visible in Figure 3. In particular, in Figure 3a, the correlations attain finite values at corners $(1,1)$–$(1,4)$ and in Figure 3b, at corners $(L,1)$–$(L,4)$, for all the three reported values of the interaction strength. The correlations vanish elsewhere in the strip geometry. Moreover, Figure 3c–f shows that the fermionic correlations rapidly decay to zero when one of the four corners of the strip is fixed and the site index m is increased from 1 to L. The aforementioned analysis suggests that, for increasing the width of the strip $N\gg 4$, these fermionic modes are expected to split into corner Majorana modes, correlating at the corners $(1,1)$–$(1,N)$ and $(L,1)$–$(L,N)$.

Let us now discuss the behavior of the entanglement spectrum. Denote by $|\Psi \rangle ={\sum }_{\alpha =1}^D{\lambda }_{\alpha }|{\varphi }_{\alpha }{\rangle }_A{\varphi }_{\alpha }{\rangle }_B$ the many-body ground state with respect to some bipartition of the system with ${\lambda }_{\alpha }$ real numbers and $D=d^{4L}$ the Hilbert space dimension. Here d denotes the dimension of the single-site Hilbert space, and the set of real eigenvalues ${\left\{{\lambda }_{\alpha }\right\}}_{\alpha =1}^D$ forms the entanglement spectrum. In the DMRG procedure, this expression can be replaced by $|\Psi \rangle ={\sum }_{\alpha =1}^m{\lambda }_{\alpha }|{\varphi }_{\alpha }{\rangle }_A{\varphi }_{\alpha }{\rangle }_B$, with the entanglement spectrum reducing to ${\left\{{\lambda }_{\alpha }\right\}}_{\alpha =1}^m$, and the dimension D replaced by a fixed bond dimension m. In Figure 3 and Figure 4, the bond dimension is fixed at the value $m=200$.

In the topological phase, the entanglement spectrum is expected to be two-fold degenerate with respect to the parity sector, this feature being a precursor of the zero energy Majorana edge states. Indeed, the presence of the Majorana edge states implies the occupation of nonlocal fermionic modes, allowing such degeneracy with respect to parity. Figure 4a–c reports the behavior of the entanglement spectrum as a function of the eigenvalue index $\alpha$ for three different values of the interaction strength U. All panels show a degenerate entanglement spectrum, a typical signature of an ordered phase.

Summing up, the findings reported in Figure 3 and Figure 4 show that Coulomb repulsive interactions do not affect significantly the topological phases of the model, featuring robust edge fermionic modes and a degenerate entanglement spectrum also in the strongly interacting regime.

4. Conclusions

We analyzed the effect of nearest-neighbor Coulomb repulsive interactions on a Majorana BBH model restricted to a strip of $N=4$ chains of length L with interactions along the chains. The simulations, performed using a DMRG procedure, provide some interesting insights on the robustness of topological phases against repulsive interactions. In particular, we observed robust fermionic modes localized at the edges of the strip, together with a degenerate entanglement spectrum. Both observations point at the existence of a robust topological order, persistent even in the regime of strong interactions.

These results can be considered preliminary to the investigation of the robustness of Majorana corner states in an interacting two-dimensional MBBH model, which can be simulated by a tree tensor network approach in the limit of large lattices.

In order to discriminate unambiguously the topological order from other types of order, we will need to go beyond the entanglement spectrum, since the latter does not discriminate between different types of order, e.g., topological order and orders associated to spontaneous symmetry breaking. We will thus need to consider entanglement measures able to quantify the nonlocal correlations between the edges and, in particular, the long-distance topological entanglement between corner Majorana modes. In fact, such a measure, the topological squashed entanglement, exists and was recently introduced and successfully applied [44] to the unambiguous characterization of topological order in some basic models of topological superconductivity, including the Kitaev chain, the two-leg Kitaev ladder [45], and the Kitaev tie [46]. We plan to report in the near future the results of a similar investigation along the same lines for the interacting Majorana BBH model.