Topological Larkin-Ovchinnikov phase and Majorana zero mode chain in bilayer superconducting topological insulator films

Topological superconductors possess a bulk superconducting gap and boundary gapless excitations, known as “Majorana fermion”. Search for new systems with topological superconductivity is of fundamental and application importance due to the potential application of Majorana fermions in topological quantum computation. Here we show that the Larkin-Ovchinnikov superconducting phase with a finite momentum pairing can emerge in a model of bilayer superconducting topological insulator films, in which superconductivity appears for both the top and bottom surface states, and can be topologically non-trivial. This “topological Larkin-Ovchinnikov phase” is induced by an in-plane magnetic field and possesses a Majorana mode chain along the edge perpendicular to the in-plane magnetic field direction due to its non-trivial Z2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Bbb Z}_2$$\end{document} topological nature. Our theoretical model can be naturally realized in superconductor/topological insulator sandwich structure or in Fe(Te, Se) film, a topological material with superconductivity, and thus provides a route to explore unconventional superconductivity in existing systems. Topological superconductors possess Majorana modes at the edge and are important in future quantum computational devices. The authors theoretically demonstrate the Larkin-Ovchinnikov superconducting phase can be topologically non-trivial in certain sandwich structure or iron-based superconductor films and possess a chain of Majorana modes.

M agnetism and superconductivity are two fundamental states of matter in condensed matter physics, and the interplay between them continues to bring us intriguing phenomena. Unlike the conventional Cooper pairs with zero momentum in the Bardeen-Cooper-Schrieffer (BCS) theory, magnetism can induce a superconducting (SCing) state with a finite momentum pairing. The pairing function of such a state can either carry a single finite momentum Q, known as Fulde-Ferrell (FF) state 1 , or multiple finite momenta Q i (i = 1,2,…), known as Larkin-Ovchinnikov (LO) state 2 . The FFLO state has been theoretically proposed [3][4][5] and experimentally tested in a variety of systems, including heavy fermion superconductors (SCs) 6,7 , organic SCs 8 , ultrathin crystalline Al films 9 , and even in cold atom systems 10 . Another recent development is to realize topological SCs by integrating magnetism, spin-orbit coupling and superconductivity into one hybrid system [11][12][13][14][15][16] , in which gapless excitations exist at the boundary or in the vortex core, dubbed "Majorana fermions" or "Majorana zero mode (MZM)". MZM possess exotic non-Abelian statistics and thus can serve as the building block for topological quantum computation 17,18 .
Since both finite momentum pairing and topological superconductivity require magnetism and superconductivity, it is natural to ask if these two SCing phenomena can coexist and if there is any interplay between them. In particular, one may ask (1) if topological SC phases can exist for FF or LO state with finite momentum pairing; (2) how to find an experimentally feasible system for a robust realization of such state; and (3) what types of boundary modes can emerge in such system. Recent theoretical developments have revealed the possibility of two-dimensional (2D) topological FF (tFF) phase with a non-zero Chern number in the bulk and chiral Majorana modes at the boundary in cold atom systems [19][20][21][22] . More recently, it was proposed that halfvortex in certain type of topologically trivial LO phase is able to host MZMs 23 .
In this work, we propose a topologically non-trivial LO (tLO) state in a model system of bilayer SCing topological insulator (TI) thin films under an in-plane magnetic field. Based on the selfconsistent gap equation and Ginzburg-Landau free energy, we construct the phase diagram of superconductivity in this model system and demonstrate the existence of LO phase in a wide parameter range. The topological nature of this LO phase is revealed by the calculation of edge modes and the non-trivial Z 2 topological invariant in the D class. In the tLO phase, a chain of numerous MZMs, dubbed "MZM chain", is predicted along the edge perpendicular to the in-plane magnetic field. The unique experimental signature for probing such MZM chain through scanning tunneling microscope (STM) is proposed to distinguish from other topological SCing phase (e.g., tFF phase). The possible material realization of tLO phase in the SC/TI sandwich structure or in the newly discovered Fe(Te, Se) films with topological surface states is also discussed.

Results
Bilayer SCing TI film. Our model system consists of a TI film under an in-plane magnetic field, in which both top and bottom surface states (bilayer) are in proximity to the conventional swave SC pairing, as shown in Fig. 1a. The low energy physics of this system is described by the Hamiltonian Here, H 0 describes the surface states at the top and bottom surfaces under an in-plane magnetic field along the x direction and is given by 24,25 wherec ¼ĉ t;" ;ĉ t;# ;ĉ b;" ;ĉ b;# T are the electron annihilation operators, ðp x ;p y Þ are in-plane momentum operators, v is the Fermi velocity, and σ and τ are Pauli matrices, representing spin and pseudo-spin (top(t) and bottom(b) surfaces), respectively. For the Hamiltonian (2), the first term m ¼ m 0 þ m 1p 2 À Á describes the tunneling between the two surface states (called inter-layer tunneling below), the second term is the Dirac Hamiltonian for two surface states, and the third term gives the Zeeman coupling between electron spin and in-plane magnetic fields. The sketch of the Fermi surface is shown in Fig. 1b. For simplicity, we have absorbed the parameters h into v and gμ B into B x . In the calculation, we define m F ¼ m 0 þ m 1 p 2 F , which gives rise to μ 2 ¼ m 2 F þ v 2 p 2 F in the absence of magnetic field. In this work, our calculation are based on realistic parameters for the model Hamiltonian H 0 , which are described in Supplementary Note 4. To include the SC pairing, we consider electron-electron attractive interaction term 26 (2) under an in-plane magnetic field B x . c Phase diagram of superconducting TI film in the absence of magnetic field, as a function of intra-layer interaction U and interlayer interaction V. It consists of one metallic phase and three superconducting phases: the intra-layer A 1g pairingĉ t"ĉt# þĉ b"ĉb# , the inter-layer A 1u pairinĝ c t"ĉb# Àĉ b"ĉt# , and the intra-layer A 2u pairingĉ t"ĉt# Àĉ b"ĉb# wheren τ ðrÞ ¼ P σ¼";#nτ;σ ðrÞ is the electron density operator on layer τ = t, b.
The phase diagram of SCing phases can be constructed by minimizing Ginzburg-Landau (GL) free energy L obtained from the microscopic Hamiltonian (1) 27 . The quadratic GL free energy reads, where the the superconductivity susceptibility is where G e ðω; pÞ and G h ðω; ÀpÞ are the single-particle Green's functions of electrons and holes, and γ represents the superconducting pairing matrix: γ ∈ {s y , τ x s y , τ y s x , τ z s y , τ y , τ y s z }. The linearized gap function is derived from ∂L 2 =∂Δ i ¼ 0 and can be used to construct the phase diagram. More technical details of constructing the phase diagram and the realistic parameters of the microscopic model 24,[28][29][30][31][32] are discussed in Supplementary Note 3. At zero magnetic field B x = 0, the Hamiltonian (1) has D 3d group symmetry, which possesses three one-dimensional (1D) irreducible representations: A 1g , A 1u , and A 2u ; and one 2D irreducible representation: E u . Consequently, the momentumindependent SCing pairing forms can be classified 26 . According to the linearized gap equation (See Supplementary Note 2), the stable SCing phases are shown as a function of interaction parameters U and V in Fig. 1c. It is found that intra-layer A 1g pairingĉ t"ĉt# þĉ b"ĉb# is favored for a strong attractive U while inter-layer A 1u pairingĉ t"ĉb# Àĉ b"ĉt# (intra-layer A 2u pairinĝ c t"ĉt# Àĉ b"ĉb# ) can appear for a strong attractive (repulsive) V. More precisely, the phase boundary between A 1g and A 1u is, μ 2 ; and the boundary between A 1g and A 2u is U = −V. However, in the realistic materials, we expect the intra-layer interaction U dominates over the inter-layer interaction V and thus gives rise to superconductivity described by A 1g pairing in the phase diagram.
Next, we focus on the case with a strong attractive U term below and set V = 0, which is most common, and discuss the T -B x phase diagram. With magnetic field B x , two intra-layer superconducting pairings, Δ A 1g $ s y and Δ A 2u $ τ z s y , can be mixed with each other. In addition, an in-plane magnetic field is possible to induce the superconducting pairing with the finite momentum q = (0, Q). Therefore, we consider the following three types of pairing forms: For a small B x , the BCS state is energetically favored. As B x increases, the self-consistently calculated momentum Q of the pairing becomes non-zero at a critical magnetic field B c , as shown in the inset of Fig. 2b, leading a first-order phase transition to a finite momentum pairing state. For B x far above the critical value B c , we find that the momentum Q of the pairing can be approximated by 2B x /v. Two possible finite momentum pairings, FF and LO states, are considered. Our calculation suggests that the LO state is energetically favored at a large B x (phase III in Fig. 2a), while the FF state only exists in a small region (phase II in Fig. 2a) between the BCS and LO states. This small region for FF state may disappear for a smaller m F ¼ m 0 þ m 1 p 2 F where p F is the Fermi momentum (see Supplementary Note 3 for the analytical calculation in the m F = 0 limit). The critical B c depends on the coupling strength m F , as shown by different color lines in the inset of Fig. 2b. Figure 2b depicts the dependence of B c on m F , from which one can see that B c can be greatly decreased when m F is reduced toward zero.
To understand the occurrence of the LO state, we may first consider the energy spectrum of the singleparticle Hamiltonian H 0 in Eq. (2), which is given by with Q = 2B x /v due to the Zeeman term, as illustrated in Fig. 1b. Spin textures of the surface states are also depicted on the Fermi surfaces, from which one can see that zero momentum pairing can only occur for electrons with the same spin and opposite layers (inter-layer equal spin triplet pairing), while intra-layer spin-singlet pairing is only possible for a finite momentum. In the limit m F → 0, the finite momentum pairing is favored and thus two FF phases with opposite momenta for each surface state appear, similar to the case of bilayer transition metal dichalcogenides (TMDs) system 33 . A finite m F can induce Josephson coupling between two FF phases at the opposite surfaces and thus drive the whole system into the LO phase. We notice that such coupling also tends to induce the inter-layer pairing (BCS type with zero momentum) between two surfaces in order to lower the free energy L. As a result, there is a competition between Josephson coupling due to finite m F , which favors BCS pairing, and the momentum shift due to in-plane magnetic fields, which favors FF or LO state. This analysis explains the dependence of the critical magnetic field B c on m F [See Fig. 2b]. We notice the similarity between our phase diagram (Fig. 2a) and that of 2D Rashba SCs 34,35 , presumably due to the same spin textures (Fig. 1b) in these two systems.
Majorana zero mode chain. We now study on topological property of the LO phase (phase III in Fig. 2a). We consider the below Bogoliubov-de Gennes (BdG) Hamiltonian, where H ee ¼ H 0 ðk x ; Ài∂ y Þ and H hh ¼ ÀH Ã 0 ðÀk x ; Ài∂ y Þ with H 0 given by Eq. (2), and the pairing Hamiltonian MZMs are hybridized with each other and the Majorana bands become dispersive, as shown in Fig. 3b, d. We notice that these two Majorana bands cross with each other at k y = Q/2. This crossing is robust against varying magnetic field B x and implies the possible non-trivial topological nature of these two Majorana bands, which will be discussed in details below.
Topological LO state. To understand the topological nature of the system, we will develop a general theoretical framework for topological LO state. We focus on the general BdG Hamiltonian (6), here H ee needs not to be specified. Due to the periodicity of the gap function Δ(y) = Δ(y + 2π/Q) in H ee , we can expand ΔðyÞ ¼ P n Δ n e inQy with the wave-vector Q and an integer n. Only one Δ n is non-zero in the FF state while multiple non-zero Δ n exist in the LO state. Consequently, H BdG can also be expanded in the momentum space as nQ À k y Þ and H eh ð Þ nm ¼ Δ nþm (n, m are integers) on the basis e n j i and h n j i with the wave-vector nQ. Here, the momentum k y is within the reduced Brillouin zone [0,Q]. The Hamiltonian (6) possesses particle-hole symmetry C ¼ t x K where the Pauli matrix t x acts on the particle-hole space and K is the complex conjugate. Based on the Hamiltonian (6), we can extract the topological property of FF state, as shown in Supplementary Note 5, which is consistent with the recent results on topological FF state in cold atom systems [19][20][21][22] .
Next, we focus on the LO state with non-zero Δ ±1 , for which the Hamiltonian (6) can be split into two decoupled blocks. All the electron part H 0 ðk x ; nQ þ k y Þ with even (odd) n is only coupled to the hole part ÀH Ã 0 ðÀk x ; nQ À k y Þ with odd (even) n. We call these two blocks as even and odd block, denoted as H even BdG and H odd BdG , respectively. Here, the even block is written on the basis e 2n j i and h 2nþ1 , while the odd block is written on the basis e 2nÀ1 j i and h 2n j i. The global particle-hole symmetry C relates these two blocks, CH even BdG C À1 ¼ ÀH odd BdG , and thus there is in general no particle-hole symmetry within one block. However, at the momentum k y = Q/2, a new particle-hole symmetry operator C can be defined asC e 2n j i ¼ h 2nþ1 andC h 2nþ1 ¼ e 2n j i for the even block H even BdG andC e 2nÀ1 j i¼ h 2n j i andC h 2n j i ¼ e 2nÀ1 j i for the odd block H odd BdG , and we haveCH even BdG ðk x ; Q=2ÞC À1 ¼ ÀH even BdG ðÀk x ; Q=2Þ andCH odd BdG ðk x ; Q=2ÞC À1 ¼ ÀH odd BdG ðÀk x ; Q=2Þ, as shown in Supplementary Note 5. The existence of this new particle-hole symmetryC suggests that the Hamiltonian H even BdG and H odd BdG at k y ¼ Q 2 can be viewed as a 1D SC chain in the D class along the x direction. A Z 2 topological invariant 17 can be defined as, where the anti-unitary matrix A is chosen as the BdG Hamiltonian H even BdG or H odd BdG at k y ¼ Q 2 in the Majorana representation. Here, M ¼ À1 is for topologically non-trivial phase with a pair of MZM present, while M ¼ þ1 is for trivial phase. Thus, we conclude the existence of tLO state characterized by Z 2 topological invariant due to the new particle-hole symmetryC that is only valid at k y = Q/2.
Based on the above argument, we directly apply this general theory of tLO phase to the Hamiltonian in Eq. (6) that describes the bilayer SCing TI film, and evaluate the Z 2 topological invariant M defined in Eq. (7). As we discussed early, there is a robust LO phase existing in this SCing system (see the phase diagram in Fig. 2a), and the corresponding Δ LO satisfies cos ðQyÞτ z σ y , and thus we can directly evaluate Z 2 topological invariant M. Direct calculation and theoretical analysis described in Supplementary Note 6 suggest that M ¼ À1 always exists in our model once the chemical potential is large enough, thus demonstrating the robust realization of the tLO phase in the bilayer SCing TI films. The non-trivial Z 2 topological invariant suggests the existence of Majorana modes at the boundary of the system. As for the two decoupled BdG Hamiltonian (subspace) H even BdG and H odd BdG at k y ¼ Q 2 , there is a non-trivial Z 2 topological invariant in each subspace. Therefore, two degenerate MZMs protected by the new particle-hole symmetryC can exist at one edge for the momentum k y ¼ Q 2 , which is consistent with numerical results in Fig. 3. Finally, we emphasize that tLO state with Z 2 topological invariant belongs to a different topological class from tFF state with non-zero Chern number discussed in literature 19,20 and thus possesses a different type of boundary mode, the MZM chain.
In the decoupling limit with m F → 0, the LO phase in our system can be viewed as two decoupled FF states, and each FF state can be adiabatically connected to the Fu-Kane model 36 by tuning magnetic field to zero. It is well-understood that although Fu-Kane model can host Majorana zero mode inside the vortex core, it does not possess chiral Majorana edge state due to the time reversal symmetry. Since there are two FF states in our system, we expect two Majorana zero modes inside the vortex core and thus the statistics is no longer non-Abelian when exchanging two vortices. Therefore, our system possesses neither chiral (or helical) Majorana edge states nor non-Abelian statistics through exchanging vortex core in the decoupling limit. This also makes the tLO phase found here different from the tFF phase with chiral Majorana edge modes discussed in literature.

Discussion
In this work, we develop a general theory of tLO phase with Z 2 classification and propose its realization in a theoretical model of bilayer SCing TI films. The tLO phase and the corresponding MZM chain also open a new route in the study of MZMs for topological quantum computation. The 1D MZM chain also provides a natural platform to study the interacting Majorana chains 37,38 . Below we will discuss the possible material realization of our model and the possible unique experimental signatures to distinguish the tLO phase from other topological SCing phase.
Our proposed model can be realized in SC/TI/SC heterostructure 30,31,39 , e.g., NbSe 2 /Bi 2 Te 3 /NbSe 2 heterostructure. With Δ 0 = 1 meV, μ = 100 meV, ℏv = 2.67 eV ⋅ Å, the g-factor g = 23, and m F = 20.1 meV, we can estimate the critical field at tricritical point to be about 0.15 Tesla according to gμ B B c /k B T c,0 ≈ 0.35 from Fig. 2b. The distance between the two MZMs is estimated as l y = πℏv/4gμ B B c ≈ 1.2 μm, which is four times larger than the localization length of MZMs ξ $ hv=Δ 0 ¼ 0:27 μm < l y . Thus, MZMs in the chain are well separated and can be resolved in a scanning tunneling microscope experiment 39 . The above estimate is based on Zeeman effect, but we emphasize that the orbital effect of in-plane magnetic fields can also plays a similar role as the Zeeman effect due to the Dirac fermion nature. Compared to the Zeeman effect, we find the orbital effect of an in-plane magnetic field can also induce the FFLO phase (the vector potential may be chosen asÃ ¼ ð0; ÀB x z; 0Þ and set the middle of layers as z = 0 resulting in the opposite momentum shift for the Fermi surface of top and bottom surface states), and it is about π hvd=2Φ 0 $ 0:91 meV=Tesla by assuming the space distance between two surfaces d = 3 nm, which is comparable to the Zeeman term with gμ B = 1.33 meV/Tesla. One can also consider SC/magnetic TI/SC heterostructure, in which the exchange coupling from magnetic doping takes a similar form as Zeeman effect, but is two orders of magnitude larger than the Zeeman effect 25 . It should be mentioned that the 1D Z 2 topological invariants at k y = Q/2 replies on the translation symmetry and the tLO phase is not robust against strong disorder scattering, thus requiring a relatively clean sample in experiments. We would like to mention a recent experiment on the evidence of finite momentum pairing induced by magnetic fields in HgTe quantum wells 40 and 3D TI Bi 2 Se 3 41 in proximity to superconductivity, and our results suggest the possibility of non-trivial topology of this finite momentum pairing. Based on the above estimate, we conclude that our proposal is feasible under the current experimental conditions.
Our proposal is also applicable to SCing TIs in which topological surface states and bulk superconductivity can coexist. Such materials include Cu doped Bi 2 Se 3 42-45 , several SCing half-Heusler compounds (e.g., YPtBi, RPdBi) 46,47 and Fe(Te, Se) films 48,49 . We emphasize that the s-wave spin-singlet superconductivity has been demonstrated for the surface states of TI Fe (Te, Se) film 49,50 , which may provide a natural high-Tc SC platform to realize our model. In addition, odd parity pairing has been proposed in several of the above compounds 26,51 , and the influence of in-plane magnetic field will be interesting for a future study. Our work provides an example of the interesting interplay between finite momentum pairing and topological physics and similar idea also exists for exciton condensate 52 .
Finally, we hope to discuss how to unambiguously distinguish the tLO phase from other topological SCing phases, particularly the tFF phase, in experiments. The key difference lies in the fact that the chiral Majorana edge mode in the tFF phase can exist at any edge of the system while the MZM chain in the tLO phase can only exist at the edge perpendicular to the in-plane magnetic field. As shown in Fig. 4, we compare the local density of states (LDOS) at the edge when the in-plane magnetic field is perpendicular (blue lines) or parallel (red line) to the edge. The calculation details of LDOS can be found in Supplementary Note 7. A zero-bias peak is found for the perpendicular magnetic field direction, demonstrating the existence of MZM chain, while a full SCing gap is found for the parallel magnetic field. This calculation suggests that a simple scanning tunneling microscopy measurement of the edge LDOS with parallel and perpendicular magnetic fields can give rise to unambiguous experimental signature to identify the tLO phase and the MZM chain.

Methods
In the Supplementary Note 1 and Supplementary Note 2, we review the Ginzburg-Landau theory and use the linearized gap equation to classify the possible SCing pairings of SCing TI thin films in the absence of magnetic field. Then, we show the details for the self-consistent calculation for the T − B x phase diagram by minimizing the Ginzburg-Landau free energy in Supplementary Note 3, show how to construct the tight-binding model in Supplementary Note 4, analyze the calculation of topological invariant in Supplementary Note 6, and discuss how to calculate LDOS in Supplementary Note 7.

Data availability
The data and codes that support the findings of this study are available from the corresponding author upon reasonable request.