Quantum Monte Carlo simulation of a two-dimensional Majorana lattice model

TL;DR

This paper employs quantum Monte Carlo simulations to investigate a two-dimensional Majorana fermion lattice, revealing insights into time-reversal symmetry breaking in topological superconductors.

Contribution

It introduces an ab-initio quantum Monte Carlo approach for Majorana systems with a semi-positive Pfaffian measure, advancing simulation techniques for topological phases.

Findings

Spontaneous time-reversal symmetry breaking observed at finite temperature.

Successful implementation of quantum Monte Carlo for Majorana fermions with a semi-positive measure.

Provides a new computational framework for studying topological superconductors.

Abstract

We study interacting Majorana fermions in two dimensions as a low-energy effective model of a vortex lattice in two-dimensional time-reversal-invariant topological superconductors. For that purpose, we implement ab-initio quantum Monte Carlo simulation to the Majorana fermion system in which the path-integral measure is given by a semi-positive Pfaffian. We discuss spontaneous breaking of time-reversal symmetry at finite temperature.

Figures

Equations

$${\mathcal{T}}^{-1}{\psi }_{\uparrow }\mathcal{T}={\psi }_{\downarrow },{\mathcal{T}}^{-1}{\psi }_{\downarrow }\mathcal{T}=-{\psi }_{\uparrow }.$$

$$\begin{array}{rl}\mathcal{H}=& \sum _{\mu }it\{{\psi }_{\uparrow }(x){\psi }_{\uparrow }(x+\hat{\mu })-{\psi }_{\downarrow }(x){\psi }_{\downarrow }(x+\hat{\mu })\}\\ & -\sum _{\mu }g{\psi }_{\uparrow }(x){\psi }_{\uparrow }(x+\hat{\mu }){\psi }_{\downarrow }(x){\psi }_{\downarrow }(x+\hat{\mu }),\end{array}$$

$$\begin{array}{rl}\mathcal{H}=& \sum _{\mu }it\{\psi (x)\psi (x+\hat{\mu })+{\psi }^{\ast }(x){\psi }^{\ast }(x+\hat{\mu })\}\\ & -\sum _{\mu }g\left(n(x)-\frac{1}{2}\right)\left(n(x+\hat{\mu })-\frac{1}{2}\right),\end{array}$$

$$\begin{array}{rl}\mathcal{H}=& -\sum _{\mu }t\{{\psi }^{\ast }(x)\psi (x+\hat{\mu })+{\psi }^{\ast }(x)\psi (x+\hat{\mu })\}\\ & +\sum _{\mu }g\left(n(x)-\frac{1}{2}\right)\left(n(x+\hat{\mu })-\frac{1}{2}\right),\end{array}$$

$$Z=\int D{\psi }_{\uparrow }D{\psi }_{\downarrow }{e}^{-S},$$

$$S=\int d\tau\sum_{x}\bigg{[}\psi_{\uparrow}\frac{\partial}{\partial\tau}\psi_{\uparrow}+\psi_{\downarrow}\frac{\partial}{\partial\tau}\psi_{\downarrow}+\mathcal{H}\bigg{]}.$$

$$Z=\int D{\psi }_{\uparrow }D{\psi }_{\downarrow }DC{e}^{-S^{\prime }-S_C},$$

$$\begin{split}S^{\prime}=&\int d\tau\sum_{x}\bigg{[}\psi_{\uparrow}(x)\frac{\partial}{\partial\tau}\psi_{\uparrow}(x)+\psi_{\downarrow}(x)\frac{\partial}{\partial\tau}\psi_{\downarrow}(x)\\ &+\sum_{\mu}it\{\psi_{\uparrow}(x)\psi_{\uparrow}(x+\hat{\mu})-\psi_{\downarrow}(x)\psi_{\downarrow}(x+\hat{\mu})\}\\ &-\sum_{\mu}iC_{\mu}(x)\{\psi_{\uparrow}(x)\psi_{\downarrow}(x)+\psi_{\uparrow}(x+\hat{\mu})\psi_{\downarrow}(x+\hat{\mu})\}\bigg{]},\end{split}$$

$$S_C=\int d\tau \sum _x\sum _{\mu }\frac{1}{2g}{C}_{\mu }^2(x).$$

$$\frac{1}{2g}C=⟨i{\psi }_{\uparrow }{\psi }_{\downarrow }⟩=\mathrm{Tr}\left(\frac{{\sigma }_2}{2}\frac{1}{M}\right).$$

$$\begin{array}{rl}{S}^{\prime }& \equiv \frac{1}{2}{\Psi }^{\top }M\Psi \\ & =\frac{1}{2}\left(\begin{array}{cc}{\psi }_{\uparrow }^{\top }& {\psi }_{\downarrow }^{\top }\end{array}\right)\left(\begin{array}{cc}{M}_{\uparrow \uparrow }& M_{\uparrow \downarrow }\\ M_{\downarrow \uparrow }& M_{\downarrow \downarrow }\end{array}\right)\left(\begin{array}{c}{\psi }_{\uparrow }\\ {\psi }_{\downarrow }\end{array}\right),\end{array}$$

$$M_{\uparrow \uparrow }$$

$$M_{\downarrow \downarrow }$$

$$M_{\uparrow \downarrow }$$

$$Z=\int D{\psi }_{\uparrow }D{\psi }_{\downarrow }DA{e}^{-S^{\prime }-S_A},$$

$$\displaystyle\begin{split}S^{\prime}=&\int d\tau\sum_{x}\bigg{[}\psi_{\uparrow}(x)\frac{\partial}{\partial\tau}\psi_{\uparrow}(x)+\psi_{\downarrow}(x)\frac{\partial}{\partial\tau}\psi_{\downarrow}(x)\\ &+\sum_{\mu}[\{it+A_{\mu}(x)\}\psi_{\uparrow}(x)\psi_{\uparrow}(x+\hat{\mu})\\ &+\{-it+A_{\mu}(x)\}\psi_{\downarrow}(x)\psi_{\downarrow}(x+\hat{\mu})]\bigg{]},\end{split}$$

$$S_A=\int d\tau \sum _x\sum _{\mu }\frac{1}{2g}{A}_{\mu }^2(x).$$

$$\frac{\partial }{\partial \tau }\Psi (\tau )\to {\hat{\partial }}_{\tau }\Psi (\tau )\equiv \frac{1}{2\delta \tau }\{\Psi (\tau +\delta \tau )-\Psi (\tau -\delta \tau )\}.$$

$$\frac{1}{M({\omega }_n,p=0)}=\frac{1}{2}\frac{\delta \tau }{i\sin ({\omega }_n\delta \tau )},$$

$$\begin{array}{rl}\frac{\partial }{\partial \tau }\Psi (\tau )& \to \left({\hat{\partial }}_{\tau }-{\sigma }_2\frac{{\hat{\Delta }}_{\tau }}{2}\right)\Psi (\tau )\\ & \equiv \frac{1}{2\delta \tau }\{\Psi (\tau +\delta \tau )-\Psi (\tau -\delta \tau )\}\\ & -{\sigma }_2\frac{1}{2\delta \tau }\{\Psi (\tau +\delta \tau )+\Psi (\tau -\delta \tau )-2\Psi (\tau )\}.\end{array}$$

$$\begin{array}{rl}& \frac{1}{M({\omega }_n,p=0)}\\ & =\frac{1}{2}\frac{\delta \tau }{i\sin ({\omega }_n\delta \tau )+{\sigma }_2\{1-\cos ({\omega }_n\delta \tau )\}}\\ & =-\frac{\delta \tau }{4}\frac{i\sin ({\omega }_n\delta \tau )-{\sigma }_2\{1-\cos ({\omega }_n\delta \tau )\}}{1-\cos ({\omega }_n\delta \tau )}.\end{array}$$

$$Z=\int DA\mathrm{pf}M{e}^{-S_A},$$

$$M=\left(\begin{array}{cc}{M}_{\uparrow \uparrow }& M_{\uparrow \downarrow }\\ M_{\downarrow \uparrow }& M_{\downarrow \downarrow }\end{array}\right),$$

$$\begin{array}{rl}{M}_{\uparrow \uparrow }=2{\hat{\partial }}_{\tau }& +\sum _{\mu }[\{it+A_{\mu }(x)\}{P}_{+\mu }\\ & -\{it+A_{\mu }(x-\hat{\mu })\}{P}_{-\mu }],\end{array}$$

$$\begin{array}{rl}{M}_{\downarrow \downarrow }=2{\hat{\partial }}_{\tau }& +\sum _{\mu }[\{-it+A_{\mu }(x)\}{P}_{+\mu }\\ & -\{-it+A_{\mu }(x-\hat{\mu })\}{P}_{-\mu }],\end{array}$$

$$M_{\uparrow \downarrow }=-M_{\downarrow \uparrow }=i{\hat{\Delta }}_{\tau }.$$

$$\begin{array}{rl}Z& =\int DA(\det M{M}^{†}{)}^{\frac{1}{4}}{e}^{-S_A}\\ & =\int D{\Phi }^{\ast }D\Phi DA{e}^{-\frac{1}{4}{\Phi }^{†}(M{M}^{†}{)}^{-1}\Phi -S_A}.\end{array}$$

Full text

I Introduction

Realization of Majorana fermions, which are exotic fermions that behave as their own anti-particles (holes), has been intensively discussed in condensed matter systems such as fractional quantum Hall (Moore-Read) states Moore and Read (1991), topological superconductors (superfluids) Read and Green (2000); Ivanov (2001); Schnyder et al. (2008); Qi et al. (2009), and surface states of topological insulators under the proximity effect between trivial superconductors Fu and Kane (2008) (For reviews, see e.g., Refs. Hasan and Kane (2010); Qi and Zhang (2011); Mizushima et al. (2016)). In these materials, Majorana fermions emerge as collective excitations localized at surfaces or defects (quantum vortices).

Majorana zero modes or Majorana bound states localized at quantum vortices cost zero energy, and lead to large ground-state degeneracy Moore and Read (1991); Read and Green (2000); Ivanov (2001). The degenerate ground-state wave function has the non-abelian Berry phase. When the position of two vortices are adiabatically rearranged, it results in the so-called braiding statistics Nayak and Wilczek (1996), which is considered to be potentially useful for quantum computation Kitaev (2003, 2001); Nayak et al. (2008).

As the number of vortices increases, the vortices form lattice structure such as the Abrikosov lattice. Low-energy dynamics of such a vortex lattice in topological superconductors is effectively described as a many-body system of Majorana zero modes. The system obeys the lattice Hamiltonian describing the inter-vortex tunnelings and interactions of them, which is similar to the Hubbard Hamiltonian. Revealing the ground-state property of such a many-body Majorana system is important not only for applications of the vortex lattices to quantum computation but also for understanding new phases of strongly correlated electron systems where the fundamental degrees of freedom are Majorana zero modes.

However, previous studies have mostly focused on one-dimensional Majorana chains, which can be solved by field-theoretical techniques such as bosonization, or by numerical techniques such as exact diagonalization and density matrix renormalization group Grover et al. (2014); Rahmani et al. (2015a, b); Chiu et al. (2015); Katsura et al. (2015). It is hard to implement these techniques to higher-dimensional systems, particularly to large two-dimensional systems although two dimensions is often critical dimension in competition between disorder by thermal/quantum fluctuation and order by mean-field dynamics. Thus ab-initio simulation is required to precisely determine the phase structure.

In this paper we study non-perturbative dynamics of interacting Majorana fermions by utilizing the quantum Monte Carlo simulations which are based on the path-integral formalism of lattice field theories. We consider an effective Majorana lattice model describing the low-energy behavior of the vortex lattice in topological superconductors with time-reversal symmetry. By applying the quantum Monte Carlo method to the effective model, we discuss spontaneous breaking of time-reversal symmetry at finite temperature. Time-reversal symmetry protects Majorana zero modes against opening a gap, and gapless (free) Majorana fermion states are realized in symmetric phase. On the other hand in broken phase, double Majorana fermions form a pair and behave as a single spinless Dirac fermion with dynamically generated mass. We note that the understanding of time-reversal symmetry and its spontaneous breaking is important for condensed matter realization of supersymmetry Grover et al. (2014); Rahmani et al. (2015a, b).

The rest of the paper is organized as follows. In Sec. II, we introduce the Hamiltonian and the path integral of the effective Majorana lattice model. By using the effective model, we first discuss the phase structure on the basis of the mean-field theory in Sec. III, and then show the result of quantum Monte Carlo simulation in Sec. IV. Finally we summarize this paper in Sec. V.

II Model

Let us introduce a low-energy effective model describing quantum vortex lattice in two-dimensional time-reversal-invariant topological superconductors. A quantum vortex supports a pair of Majorana zero modes $\psi_{\uparrow,\downarrow}$ with spin up and down Qi et al. (2009), which are related by time-reversal symmetry as

[TABLE]

Time-reversal symmetry is anti-unitary operation, namely, ${\cal T}^{-1}i{\cal T}=-i$.

We consider the model Hamiltonian

[TABLE]

where $\hat{\mu}$ is the unit lattice vector in the $\mu$ direction ($\mu=1$, 2). The first term describes the inter-vortex hopping and the second term describes the nearest-neighbor interaction between the spin-up and spin-down Majorana fermions. Both the hopping term and the interaction term keep time-reversal symmetry. The unique local term $i\psi_{\uparrow}(x)\psi_{\downarrow}(x)$, which can open a gap in energy spectrum, is prohibited by time-reversal symmetry since this term is time-reversal odd. We consider attractive interaction $g>0$, which favors pairing of Majorana fermions. For the attractive interaction, this model is free from the fermion sign problem for any lattice structure in any dimension, as shown in Sec. IV. We study this model on a two-dimensional square lattice.

A pair of two Majorana fermions can be rewritten by a single spinless Dirac fermion as $\psi(x)=\{\psi_{\uparrow}(x)+i\psi_{\downarrow}(x)\}/\sqrt{2}$. The Hamiltonian (2) is rewritten as

[TABLE]

where $n(x)=\psi^{*}(x)\psi(x)$. For a bipartite lattice such as a square lattice, the Hamiltonian can be further rewritten by the particle-hole transformation, $\psi(x)\to i^{x_{1}+x_{2}}\psi(x)$ on even sites $x$ and $\psi(x)\to(-i)^{x_{1}+x_{2}}\psi^{*}(x)$ on odd sites $x$. The model is equivalent to the spinless repulsive Hubbard model at half-fillingScalapino et al. (1984); Gubernatis et al. (1985)

[TABLE]

which was shown to be Majorana positive Li et al. (2015). The pairing term $i\psi_{\uparrow}(x)\psi_{\downarrow}(x)$ in the original Hamiltonian (2) corresponds to the Dirac mass term $\psi^{*}(x)\psi(x)-1/2$ in the Hamiltonian (3) and the staggered mass term $(-1)^{x_{1}+x_{2}}\left(\psi^{*}(x)\psi(x)-1/2\right)$ in the Hamiltonian (4).

We introduce the Euclidean path integral

[TABLE]

with the Euclidean action

[TABLE]

The integral of imaginary time $\tau$ is anti-periodic with the period $1/T$. In the path-integral formalism, the Majorana fermions $\psi_{\uparrow,\downarrow}$ are real Grassmann fields. In the following sections, we evaluate this path integral by using the mean-field theory and the quantum Monte Carlo simulation.

III Mean field

First we study this model in the mean-field approximation. By the Hubbard-Stratonovich transformation, we obtain

[TABLE]

with the bilinear fermion action

[TABLE]

and the auxiliary field action

[TABLE]

Assuming that the mean field is homogeneous and isotropic, $C_{\mu}(x)=C$, we obtain the gap equation

[TABLE]

The fermion matrix $M$ is defined by the fermion action in the matrix notation

[TABLE]

with the matrix elements

[TABLE]

and $P_{\pm\mu}\psi_{\uparrow,\downarrow}(x)=\psi_{\uparrow,\downarrow}(x\pm\hat{\mu})$.

Solving the gap equation (10), we obtain the pair condensate $\langle i\psi_{\uparrow}\psi_{\downarrow}\rangle$. As mentioned in Sec. II, the operator $i\psi_{\uparrow}\psi_{\downarrow}$ is time-reversal odd, so that the pair condensate $\langle i\psi_{\uparrow}\psi_{\downarrow}\rangle$ is an order parameter of spontaneous breaking of time-reversal symmetry. Its temperature dependence is shown in Fig. 1. At high temperature, time-reversal symmetry is preserved and there appear gapless Majorana fermions. At low temperature, time-reversal symmetry is broken and double Majorana fermions form a massive Dirac fermion. The finite-size-scaling analysis shows that this is a second-order phase transition. It corresponds to the phase transition to a density wave state in the Hamiltonian (4) Scalapino et al. (1984); Gubernatis et al. (1985).

IV Quantum Monte Carlo

Next, we study the phase structure in the quantum Monte Carlo method. We introduce another type of the Hubbard-Stratonovich transformation to satisfy the Majorana positivity condition Jaffe and Janssens (2016); Wei et al. (2016); Li et al. (2016); Hayata and Yamamoto (2017). After the Hubbard-Stratonovich transformation, the path integral is

[TABLE]

with

[TABLE]

In the path-integral Monte Carlo simulation, imaginary time $\tau$ is discretized. Naively, a temporal derivative is replaced by a central difference,

[TABLE]

However, the central difference has the doubling problem. To understand this problem, let us consider the non-interacting propagator in momentum space

[TABLE]

where $\omega_{n}$ is the fermionic Matsubara frequency. The propagator has two poles at $\omega_{n}\delta\tau\to 0$ and $\pi$ (mod $2\pi$). The pole at $\omega_{n}\delta\tau\to\pi$ is unphysical and thus called the doubler. Although one well-known solution in non-relativistic theory is the use of a forward or backward difference, it does not solve the problem now. In the case of the Majorana fermion, only the anti-symmetric part of the fermion matrix contributes to the path integral. Because the anti-symmetric part of a forward or backward difference is a central difference, the calculation with a forward or backward difference is equivalent to that with a central difference. To avoid the doubling problem, we discretize a temporal derivative as

[TABLE]

This is inspired by the Wilson fermion formalism in relativistic lattice field theory Wilson (1977). Since the second term in Eq. (20) is a higher order of $\delta\tau$, it does not contribute to the $\delta\tau\to 0$ limit. Now the non-interacting propagator becomes

[TABLE]

Only the physical pole at $\omega_{n}\delta\tau\to 0$ survives and the unphysical pole at $\omega_{n}\delta\tau\to\pi$ gets massive. This resolves the doubling problem. On behalf of it, time-reversal symmetry is explicitly broken at $\delta\tau\neq 0$. This is similar to the explicit chiral symmetry breaking of the Wilson fermion formalism Wilson (1977).

After the Majorana fermions are integrated out, the path integral (15) becomes

[TABLE]

where $\mathrm{pf}M$ is the Pfaffian of the anti-symmetric fermion matrix $M$. The fermion matrix is given by

[TABLE]

with

[TABLE]

This matrix satisfies the Majorana positivity condition Hayata and Yamamoto (2017). The Pfaffian $\mathrm{pf}M$ is real and semi-positive. Thus the path integral is free from the sign problem. We remark here that the sub-lattice symmetry is not essential for the Majorana positivity of the model Hamiltonian (2). It is still Majorana positive and does not suffer from the sign problem even for non-bipartite lattices such as a triangular lattice. (On the other hand, for the case of repulsive interaction $g<0$, the auxiliary field $A_{\mu}$ is replaced by $iA_{\mu}$. Then the fermion matrix no longer satisfies the Majorana positivity condition. The Pfaffian becomes indefinite and causes the sign problem.)

We performed the path-integral Monte Carlo simulation by the Hybrid Monte Carlo algorithm Duane et al. (1987). For applying the Hybrid Monte Carlo algorithm, we transformed the path integral (22) to

[TABLE]

The complex scalar field $\Phi$ is called the pseudo-fermion field. We note that the transformation (27) is invalid if the Pfaffian is indefinite. We extrapolated the infinite-volume limit $V\to\infty$ from three lattice volumes $V=8^{2}$, $12^{2}$, and $16^{2}$. We fixed the hopping parameter $t$ and the temporal discretization $\delta\tau$ at $t\delta\tau=0.1$. Taking the continuum limit $\delta\tau\to 0$ is a future work.

The simulation results are shown in Fig. 2. The non-interacting value $\langle i\psi_{\uparrow}\psi_{\downarrow}\rangle_{0}$, which is induced by the explicit symmetry breaking (20), is subtracted from the Monte Carlo data $\langle i\psi_{\uparrow}\psi_{\downarrow}\rangle$. The results are qualitatively consistent with the mean-field results. The symmetry is broken at low temperature and restored at high temperature. However, the phase transition seems crossover, not second order. This is the artifact of the explicit symmetry breaking (20). To determine the order of the physical phase transition, we need to take the continuum limit $\delta\tau\to 0$ or formulate a discretization scheme to preserve the symmetry.

The phase diagram in the $(g,T)$ plane is shown in Fig. 3. The second-order phase transition line in the mean-field approximation is also shown. The Monte Carlo data and the mean-field results are qualitatively consistent. Note that the temperatures of the path-integral Monte Carlo simulations are nonzero and the zero-temperature limit is obtained by the extrapolation of them. The extrapolated data are shown at $T=0$ in Fig. 3.

V Summary

We have introduced the effective Majorana lattice model describing the low-energy physics of the vortex lattice in time-reversal-invariant topological superconductors. We have studied spontaneous breaking of time-reversal symmetry at finite temperature. To perform ab-initio analysis, we implemented the path-integral Monte Carlo simulation for Majorana fermion systems. For attractive interaction, our model is free from the sign problem for any lattice structures including non-bipartite lattices in any dimensions, and can be evaluated by using the standard importance sampling method. We have shown that the phase diagram of Majorana fermion systems can be investigated in the same way as that of complex fermion systems.

Our formulation is applicable to general Majorana fermion systems. For example, it can be applied to variants of the Kitaev spin model Kitaev (2006); Yao and Kivelson (2007); Mandal and Surendran (2009); Katsura et al. (2015) by expressing them with Majorana fermions through the Jordan-Wigner transformation Feng et al. (2007); Chen and Hu (2007); Nasu et al. (2014). It will be interesting to study the spin liquid state by using our method. Another application is the Majorana fermion in particle physics, such as neutrinos, supersymmetry, and quantum vortices in color superconductors Nishida (2010); Yasui et al. (2010); Eto et al. (2014).

Acknowledgements.

T. H. was supported by JSPS Grant-in-Aid for Scientific Research (Grant No. JP16J02240). A. Y. was supported by JSPS KAKENHI (Grant No. JP15K17624). The numerical simulations were carried out on SX-ACE in Osaka University.