Pairing symmetry and topological surface state in iron-chalcogenide superconductors

The symmetries of superconducting gap functions remain an important question of iron-based superconductivity. Motivated by the recent angle-resolved photoemission spectroscopic measurements on iron-chalcogenide superconductors, we investigate the influence of pairing symmetries on the topological surface state. If the surface Dirac cone becomes gapped in the superconducting phase, it implies magnetization induced from time-reversal symmetry breaking pairing via spin-orbit coupling. Based on the crystalline symmetry constraints on the Ginzburg-Landau free energy, the gap function symmetries are among the possibilities of $A_{1g(u)}\pm iA_{2g(u)}$, $B_{1g(u)}\pm iB_{2g(u)}$, or, $E_{g(u)}\pm i E_{g(u)}$. This time-reversal symmetry breaking effect can exist in the normal state very close to $T_c$ with the relative phase between two gap functions locked at $\pm \frac{\pi}{2}$. The coupling between magnetization and superconducting gap functions is calculated based on a three-orbital model for the band structure of iron-chalcogenides. This study provides the connection between the gap function symmetries and topological properties of the surface state.

The symmetries of superconducting gap functions remain an important question of iron-based superconductivity. Motivated by the recent angle-resolved photoemission spectroscopic measurements on iron-chalcogenide superconductors, we investigate the influence of pairing symmetries on the topological surface state. If the surface Dirac cone becomes gapped in the superconducting phase, it implies magnetization induced from time-reversal symmetry breaking pairing via spinorbit coupling. Based on the crystalline symmetry constraints on the Ginzburg-Landau free energy, the gap function symmetries are among the possibilities of A 1g(u) ± iA 2g(u) , B 1g(u) ± iB 2g(u) , or, E g(u) ± iE g(u) . This time-reversal symmetry breaking effect can exist in the normal state very close to Tc with the relative phase between two gap functions locked at ± π 2 . The coupling between magnetization and superconducting gap functions is calculated based on a three-orbital model for the band structure of iron-chalcogenides. This study provides the connection between the gap function symmetries and topological properties of the surface state.
However, the non-trivial topology of the surface superconductivity in FeSe 0.5 Te 0.5 is mostly a property inherited from the normal state band structure in a similar way to the Fu-Kane proposal of two-dimensional topological superconductivity via the proximity effect [58]. It does not directly reveal the symmetry properties of the superconducting gap functions. It would be highly desirable if the topological surface states could be used for phase-sensitive detections to unconventional pairing symmetries [3,4]. In contrast, in a very recent laser-based angle-resolved photoemission spectroscopy (APERS) experiment on FeSe 0.3 Te 0.7 [59], the surface Dirac cone is observed to develop a gap as the system enters the superconducting state. The surface Dirac point is well-below the Fermi energy. Thus the splitting cannot be the superconducting gap, but implies TR symmetry breaking in the spin channel, directly correlated with the superconducting transition. In earlier literature, both TR symmetry breaking pairing of the types s + id [34,[60][61][62][63] and s + is [64][65][66][67] have been proposed. However, in both cases magnetization only appears around impurities.
In the present work, we investigate how the topological surface states are affected by TR breaking gap functions, which in turn constrains the possible pairing symmetries. By employing the Ginzburg-Landau formalism, we explore possible TR breaking gap functions which can induce magnetization via spin-orbit coupling to split the degeneracy at the surface Dirac point. Based on crystalline symmetry analysis, the superconductivity gap symmetries include the possibilities of A 1g(u) ± iA 2g(u) , B 1g(u) ± iB 2g(u) , and E g(u) ± iE g(u) . In the normal state sufficiently close to T c , the relative phase between two gap functions can still be locked at ± π 2 even though neither of them is long-range ordered. Calculations based on a three-orbital model are performed to derive the coupling between magnetization and superconducting gap functions. In doing so, our study bridges the topological properties of the surface state and the pairing symmetries of the superconducting gap functions.
We begin with a discussion of the splitting of the surface Dirac cone in the superconducting state. In the FeSe 1−x Te x materials, the degeneracy of the surface Dirac point is protected by the band topology if the normal state maintains TR symmetry. Consequently a surface Dirac cone appears at the Γ-point, as shown in Fig. 1 (a), described by an effective k · p Hamiltonian where σ's are Pauli matrices defined for the Kramers doublet at the Γ-point, and µ is the chemical potential. A superconducting gap ∆ by itself, i.e., H Γ = −µ + ∆w − + ∆ * w + where w ± = w x ± iw y are Pauli matrices in the Nambu space of the particleparticle channel, does not lift the degeneracy. To split the degeneracy, a mass term breaking TR symmetry is necessary, i.e., ∆H Γ = −m z σ z , and the associated Bogoliubov spectra become ±m z + |∆| 2 + µ 2 around the Γ-point. The sketch of the Bogoliubov dispersion is illustrated in Fig. 1(b) with both particle and hole branches. The mass term corresponds to a splitting between two eigenstates of σ z . Since no magnetic field is applied, m z should arise from the Weiss field of a ferromagnetic ordering along the z-axis induced by superconductivity. The surface state spectra in the superconducting state (T < Tc). Two gaps appear: the superconducting gap at the chemical potential µ, and the splitting of the Dirac cone is due to the magnetic ordering breaking TR symmetry.
Now we examine how the ferromagnetic order m z can be induced in the superconducting state. Apparently, this requires the spontaneous breaking of TR symmetry. Indeed, within the Ginzburg-Landau (GL) formalism [2,34,37,68,69], it has been shown that the mixing between two gap functions ∆ 1,2 , which are TR invariant by themselves and possess different pairing symmetries, leads to the spontaneous TR symmetry breaking. The corresponding physical consequences were studied in the case of iron-based superconductors. ∆ 1,2 cannot form a symmetry invariant at the quadratic level, but they do at the quartic level via The β ′ -term locks the relative phase between two gap functions, which equals β ′ |∆ 1 | 2 |∆ 2 | 2 cos 2∆ϕ, where ∆ϕ = ϕ 1 − ϕ 2 , and ϕ 1,2 are the phases of two gap functions. When β ′ > 0, ∆ϕ is pinned at ± π 2 , giving rise to complex gap functions ∆ 1 ± i∆ 2 , which break TR symmetry spontaneously. This formalism also applies to the case that ∆ 1,2 form a two-dimensional (2D) irreducible representation. The complex gap functions ∆ 1 ± i∆ 2 distribute more evenly over the Fermi surface than the real ones ∆ 1 ± ∆ 2 , hence, they are energetically more preferable at the mean-field level [37]. The corresponding Cooper pairs carry non-zero orbital moments, which could generate magnetic fields at boundaries as shown earlier [68,69]. However, these magnetic fields are typically of the order of 1 Gauss, for which the Zeeman energy is negligible. Instead, here we consider the spin magnetization m z coupling to ∆ 1,2 through a cubic term as This satisfies both the U (1) and TR symmetry. α > 0 is assumed in Eq. 2, and hence there is no spontaneous magnetic ordering by itself, rather, the magnetization is induced, m z = γ α |∆ 1 ∆ 2 | sin ∆ϕ, i.e, by coupling to the TR breaking superconducting orders. The sign of m z is determined by the relative phase ∆ϕ between ∆ 1,2 .
The free energy density F M of Eq. 2 is further required to satisfy all crystalline symmetries. At elevated temperatures the pristine FeSe and FeTe crystals are layered quasi-2D systems, whose Bravais lattices are primitive tetragonal. They exhibit a tri-layer structure with each unit cell consisting of two Fe cations and two Se(Te) anions: A square lattice of Fe cations in the middle layer sandwiched between two layers of Se(Te) anions in a √ 2 × √ 2 structure. The Se(Te) lattices above and below the iron planes are off-set by one Fe-Fe bond length, and their projections are at iron plaquette centers. The crystalline space group is the non-symmorphic one P 4/nmm [16,42,49,50], which is reviewed in the Supplemental Material (S.M.) I. It can be decomposed into 16 cosets: 8 of them are denoted as g i T where T is the translation group of the primitive tetragonal lattice and g i (i = 1 ∼ 8) span the point group C 4v centering at Se(Te) anions, and the other 8 cosets are Ig i T by further applying the inversion I with respect to the Fe-Fe bond center. In the actual experimental systems of FeSe x Te 1−x , the distribution of Se and Te breaks the P 4/nmm symmetry, nevertheless, this effect is weak after averaging over random configurations and will be neglected below. As shown in Fig. 2, the rotations with respect to Se/Te and Fe are 4-fold and 2-fold denoted as C 4 (z) and C 2 (z), respectively. The point group symmetry centering around the Fe cations is D 2d . The vertical reflection planes are along x, y denoted as σ x , σ y , and are along the diagonal lines x ′ and y ′ denoted as σ x ′ and σ y ′ , respectively.
We consider the order parameter properties under the crystalline symmetries. In multi-orbital systems, the su- perconducting gap function is expressed as where repeated indices mean summation; a, b refer to the orbital band components, and i, j are the sublattice indices of two Fe-cations in one unit cell; τ is a 2 × 2 matrix representing the sublattice channel; f l (k) is the angular form factor of momentum k, and M l is the pairing matrix in the orbital channel; l is the index for multiple combinations between f l (k) and M l . The pairing matrix is defined as ∆ bj,ai (k) = k iσ y,αβ |c † k,αbj c † −k,βai | , where iσ y projects out the singlet pairing with Greek indices representing spin components, and || represents averaging over the thermal equilibrium state.
As required by Fermi statistics, for singlet pairing, the product τ ij f (k)M ab in Eq. 3 needs to be even under the combined operations of k → −k and the transposes of M and τ . f (k) can be an even function taking the forms of 1, cos k x ± cos k y , cos k x cos k y , sin k x sin k y , cos(k x ± k y ). or, an odd function among sin k x , sin k y and sin(k x ± k y ). We choose the three t 2g -orbital bases, Hermitian matrix which is expanded in terms of the Gellmann matrices λ i (i = 1 ∼ 8) under the basis in the sequence of (d x ′ z , d y ′ z , d xy ), and the 3 × 3 identity matrix λ 0 , whose expressions are presented in S. M. II.
The representation of a gap function under the crystalline symmetry group is determined by the symmetry properties of f (k), M , and τ as analyzed and presented in the S. M. II. Their possible symmetries are denoted as A 1g(u) , A 2g(u) , B 1g(u) , B 2g(u) , and E g(u) , respectively, where g and u represent even and odd parities, respectively. The A, B, and E symbols represent the discrete angular momentum, loosely speaking, they are analogues to the s, d, and p-wave symmetries, respectively. A 1 and A 2 exhibit even and odd parities under vertical reflection planes, for example, the ferromagnetic order m z carries the A 2g symmetry. B 1 and B 2 are analogous to the d xy and d x 2 −y 2 symmetries, respectively, exhibiting opposite parities under the σ x(y) and σ x ′ (y ′ ) operations. Symmetries of singlet channel gap functions are classified accordingly: The next-nearest neighbor (NNN) pairings are summarized in Tab. III and the neighbor(NN) pairings in Tab. IV in S. M. II. Different combinations of f (k), M , and τ often lead to the equivalent symmetries, and in general, the existence of one can induce others in the same symmetry class. Many orbital-dependent gap functions have been proposed in the literature [70][71][72][73], and their importance have been analyzed in recent experiments [74,75]. Due to the multi-orbital nature, the gap functions rigorously speaking cannot be intuitively represented by the partial-wave channels alone, i.e., the symmetry of the angular form factor f (k). For example, consider the following two gap functions with even parity, which carry the d xy and d x 2 −y 2 -like symmetries, or, more precisely, B 1g and B 2g symmetries, respectively, although their angular form factors are s-wave like. They involve the intra-and inter-orbital pairings between the d x ′ z and d y ′ z -orbitals as shown in Fig. 3 (a) and (b). Since d xz → d yz and d yz → −d xz under the 90 • rotation, λ 1,3 transform analogously in the d-wave way. By examining their reflection symmetries, they belong to the B 1g and B 2g symmetries. By similar analysis, the following gap functions, exhibit the A 1g and A 2g symmetries, respectively, or, loosely speaking, the s-wave symmetry, in spite of their d x 2 −y 2 angular form factor: Their orbital configurations are shown in Fig. 3 (c) and (d). Furthermore, there can exist p-wave like pairing symmetry, or, the E g -symmetry, in the singlet pairing channel, The former (latter) describes the pairing between the d x ′ z -orbital (d y ′ z ) with the d xy one.
The crystalline symmetries impose stringent constraints to the superconducting gap functions. According to Eq. 2, the direct product of the irreducible representations of ∆ 1 and ∆ 2 should contain that of m z , i.e., A 2g . This yields the following possibilities of pairing symmetries: . Examples of the above pairing symmetries with even parity are provided in Eqs. 4, 5, and 6.
An important issue is that spin-orbit coupling is necessary to break the SU(2) symmetry such that the ferromagnetic order m z and superconducting orders ∆ 1,2 can couple, since the former and latter lie in the spin triplet and singlet channels, respectively.
We employ a widely used three-band model for the topological band structure of FeTe 1−x Se x around the Γpoint, which consists of the t 2g -orbitals d x ′ z , d y ′ z and d xy [44,51]. Neglecting the small dispersion along the z-axis, the 3-band tight-binding Hamiltonian is expressed as [76], T and the matrix kernel H 0 (k) is given by, where H N N N and H N N represent the NNN and NN hoppings, respectively, with detailed forms presented in S. M. III. H soc = λ soc τ 0 L · σ is the atomic spin-orbit coupling with λ soc the coupling strength and L representing the onsite orbital angular momentum projected to the t 2g -basis. Explicitly, L = (λ 5 − λ 7 )/ √ 2, −(λ 5 + λ 7 )/ √ 2, −λ 2 . Based on the band structure Eq. 7, the coupling coefficient γ in Eq. 2 can be evaluated as where β = 1/(k B T ) is the inverse of temperature; f li i and M li i with i = 1, 2 are the angular form factors and orbital pairing matrix kernels of ∆ 1,2 , respectively; G e (k, iω n ) = (iω n − H 0 (k)) −1 is the Matsubara Green's function with ω n = (2n + 1)π/β, and G h (k, iω n ) = G * e (−k, −iω n ). As shown in S. M. III, the hole-like Fermi pockets around the Γ-point are mainly from the bonding states between two Fe-sublattices, i.e., they are approximately eigenstates of τ 1 with the eigenvalue of 1. Hence, only tracing over the spin and orbital channels are needed, and only these gap functions characterized by τ 0,1 are considered. Gap functions with τ 2,3 are pairing between bonding and anti-bonding states between two-sublattices, which will be neglected below.
Next we present the examples of gap functions ∆ 1,2 leading to the spontaneous magnetization m z . We begin with the cases of B 1g(u) ± iB 2g(u) . For parity even, i.e., B 1g ± iB 2g , we take ∆ 1,2 in the form of Eq. 4. As shown in S. M. IV, after further reducing the band Hamiltonian Eq. 7 to a two-band model only based on d x ′ (y ′ )z , Eq. 8 yields an analytic expression as with N 0 the density of states at the Fermi surface.
Since m z is induced by the TR breaking pairings via spin-orbit coupling, the coupling coefficient γ is proportional to the spin-orbit coupling strength. A calculation based on the 3-band Hamiltonian is also performed numerically, which yields consistent results (see S. M. V). The parity odd case, i.e., B 1u ± iB 2u , is also numerically checked to yield a nonzero γ, for example, with ∆ 1 : τ 0 (sin(k x + k y )λ 7 + sin(k x − k y )λ 5 ) and ∆ 2 : τ 0 (sin(k x + k y )λ 5 − sin(k x − k y )λ 7 ). The above two cases also break the mirror symmetries of σ x(y) and σ x ′ (y ′ ) spontaneously. They are topologically non-trivial belonging to the C-class supporting the chiral Majorana edge modes [77][78][79][80]. We have also studied both cases of A 1g(u) ± iA 2g(u) , which also yield non-zero γ's as shown in S. M. V. The nodal pairing gap functions presented in Eq. 5 are used for the even parity case, and the nodeless pairing with ∆ 1 : τ 0 (sin(k x + k y )λ 5 + sin(k x − k y )λ 7 ) and ∆ 2 : τ 0 (sin(k x + k y )λ 7 − sin(k x − k y )λ 5 ) are used for the odd parity case. For the case of E g(u) ± iE g(u) , we take the gap functions presented in Eq. 6 as an example.
Since E g(u) ⊗E g(u) = A 1g ⊕A 2g ⊕B 1g ⊕B 2g , it also yields a nonzero γ as calculated in S. M. IV, and consequently induces magnetization.
In strongly correlated superconductors, there exist strong superconducting phase fluctuations in the normal state close to T c [81]. In this case, the phases ϕ 1 and ϕ 2 of gap functions are disordered such that ∆ 1,2 = 0, but |∆ 1,2 | remains finite. The β ′ term in Eq. 1 can still pin the relative phase ∆ϕ = ± π 2 . This transition breaks TR symmetry and its critical temperature T ′ > T c . We still expect a weaker but still finite m z in the temperature window between T c and T ′ .
The bulk spin magnetization induced by TR breaking pairing qualitatively explains the gap opening of the surface Dirac cone observed in the FeSe 0.3 Te 0.7 superconductor [59]. We next briefly discuss the issue of the possible magnetic field generated by the bulk magnetization. Its upper bound is estimated to be 15 Gauss (See S. M. VI), which is still smaller than the lower critical magnetic field of the FeSe 1−x Te x superconductors (H c1 ∼ 30 Gauss) [82,83]. The actual field due to spin magnetization should be much smaller than this bound, hence, it can be offset by the orbital magnetization such that the total magnetic field remains zero in the Meissner state [84,85]. This is consistent with the observation that neutron spectroscopy does not detect a bulk magnetic field [86]. A very weak but finite internal magnetic field around 0.15 Gauss is detected in the FeSe superconductor by the muon spin rotation measurement (µSR) [87,88], which may arise from the imperfect screening due to impurities and domains.
Discussion and Conclusion-. In this article, we have studied how the TR symmetry breaking superconducting states can gap out the topological surface modes in the iron-chalcogenide superconductors. Spin-orbit coupling is necessary to break the SU(2) symmetry in the spin channel, such that it bridges the magnetic ordering and the TR breaking pairing states. Three classes of gap function symmetries can lead to such an effect based on group theory analyses: A 1g(u) + iA 2g(u) , B 1g(u) + iB 2g(u) , E g(u) + iE g(u) . In strongly correlated superconductors, the superconducting phase fluctuations can also lock their relative phase at ± π 2 breaking TR symmetry in the normal state. This work builds connections between novel pairing symmetries in iron-based superconductors and their topological band structures. In particular, it is helpful for understanding the superconductivity induced gap opening for the topological surface state recently discovered in the FeTe 0.7 Se 0.3 [59].

P4/nmm non-symmorphic group
Here we review the P 4/nmm group for the FeSe crystal.
If the Fe site is chosen as the origin, there exists the D 2d point group symmetry centering at the Fe site, consisting of operations of E, S 4 , C ′ 2 (z), S 3 4 , C 2 (x), C 2 (y), σ x ′ , σ y ′ , where C ′ 2 (z) is denoted to distinguish the C 2 (z) centering at Se/Te site below. Each unit cell contains two Fe and Se atoms exhibiting a square lattice with the lattice constant √ 2a where a measures the nearest neighbor Fe-Fe bond in the x-y plane. It is also convenient to choose the Se/Te site as the origin, then the corresponding point group is C 4v = E, C 4 (z), C 2 (z), C 3 4 (z), σ x , σ y , σ x ′ , σ y ′ , which is used in this work to specify the irreducible representations of gap functions. The configurations of the reflection planes of σ x , σ y , σ x ′ and σ y ′ are given in Fig. 2 in the main text.
The P 4/nmm space group of the FeSe crystal can be decomposed into the cosets denoted as gT , where g is a symmetry operation; T is the lattice translation group consisting of translations of i=1,2,3 l i a i with integers l 1,2,3 , and the lattice vectors are defined as a 1 = (a, a, 0), a 2 = (a, −a, 0) and a 3 = (0, 0, c) with c the length of Fe-Fe bond along the z-direction. There are 8 cosets generated by g ∈ C 4v point group: where C 1,3 4 (z) and C 2 (z) are rotations around the z-axis passing the Se atom, and σ x(y) and σ x ′ (y ′ ) are reflections with respect to the vertical planes passing the Se atoms as shown in Fig. 2 in the main text. The rest cosets are constructed by applying the inversion operation I with respect to the Fe-Fe bond centers to the previous 8 ones, which are where S 4 and S 3 4 are the rotary reflection centering around the Fe cation; C 2 (x) and C 2 (y) are the 2-fold rotations around the nearest Fe-Fe bonds. There are three class of non-symmorphic operations: g(σ h , τ ) is the glide reflection with σ h the reflection with respect to the xyplane and τ = a(1, 0, 0); g(C 2 (x ′′ ), τ x ) and g(C 2 (y ′′ ), τ y ) are screw rotations with their axes x ′′ and y ′′ along the 45 • and 135 • degrees passing the Fe-Fe bond centers; τ x = 1 2 (−a, a, 0) and τ y = 1 2 (a, a, 0). The character table of C 4v is given in Table. I, where there are four one-dimensional irreducible representations A 1 , A 2 , B 1 , B 2 and one two-dimensional representation E. Taking the inversion symmetry of P 4/nmm into account, the total representations are: A 1g/1u , A 2g/2u , B 1g/2u , B 2g/2u and E g/u , where g and u mean even and odd under the inversion symmetry, respectively. According to the crystal symmetry, the d x ′ z , d y ′ z , d xyorbitals on the A and B sublattices of Fe are constructed shown in Fig. 4, where the d x ′ z and d y ′ z are along the diagonal Fe-Fe lines. For the three dimensional orbital space, we define the 3 × 3 Gell'mann matrices according to the basis of (d x ′ z , d y ′ z , d xy ) T as In addition, λ 0 is the 3 × 3 identity matrix. We also denote the Pauli matrices in the channel of the Fe A/B sublattice as τ 1,2,3 and τ 0 the 2 × 2 identity matrix. Similarly, the Pauli matrices in the spin channel is denoted as σ x,y,z , and σ 0 is also the 2 × 2 identity matrix.
The classification for matrices in the sublattice channel (τ ), the orbital channel (λ), and the spin channel (σ) according to the C 4v and inversion symmetries are summarized in the following Table. II.

Irreps. A/B Sub-lattice
Orbital τ 's in the channel of the iron A/B sublattice, λ's in the orbital channel, and σ's in the spin channel. For the Eg representation, the pair of the λ-matrices inside the braces form a two-dimensional basis.
The pairing gap functions in the spin singlet channels and their symmetry properties under the C 4v group and the inversion operations are systematically analyzed. Symmetries for pairings across the NNN Fe-Fe bonds are summarized in Tab. III, and those for pairings along the NN Fe-Fe bonds are summarized in Tab. IV.

Three-orbital model
The corresponding tight-binding model can be constructed below,

and the matrix kernel is defined as
The H N N N term describes the next-nearest neighbor (NNN) hopping Hamiltonian, ) cos k x cos k y 2λ 0 + √ 3λ 8 /3 + (4t nnn xy cos k x cos k y + ∆ xy ) λ 0 − √ 3λ 8 /3 where ∆ xy is the the on-site energy difference between the d xy -orbital and the two degenerated d x ′ z/y ′ z ; the hopping integrals t nnn σ,π,η,xy here describe the bonding along the NNN Fe-Fe bond (see Fig. 5(a-d)). t nnn σ and t nnn π are the σ and π-bonding strengths of the d x ′ (y ′ ),z orbitals; t nnn xy is the bonding strength between two d xy -orbitals; t nnn η is the bonding between d x ′ (y ′ )z and d xy -orbitals. The H N N term represents the nearest neighbor (NN) hopping Hamiltonian, where t nn 1,2,3,xy describe the bonding along the NN Fe-Fe bond (see Fig. 5(e-h)). t nn 1 and t nn 2 describe the intraand inter-orbital bonding between d x ′ z and d y ′ z -orbitals; t nn xy is the bonding between d xy -orbitals; t nn 3 is the bonding between d xy and d x ′ (y ′ )z orbitals.
The last term H soc is the spin-orbit coupling Hamiltonian, where L = (λ 5 − λ 7 )/ √ 2, −(λ 5 + λ 7 )/ √ 2, −λ 2 . Since τ 1 is conserved in Eq.(17), we can label the eigenstates of H 0 by τ 1 's eigenvalues of ±1. There are mainly three hole-like pockets around the Γ-point, and all of them carry the eigenvalue of 1 of τ 1 . Projecting H 0 into this sector, we arrive at a 6 × 6 Hamiltonian matrix as The corresponding band structure is calculated along the high symmetry lines Γ-X-M -Γ as shown in Fig. 6 with the parameters given in the figure caption. Each band is doubly degenerate due to the TR symmetry and inversion symmetry.
The reduced two-orbital model Due to the relatively large value of ∆ xy , we can project out the d xy -orbital and arrive a reduced two-orbital model with only d x ′ z and d y ′ z -orbitals, whose matrix kernel readH with H N N N = 2t 1 cos k x cos k yλ0 + 2t 2 sin k x sin k yλ3 , H N N = 2t 5 (cos k x + cos k y )λ 0 + 2t 7 (cos k x − cos k y )λ 1 , wheret nnn σ,π are hopping elements for NNN Fe-Fe bonding; t nn 1,2 are those for the NN Fe-Fe bonding; theλ 1,2,3 are Pauli matrices defined in the orbital channel (d x ′ z and d y ′ z ), andλ 0 is the 2 × 2 identity matrix. Then there will be two hole-like pockets around the Γ-point and all the states on the Fermi surface are eigenstates of τ 1 with the eigenvalues +1. Within this sector, theH 0 can be further simplified as whereH soc = λ socλ2 σ z . Then the Green's function can be analytically solved as, where β = 1/T (T is temperature) and the fermion Matsubara frequency ω n = (2n + 1)π/β. The projection op-erators are defined as with where ǫ 0 = 2t 1 cos k x cos k y + 2t 5 (cos k x + cos k y ). This simplification helps to calculate the Feymann diagram analytically discussed in the main text.

Numerical calculations of γ
Based on the Hamiltonian in Eq. (21), we numerically calculate the Green's function, and then evaluate the Feymann diagram discussed in the main text. A very low temperature is used with T = 0.004t 0 , which is two orders smaller than the Fermi energy relative to the band top. The results are shown in Fig. 7. In the weak spinorbit coupling λ soc region, γ is linearly proportional to λ soc . The five examples for ∆ 1,2 are studied with the following angular form factors and pairing matrices.
2.) Example of A 1u ± iA 2u case: 3.) Example of B 1g ± iB 2g case: 4.) Example of B 1u ± iB 2u case: 5.) Example of E g ± iE g case: However, we find that for the E u ± iE u case, the trace operation in Eq. 4 in the main text already yields zero based on the above three-band model, yielding γ = 0, although a non-zero γ is allowed by symmetry.   The effect of the surface Neel order We discuss whether the antiferromagnetic Neel ordering at the Fe-sites along the z-direction can gap out the surface Dirac cone, whose symmetry belongs to the B 2u representation. It cannot be induced by mixing pairing order parameters with the same parity in the bulk, for example, for the s + id-pairing, or, more precisely, A 1g + iB 2g -pairing [34]. Nevertheless, the inversion symmetry is broken at the surface, it could appear at the surface for the above pairing. Even though, the surface Dirac cone would still remain gapless due to the protection from the mirror symmetries σ x and σ y along two orthogonal directions which remains in the presence of the Neel ordering. Hence, the ferromagnetic ordering m z is necessary to gap out the Dirac cone.
Next we present a detailed calculation. We define the orbital band bases of {|i, σ , } with i = 1 ∼ 4 and σ =↑, ↓ as follows, where φ A and φ B are the three t 2g d-orbitals of the Fe atoms at A and B-sublattices, respectively. The Neel ordering along the z-direction does not have the diagonal matrix elements in the above bases due to the vanishing of the overall magnetization. It does have the off-diagonal matrix elements between |1 σ and |4 σ , and its Hamiltonian matrix is given by H We find that the surface Dirac cone is still gapless in the presence of N z by carrying out a numerical calculation for the spectra for a 3D lattice system with the open boundary condition, shown in Fig. 8. In fact, this twofold degeneracy is protected by two reflection symmetries with respect to two perpendicular mirror planes, i.e., σ x and σ y . When applying to half-integer spin fermions, actually σ x and σ y anti-commute with each other, giving rise to the protected double degeneracy.

Estimation of the internal magnetic field
In the time-reversal breaking superconducting states, there is a spontaneous magnetization in the spin channel.
Assuming each Fe site has a magnetization of bµ B , where µ B = e /2m e c is the Bohr magneton, and b is the relative magnetization, then the magnetization strength M s is given by, where a 0 = 2 /m e e 2 is the Bohr radius, (r s a 0 ) 3 is average volume containing one iron cation. Therefore, the induction magnetic field B s is, where Φ 0 = hc/2e is the magnetic flux quantum and α is the fine-structure constant.
In the FeSe 1−x Te x superconductors, the lattice constant in the a-b plane is about a = 0.379 nm, and c = 0.596 nm in the c-direction, hence r 3 s is estimated as 680. The upper bound of the relative magnetization b can be estimated as the total hole density, which is 6πk 2 F /(2π) 2 ≈ 3π/200 ≈ 4.7% with k F ≈ π/10 as the average Fermi momentum of the three hole pockets (see Fig. 6). Therefore, B s ≪ B s,upper = 15 Gauss.
This estimation of the upper bound of B s is at the same order as the saturation level of the magnetization curve of FeSe system, which is still smaller that the lower critical field (30 Gauss) of the superconducting FeSe.