Spontaneous breaking of time-reversal symmetry in the orbital channel for the boundary Majorana flat bands

We study the boundary Majorana modes for the single component p-wave weak topological superconductors or superfluids, which form zero energy flat bands protected by time-reversal symmetry in the orbital channel. However, due to the divergence of density of states, the band flatness of the edge Majorana modes is unstable under spontaneously generated spatial variations of Cooper pairing phases. Staggered current loops appear near the boundary and thus time-reversal symmetry is spontaneously broken in the orbital channel. This effect can appear in both condensed matter and ultra-cold atom systems.

The boundary Majorana zero modes also appear in unconventional superconductors in 2D and 3D on edges or surfaces with suitable orientations determined by Cooper pairing symmetries. This is because along the directions of incident and reflected wavevectors of the Bogoliubov quasi-particles, the pairing gap functions have opposite signs [33,34,35]. These unconventional superconductors can be viewed as weak topological states. They are topologically non-trivial (or trivial) along the directions perpendicular to boundaries with (or without) zero energy modes, respectively.
Interactions can significantly change topological properties of edge and surface in topological insulators and superconductors. For non-chiral systems, these boundary modes are low energy midgap states, thus they are more susceptible to interactions than the gapped bulk states. In the helical edge Luttinger liquids of 2D topological insulators, the edge magnetic fluctuations are much stronger than those in the bulk, which can lead to spontaneous time-reversal symmetry breaking and destroy edge modes [36,37,38]. For the zero energy boundary Majorana modes along the coupled chains of the p xtopological superconductors, it has been found that interactions can even generate gaps without breaking time-reversal symmetry in the orbital channel [39]. Similar effects are also found for the helical edge Majorana modes of time-reversal invariant topological pairing states [40,41]. In both cases, the new topological class under interactions is classified by a Z 8 periodicity.
In this article, we consider the most natural interaction effects in the boundary states in the p-wave weak topological superconductors or paired superfluids: the coupling between Cooper pairing phases and the zero energy Majorana modes. The degeneracy of these boundary modes is protected by time-reversal symmetry, but is vulnerable under spontaneous time-reversal symmetry breaking in the orbital channel. Spatial variations of phases of Cooper pairing order parameters induce bonding among these boundary modes, and thus lower the energy. Staggered current loops are generated near the boundary, which split the zero energy Majorana peaks. Due to the divergence of the surface density of states, this time-reversal symmetry breaking mechanism is robust.
The rest of this paper is organized as follows. In Sect. 2, we explain the mechanism of spontaneous time-reversal symmetry breaking in the orbital channel wtih coupled superconducting quantum wires through weak links. In Sect. 3, a 2D pwave weak topological superconductor is studied with the open boundary condition. The spontaneous staggered orbital current loops appear on the topological non-trivial boundaries, and the consequential splitting of the Majorana zero bias peaks are calculated self-consistently. Experimental realizations of the above effects are discussed in Sect. 4. The summary is presented in Sect. 5.

Bonding between boundary Majorana modes
We begin with a heuristic example of two parallel quantum wires along the x-direction. Each of them is a 1D topological superconductor of single component fermion with p xpairing symmetry. As explained in Ref. [5], zero energy Majorana modes exist near the two ends of each wire. For example, for the end at x = L of the i-th wire (i = 1, 2), the operator for the Majorana mode can be expressed as where θ i is the Cooper pairing phase in the i-th chain. We denote ∆ the magnitude of the bulk gap, E f and k f are the Fermi energy and Fermi momentum, respectively. In the case of ∆ ≪ E f , the zero mode wavefunction is approximated as where ∆ is the coherence length along the x-direction. Now let us connect two right ends with a weak link where w is the width of the link. We assume w ≪ ξ x such that the uniform distribution of the pairing phase within each wire is not affected by the link. By expressing .." represents for Bogoliubov eigenstates outside the gap ∆, we arrive at two different contributions for the Josephson couplings: the usual one through second-order perturbation process, and the fractionalized one through the zero Majorana mode as Lx−∆L dxu 0 (x)u 0 (x + ∆y). The fractionalized Josephson coupling of the weak link between two wires with the "head-to-tail" configuration [7,42] gives rise to iγ 1 γ 2 cos ∆θ 2 , where γ 1,2 are Majorana modes at two sides of the link. In contrast, for two parallel wires, the corresponding J 1 term in Eq. 4 is proportional to sin ∆θ 2 . This difference can be intuitively understood as illustrated in Fig. 1. We rotate the 2nd wire by 180 • , which corresponds to a phase shift of the 2nd wire by π due to the p x symmetry, and thus ∆θ is shifted to ∆θ + π. This changes the dependence of fractionalized Josephson coupling from cosine to sine in Eq. 4. As a result, the Majorana fermion bonding strength vanishes at ∆θ = 0 which reflects the fact that these Majorana zero modes are protected by the TR symmetry in the orbital channel. The bonding strength varies linearly with the phase difference as ∆θ → 0, thus Eq. 4 is minimized at nonzero values of ∆θ 0 satisfying sin ∆θ 0 2 = ± min( J 1 4J 0 , 1). The energy is gained from the bonding of Majorana fermions. For the general case that ∆θ 0 = ±π, the orbital channel time-reversal symmetry is spontaneously broken. However, in the case of J 1 4J 0 ≥ 1, ∆θ 0 = ±π, the Majorana fermion bonding strength is maximal. Because of the 2π periodicity of phase angles, this case is also TR invariant.
We can further extend Eq. 4 to a group of parallel wires whose ends at x = L are weakly connected, which forms 1D or 2D lattices of Majorana modes coupled to the superfluid phases. The effective Hamiltonian is where ii ′ refers to the nearest neighbor bonds. The J 0 term favors a global phase coherence with a constant phase value in each wire. However, this results in a flat band of zero energy Majorana fermions which is highly unstable.
In bipartite lattices, a staggered phase configuration θ i = (−) i θ ′ 0 can minimize the ground state energy by generating a uniform bonding strength, where the value of θ ′ 0 can be obtained self-consistently. For a 2D lattice of surface Majorana fermions, if the superfluid phase distribution of θ i forms a vortex, the plaquette in which the vortex core is located exhibit a Z 2 vortex for the Majorana fermion. The motion of the superfluid vortex introduces the dynamic Z 2 flux for the Majorana fermions.

Spontaneous staggered orbital currents near the boundary
Having presented the spontaneous orbital channel time-reversal symmetry breaking of Majorana fermions with weak links, now let us consider a 2D superconductor of a single component fermion with the p x pairing. The bulk of the system with p x -pairing maintains time-reversal symmetry in the orbital channel unlike the p x + ip y one in which time-reversal symmetry is broken. We use the tight-binding model of the following meanfield Bogoliubov-de Gennes (B-de G) Hamiltonian with the open boundary condition where ∆ i,i+êx = c i+êx c i . In conventional superconductors, the open boundary only affects the magnitude of the pairing order parameter not on the phase distribution. The p x -type superconductors are weakly topological in 2D, and as is well-known, if the phase distribution of ∆ i,i+êx is uniform such that time-reversal symmetry is maintained in the oribtal channel, a flat band of Majorana surface modes appear on the boundary along the y-direction. However, this Majorana band flatness is unstable when the spatial variations of pairing phases are considered. It spontaneously breaks time-reversal symmetry and generates finite bonding strengths among the Majorana modes. Consequently, the macroscopic degeneracy of the surface zero modes is removed and the total energy is lowered. To confirm this intuitive picture, we perform the self-consistent solutions to the B-de G equation Eq. 6 at zero temperature. The mean-field Hamiltonian is diagonalized to obtain the Bogoliubov eigenstates whose eigen-operators are defined through c i = n u n (i)γ n + v * n (i)γ † n , where (u n , v n ) T are the eigenvectors with positive energy eigenvalues E n > 0. The gap equations for ∆ i,i+êx read The current along the bond between i and i +ê a (a = x, y) is where f (E) is the Fermi distribution function. The spontaneous time-reversal symmetry breaking effect in the orbital channel appears near the boundary normal to the x-direction: circulation current loops as depicted in Fig. 2 a). As presented in Fig. 2 b), this edge induces non-uniform phase distribution of the pairing order parameter ∆ i,i+êx . However, it is non-singular which does not exhibit the vortex configuration and thus cannot give rise to the current loops by itself. Loosely speaking there are two different contributions to the current expression of Eq. 8: the Bogoliubov modes whose energies outside the bulk gap, and the in-gap states of Majorana modes. The former contribution can be captured by the Ginzburg-Landau formalism as spatial phase variations of ∆ i,i+êx , while the latter arises from Majorana modes cannot. The circulation pattern of these induced edge currents are staggered, which is natural, since in the background of superfluidity the overall vorticity should be neutral in the absence of external magnetic fields. The size of current loops should be at the order of coherence lengths which can be estimated as ξ a /a 0 ≈ 2t a /∆ with a = x, y where a 0 is the lattice constant. Currents decay exponentially along the x-direction at the scale of ξ x , and change directions along the y-direction at the scale of 2ξ y . Next we calculate the local density of states (LDOS) near the boundary as depicted in Fig. 3. The expression for LDOS reads For comparison, we first present the result without self-consistency in Fig. 3 a) by  Figure 3. The edge LDOS spectra for the system in Fig. 2 δ(ω ± E n ) in Eq. 9 are approximated by 1 π η/[(ω ± E n ) 2 + η 2 ] with the broadening parameter η = 0.04. a), b), c) and d) are the LDOS spectra at points 1, 2, 3, and 4 marked in Fig. 2 a), respectively. a) LDOS without self-consistency by fixing ∆ i,i+êx at a uniform bulk value. b), c) and d) are self-consistent results. The corresponding parameter values are presented in the caption of Fig. 2. fixing order parameters ∆ i,i+êx uniform with the bulk value. The central peak at zero energy represents the edge Majorana states. However, with self-consistency, spatial variations of ∆ i,i+êx couple to the edge Majorana modes. The zero energy peaks are split as depicted in Fig. 3 b), c) and d). Sites 2, 3 (Fig. 3 b) and c), respectively) are right on the boundary normal to x-direction, and thus their coherence peaks are strongly suppressed. Site 4 (Fig. 3 d) is relatively inside, and thus the midgap peaks are suppressed. Site 2 locates in the middle of current loops, and thus the splitting of the zero bias peak is the largest.
Due to the divergence of density of states at zero energy, the above time-reversal symmetry breaking mechanism in the orbital channal is general. It also applies to the 3D case with the p z -pairing. Then the xy-surface is non-trivial along which Majorana zero energy flat bands appear with the assumption of time-reversal symmetry. Again due to the same reasoning, this degeneracy will be lifted by spontaneous TR symmetry breaking. In this case, we expect that the 2D vortex-anti-vortex pattern in Fig. 2 will change to that of 3D closed vortex rings, which may further form a lattice structure parallel to the xy-surface. Further study of this lattice structure will be deferred to a later publication. Moreover, the above physics also applies in the continuum not just in the tight-binding model. Staggered vortices or closed vortex rings will appear near the edge or surfaces within the size of coherence lengths, which will also split the zero energy peaks of Majorana fermions.
Next we briefly discuss the boundary Majorana flat bands for the TR invariant spinful p x -wave Cooper pairing. For example, let us consider the pairing order parameter in the bulk as ∆ i,i+ex = i|∆|σ 2 σ · d. Without loss of generality, we can assume that the d-vector lies in the xy-plane with the azimuthal angle φ, i.e., d x + id y = |d|e iφ . This corresponds to that spin-up fermions are paired with spin-up, and spin-down ones are paired with spin-down, with the relation ∆ ↑↑ = |∆|e iφ and ∆ ↓↓ = |∆|e −iφ . Then the boundary Majorana modes consist two branches γ ↑ and γ ↓ forming Kramer pairs. The bonding among Majorana modes can occur within each branch of Kramer pair by developing currents as explained above. If a residue magnetic interaction exists, another way to remove the degeneracy is to develop magnetism with the order parameter iγ ↑ γ ↓ . The spin polarization for such an order is in the xy-plane. In general, these two different bonding mechanisms can coexist.

Discussions
The above real p-wave Cooper pairing of single component fermion can be realized in ultra-cold dipolar fermionic molecular systems. For example, the 40 K-87 Rb systems have been cooled down below quantum degeneracy [43,44]. Moreover, the chemically stable dipolar fermion molecules of 23 Na-40 K have also been laser-cooled [45]. If the dipole moments are aligned along the z-axis by external electric fields, the interaction between two dipole moments exhibits the d r 2 −3z 2 -type anisotropy. It means that the interaction is attractive if the displacement vector between two dipoles is along the z-axis, and repulsive if it lies in the xy-plane. Naturally, this leads to the p z -type Cooper pairing in the strong coupling limit. Even in the weak coupling limit, partial wave analysis also shows that the dominant pairing symmetry is the p z -type. Due to the anisotropy of the interaction, it slightly hybridizes with other odd partial wave channels as predicted in previous works [46,47,48]. For the two-component dipolar fermions, the triplet p z Cooper pairing has also been theoretically predicted [49,50]. The above predicted spontaneous time-reversal symmetry breaking effects in the orbital channel will appear on the surfaces perpendicular to the z-axis.
In the condensed matter systems, the boundary Majorana fermion has been detected in the 1D quantum wire with SO coupling [16,17,18,19,20]. A magnetic field is also applied in parallel to the wire to lift the degeneracy at zero momentum such that the system is effectively one component if the Fermi energy lies inside the gap. In this case, the major effect of the magnetic field is the Zeeman effect, and the time-reversal symmetry is not explicitly broken in the orbital channel. The proximity effect from the superconducting electrodes induces superconductivity in the 1D wire, which leads to the zero energy Majorana ending states. If we put an array of these 1D wires together and their ends are weakly connected, these Majorana modes will also spontaneously build up bonding strength by developing staggered orbital currents patterns, the consequential splitting of the zero bias peaks will be exhibited in the tunneling spectroscopy.

Conclusion
In summary, we have found that the zero energy boundary Majorana flat bands of the p-wave superconductors are unstable towards spontaneous time-reversal symmetry breaking in the orbital channel, which results from the coupling between Majorana modes and the Cooper pairing phases. The divergence of density of states at zero energy leads to spontaneous formation of staggered current loops, which generates bondings among the Majorana modes and lifts the degeneracy. This is a robust mechanism which appears even without introducing new instabilities.