An exactly solvable model of 2D Majorana zero modes

The Kitaev's Majorana Model is an exactly solvable model of p-wave superconductor with localized Majorana zero modes in 1D quantum wire. For the 2D case, the general theory of Majorana zero modes near edges and in the core of vortices was set up in a paper by N.Read and D.Green in 2000.

Is there an exactly solvable model of 2D chiral p-wave superconductor with Majorana zero modes? (I mean, the explicit form of the $2^M$ degenerate ground states in the presence of $2M$ vortices should be analytically solvable. I'm eager to see how braiding of Majorana fermions lead to Non-Abelian unitary transformation of ground states in a manifest manner.)

PS: If anyone working in this area tells me he/she has never heard of such an exactly solvable model, then it would be enough for me--I recently constructed a fancy one, just come here to make sure it hasn't been studied before. I promise to come back to answer my own question as soon as my paper is published:)

Well it depends on what you mean by "exactly solvable". Many people would consider any quadratic Hamiltonian(i.e. non-interacting) exactly solvable, so the BCS Hamiltonian of a p+ip superconductor already falls into this category. One can demand more -- like at a special point of Kitaev chain ($t=\Delta, \mu=0$), the Hamiltonian consists of commuting terms. This is in my opinion the strongest form of exact solubility. However we believe that this is not possible for chiral p+ip superconductors. There are other exactly solvable models (in the sense of commuting terms), in which certain types of vortices/fluxes bind Majorana zero modes.

• Hi, Thanks! Could you please give me some references for your last sentence? I'm very interested in such models! – Lagrenge Jan 17 '17 at 12:17
• – Meng Cheng Jan 17 '17 at 14:10
• And I should also add Kitaev's famous honeycomb lattice model. But the exact solvability lies in the fact that the spin model can be mapped to quadratic fermion Hamiltonian, so essentially it is just p+ip superconductor, in this regard. – Meng Cheng Jan 17 '17 at 14:12

OK, as my paper is already on arxiv, it's time to come back to answer my own question. In our paper we found a special point of the Read-Green's mean-field Hamiltonian at which all ground state wavefunctions can be obtained explicitly even in translation non-invariant geometries~(e.g. with boundaries and vortices), analogous to the special point $\Delta=t,\mu=0$ in Kitaev's Majorana chain.

As Meng has mentioned, it is difficult to construct a chiral p+ip superconducting Hamiltonian as a sum of commuting local operators. However, our exact solutions does not depend on commuting operators. It is based on the "annihilator construction" that lays the foundation of frustration-free Hamiltonians: if a many-body Hamiltonian can be written in the form $$\label{eq:anni} \hat{H}=\sum_{i}A^\dagger_iA_i,\tag{1}$$ where $i$ is a collective spatial index, and if there exist a state $|G\rangle$ satisfying $A_i|G\rangle=0$ for $\forall i$, then $|G\rangle$ must be a ground state of $\hat{H}$ and consequently all ground states of $\hat{H}$ should be annihilated by $A_i$ for $\forall i$.

Let's start from the Read-Green's mean-field Hamiltonian $$\label{eq:K_2Dpp} \hat{K}=\int_{S}\left[ \frac{\nabla \psi_{z}^{\dagger }\cdot \nabla \psi _{z}}{2m}-(\Delta \psi_{z}\partial_{\bar{z}}\psi_{z}+\mathrm{H.c.})+\mu \psi_{z}\psi_{z}^{\dagger }\right] d^{2}z, \tag{2}$$ where $S$ denotes a region in the 2D plane with complex coordinates $z=x+iy$, ${\bar z}=x-iy$, $\partial_{z}=(\partial_{x}-i\partial_{y})/2$, $d^{2}z=dxdy$, and $\psi_{z}$ is the fermionic annihilation operator at position $z$. At $\mu=m\Delta^2/2$, we can factorize $\hat{K}$ as \begin{eqnarray}\label{eq:K_bb} \hat{K} =\int_{S}\frac{1}{2m}(2\partial_{z}\psi_{z}^{\dagger }-m\Delta^*\psi _{z})(2\partial_{\bar{z}}\psi_{z}-m\Delta\psi_{z}^{\dagger })d^{2}z, \tag{3} \end{eqnarray} where $\partial S$ denotes the boundary of $S$ and an irrelevant boundary term is neglected. In the following we use natural units $2m=m\Delta =1$ for simplicity.

Eq.~(\ref{eq:K_bb}) has essentially the same form as Eq.~\eqref{eq:anni}, so the remaining question is whether we can find states $|G\rangle$ satisfying $(2\partial_{\bar{z}}\psi_{z}-\psi_{z}^{\dagger })|G\rangle=0$. Let's try an ansatz state $$\label{eq:Gg} |G\rangle=|G_{[g]}\rangle\equiv\exp \left[ \frac{1}{2}\int_{S}g(z,z^{\prime })\psi _{z}^{\dagger }\psi_{z^{\prime }}^{\dagger }d^{2}zd^{2}z^{\prime }\right] |0\rangle , \tag{4}$$ where the two-particle wavefunction $g(z,z^{\prime })$~(to be determined) satisfies $g(z,z^{\prime })=-g(z^{\prime },z)$. In order to have $(2\partial_{\bar{z}}\psi_{z}-\psi_{z}^{\dagger })|G\rangle=0$, we need \begin{eqnarray} 0&=&(2\partial_{\bar{z}}\psi_{z}-\psi_{z}^{\dagger })|G_{[g]}\rangle= (2\partial_{\bar{z}}\frac{\delta}{\delta\psi^\dagger_{z}}-\psi_{z}^{\dagger })|G_{[g]}\rangle\nonumber\\ &=&[\int_S2\partial_{\bar{z}}g(z,z')\psi^\dagger_{z'}d^2z'-\psi^\dagger_z]|G_{[g]}\rangle.\tag{5} \end{eqnarray} Therefore, $g(z,z')$ has to satisfy the differential equation $$\label{eq:gzz} 2\partial_{\bar{z}}g(z,z') = \delta^{2}(z-z'). \tag{6}$$ The solution to Eq.~\eqref{eq:gzz} in an unbounded 2D plane takes a very simple form \begin{eqnarray}\label{eq:G_TIG1} g(z,z')=\frac{1}{2\pi(z-z')},\tag{7} \end{eqnarray} in virtue of the identity $$\label{eq:deltafunction} \frac{1}{\pi}\partial_{\bar{z}}\frac{1}{z-z'}=\delta^2(z-z').\tag{8}$$ It is easy to prove that Eq.~\eqref{eq:deltafunction} is the unique solution to Eq.~\eqref{eq:gzz} that vanishes as $z\to \infty$~($z'$ fixed), using Liouville's theorem. So the unique ground state in an unbounded 2D plane is $$|G\rangle =\exp \left[ \frac{1}{4\pi}\int_{S}\frac{1}{z-z^{\prime }}\psi _{z}^{\dagger }\psi_{z^{\prime }}^{\dagger }d^{2}zd^{2}z^{\prime }\right] |0\rangle.\tag{9}$$

Now let's try to obtain all degenerate ground states in an unbounded 2D plane with 4 vortices at positions $\eta_j,j=1,2,3,4$~(This is the simplest case where non-Abelian statistics happens). The Hamiltonian of the system is formally the same as Eq.~\eqref{eq:K_bb}, with the understanding that fermion operators $\psi_z,\psi^\dagger_z$ are multivalued functions of $z$--they acquire a minus sign when $z$ winds around any vortex $\eta_j$~(Read and Green's gauge convention). We can still use the ansatz state in Eq.~\eqref{eq:Gg} where $g(z,z')$ satisfies the same equation Eq.~\eqref{eq:gzz}, however, in order for the integral in Eq.~\eqref{eq:Gg} to be well-defined, $g(z,z')$ must have the same multivalued behavior with $\psi^\dagger_z$ when considered as a function of $z$ for any fixed $z'$. Thus we're searching for solutions to Eq.~\eqref{eq:gzz} with four branch points at $\eta_1,\eta_2,\eta_3,\eta_4$, in both $z$ and $z'$. A good candidate is \begin{eqnarray} \label{eq:g_1234} g_{12,34}(z,z^{\prime })&=&\frac{1}{4\pi(z-z^{\prime })}\left[\sqrt{\frac{(z-\eta_1)(z-\eta_2)(z^{\prime }-\eta_3)(z^{\prime }-\eta_4)}{(z^{\prime }-\eta_1)(z^{\prime }-\eta_2)(z-\eta_3)(z-\eta_4)}}+(z\leftrightarrow z^{\prime })\vphantom{\frac{1}{2}}\right].\tag{10} \end{eqnarray} One can easily verify that $g_{12,34}(z,z^{\prime })$ is a solution to Eq.~\eqref{eq:gzz}, since it is essentially the same as the vortex free case Eq.~\eqref{eq:G_TIG1} in the sense that it is also a meromorphic function of $z,z'$ having a first order pole at $z=z'$ with same residue $1/2\pi$, and has the same asymptotic behavior as $z\to \infty$. Plus, it has the four branch points at required positions $\eta_1,\eta_2,\eta_3,\eta_4$. This gives us one of the ground states $|G_{12,34}\rangle=|G_{[g_{12,34}]}\rangle$. By permuting the indices $1,2,3,4$ we get two other solutions $g_{13,24}(z,z')$ and $g_{14,23}(z,z')$, however, one can prove the following relations~(where $\eta_{ij}=\eta_i-\eta_j$,) \begin{eqnarray} \eta_{12}\eta_{34}&=&\eta_{13}\eta_{24}-\eta_{14}\eta_{23},\nonumber\\ \eta_{12}\eta_{34}~g_{12,34}&=&\eta_{13}\eta_{24}~g_{13,24}-\eta_{14}\eta_{23}~g_{14,23},\nonumber\\ \eta_{12}\eta_{34}|G_{12,34}\rangle&=&\eta_{13}\eta_{24}|G_{13,24}\rangle-\eta_{14}\eta_{23}|G_{14,23}\rangle,\tag{11} \end{eqnarray} i.e. the three different solutions are linearly dependent. Therefore, the ground state subspace with even fermion parity is two-dimensional. This is consistent with the fact that with $4$ vortices there are $4$ independent Majorana edge modes $\gamma_{1}\ldots \gamma_{4}$, which could be combined to $2$ Dirac fermion modes $\chi_1=(\gamma_{0,1}+i\gamma_{0,2})/2,~\chi_2=(\gamma_{0,3}+i\gamma_{0,4})/2$, resulting in $4$ degenerate ground states, $2$ of which have even fermion parity. We can construct the occupation number basis for the nonlocal fermions $\chi_{1},\chi_2$ \begin{eqnarray}\label{eq:basis0011} |00\rangle&=&\sqrt{N_{00}}(\lambda_{00}|G_{13,24}\rangle+\bar{\lambda}_{00}|G_{14,23}\rangle),\nonumber\\ |11\rangle&=&\sqrt{N_{11}}(\lambda_{11}|G_{13,24}\rangle+\bar{\lambda}_{11}|G_{14,23}\rangle),\tag{12} \end{eqnarray} where $N_{00(11)}=\sqrt{\eta_{13}\eta_{24}}\pm\sqrt{\eta_{14}\eta_{23}}$, $\lambda_{00}=\sqrt{\eta_{13}\eta_{24}}/N_{00}$, $\lambda_{11}=\sqrt{\eta_{13}\eta_{24}}/N_{11}$, $\bar{\lambda}_{00}=1-\lambda_{00}$ and $\bar{\lambda}_{11}=1-\lambda_{11}$. The two states satisfy $\chi_1|00\rangle=\chi_2|00\rangle=0,~|11\rangle=\chi^\dagger_1\chi^\dagger_2|00\rangle$. Therefore, they are orthogonal and have equal norm. It is interesting to point out that the non-Abelian statistics of vortices can be directly calculated from these wavefunctions, by braiding the vortex positions $\eta_j$ and study the multivalued behavior of the wavefunctions when analytically continued. For example, if we do a counterclockwise exchange of vortices $\eta_2,\eta_3$, the states transform to \begin{eqnarray}\label{eq:B_23_f} |00\rangle\to\frac{e^{\frac{\pi}{4}i}}{\sqrt{2}}(|00\rangle-i|11\rangle),~~|11\rangle\to\frac{e^{\frac{\pi}{4}i}}{\sqrt{2}}(|11\rangle-i|00\rangle).\tag{13} \end{eqnarray} The result agrees with the previous approach based on universal arguments~(up to an Abelian phase factor).

An interesting application of the special point presented here is that we can use it to build a particle number conserving, interacting Hamiltonian whose ground states are exactly projected mean-field states (see my paper and an earlier paper for details). For example, for the Hamiltonian \begin{eqnarray}\label{eq:H_pipsimple} \hat{H}=\int_S\nabla\psi^\dagger_{z}\cdot\nabla \psi_{z}d^2 z+4\int_S e^{-\lambda|z-z'|}\psi^{\dagger}_{z}(\partial_{z^{\prime }}\psi^\dagger_{z^{\prime }})\psi_{z^{\prime }}\partial_{\bar{z}}\psi_z d^2z d^2z^{\prime },\tag{14} \end{eqnarray} the ground states in the $2N$-particle sector are $$|G_{2N}\rangle=\hat{P}_{2N}|G_{[g]}\rangle= \left[ \frac{1}{2}\int_{S}g(z,z')\psi_{z}^{\dagger }\psi_{z^{\prime }}^{\dagger }d^{2}zd^{2}z^{\prime }\right]^N |0\rangle,\tag{15}$$ where $g(z,z')$ is the same two-particle wavefunction as before, given in Eq.~\eqref{eq:G_TIG1} for the unbounded 2D plane and in Eq.~\eqref{eq:g_1234} for the case with vortices.