Pairing symmetry and topological surface state in iron-chalcogenide superconductors

TL;DR

This paper explores how different pairing symmetries in iron-chalcogenide superconductors influence topological surface states, especially regarding time-reversal symmetry breaking and the resulting gapped Dirac cones.

Contribution

It identifies possible pairing symmetries that break time-reversal symmetry and links these to topological surface state properties in iron-chalcogenide superconductors.

Findings

Time-reversal symmetry breaking pairing can gap the surface Dirac cone.

Specific pairing symmetries are constrained by crystalline symmetry.

Coupling between magnetization and gap functions is modeled for iron-chalcogenides.

Abstract

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.

Tables

Table 1. Table 1: The character table for the point group C 4 ​ v subscript 𝐶 4 𝑣 C_{4v} .

$C_{4v}$	$E$	$2C_4\left(z\right)$	$C_2\left(z\right)$	$2{\sigma }_{x^{\prime }}$	$2{\sigma }_x$
$A_1$	$1$	$1$	$1$	$1$	$1$
$A_2$	$1$	$1$	$1$	$-1$	$-1$
$B_1$	$1$	$-1$	$1$	$1$	$-1$
$B_2$	$1$	$-1$	$1$	$-1$	$1$
$E$	$2$	$0$	$-2$	$0$	$0$

Table 2. Table 2: Representation for matrices in different channels: τ 𝜏 \tau ’s in the channel of the iron A/B sublattice, λ 𝜆 \lambda ’s in the orbital channel, and σ 𝜎 \sigma ’s in the spin channel. For the E g subscript 𝐸 𝑔 E_{g} representation, the pair of the λ 𝜆 \lambda -matrices inside the braces form a two-dimensional basis.

Irreps.	A/B Sub-lattice	Orbital	Spin
$A_{1g}$	${\tau }_0$, ${\tau }_1$	${\lambda }_0$, ${\lambda }_8$	${\sigma }_0$
$A_{2g}$	/	${\lambda }_2$	${\sigma }_z$
$B_{1g}$	/	${\lambda }_1$	/
$B_{2g}$	/	${\lambda }_3$	/
$B_{2u}$	${\tau }_2$, ${\tau }_3$	/	/
$E_g$	/	$\{{\lambda }_4,{\lambda }_6\}$, $\{{\lambda }_5,{\lambda }_7\}$	$\{{\sigma }_x,{\sigma }_y\}$

Table 3. Table 3: Classifications of the singlet superconducting order parameters for the NNN pairings according to the symmetries of C 4 ​ v subscript 𝐶 4 𝑣 C_{4v} and inversion. k ± = k x ± k y subscript 𝑘 plus-or-minus plus-or-minus subscript 𝑘 𝑥 subscript 𝑘 𝑦 k_{\pm}=k_{x}\pm k_{y} ; u 𝑢 u and g 𝑔 g represent the odd and even parities for the inversion with respect to the Fe-Fe bond centers, respectively. Each line for the 2D E 𝐸 E -representation contains a pair of degenerate gap functions.

	$u$	$g$
$A_1$	$\sin k_x\sin k_y{\tau }_3{\lambda }_{0,8}$	$\cos k_x\cos k_y{\tau }_0{\lambda }_{0,8}$
	$\cos k_x\cos k_y{\tau }_3{\lambda }_3$	$\sin k_x\sin k_y{\tau }_0{\lambda }_3$
	$i{\tau }_0\left(\sin k_{+}{\lambda }_5+\sin k_{-}{\lambda }_7\right)$	$i{\tau }_3\left(\sin k_{+}{\lambda }_5+\sin k_{-}{\lambda }_7\right)$
$A_2$	$\cos k_x\cos k_y{\tau }_3{\lambda }_1$	$\sin k_x\sin k_y{\tau }_0{\lambda }_1$
	$i{\tau }_0\left(\sin k_{+}{\lambda }_7-\sin k_{-}{\lambda }_5\right)$	$i{\tau }_3\left(\sin k_{+}{\lambda }_7-\sin k_{-}{\lambda }_5\right)$
$B_1$	$\sin k_x\sin k_y{\tau }_3{\lambda }_1$	$\cos k_x\cos k_y{\tau }_0{\lambda }_1$
	$i{\tau }_0\left(\sin k_{+}{\lambda }_7+\sin k_{-}{\lambda }_5\right)$	$i{\tau }_3\left(\sin k_{+}{\lambda }_7+\sin k_{-}{\lambda }_5\right)$
$B_2$	$\cos k_x\cos k_y{\tau }_3{\lambda }_{0,8}$	$\sin k_x\sin k_y{\tau }_0{\lambda }_{0,8}$
	$\sin k_x\sin k_y{\tau }_3{\lambda }_3$	$\cos k_x\cos k_y{\tau }_0{\lambda }_3$
	$i{\tau }_0\left(\sin k_{+}{\lambda }_5-\sin k_{-}{\lambda }_7\right)$	$i{\tau }_3\left(\sin k_{+}{\lambda }_5-\sin k_{-}{\lambda }_7\right)$
$E$	$(\sin k_{+},\sin k_{-})i{\tau }_0{\lambda }_2$	$(\sin k_{+},\sin k_{-})i{\tau }_3{\lambda }_2$
	${\tau }_3(\cos k_{+}{\lambda }_4,\cos k_{-}{\lambda }_6)$	${\tau }_0(\cos k_{+}{\lambda }_4,\cos k_{-}{\lambda }_6)$
	${\tau }_3(\cos k_{+}{\lambda }_6,\cos k_{-}{\lambda }_4)$	${\tau }_0(\cos k_{+}{\lambda }_6,\cos k_{-}{\lambda }_4)$

Table 4. Table 4: Classifications of the singlet superconducting order parameters for the NN pairings. λ ± 46 = λ 4 ± λ 6 superscript subscript 𝜆 plus-or-minus 46 plus-or-minus subscript 𝜆 4 subscript 𝜆 6 \lambda_{\pm}^{46}=\lambda_{4}\pm\lambda_{6} and λ ± 57 = λ 5 ± λ 7 superscript subscript 𝜆 plus-or-minus 57 plus-or-minus subscript 𝜆 5 subscript 𝜆 7 \lambda_{\pm}^{57}=\lambda_{5}\pm\lambda_{7} .

	$u$	$g$
$A_1$	$\left(\cos k_x-\cos k_y\right){\tau }_2{\lambda }_2$	$\left(\cos k_x+\cos k_y\right){\tau }_1{\lambda }_{0,8}$
		$\left(\cos k_x-\cos k_y\right){\tau }_1{\lambda }_1$
	$i{\tau }_1\left(\sin k_x{\lambda }_{+}^{57}+\sin k_y{\lambda }_{-}^{57}\right)$	$i{\tau }_2\left(\sin k_x{\lambda }_{-}^{46}+\sin k_y{\lambda }_{+}^{46}\right)$
$A_2$		$\left(\cos k_x-\cos k_y\right){\tau }_1{\lambda }_3$
	$i{\tau }_1\left(\sin k_x{\lambda }_{-}^{57}+\sin k_y{\lambda }_{+}^{57}\right)$	$i{\tau }_2\left(\sin k_x{\lambda }_{+}^{46}+\sin k_y{\lambda }_{-}^{46}\right)$
$B_1$	$\left(\cos k_x+\cos k_y\right){\tau }_2{\lambda }_2$	$\left(\cos k_x-\cos k_y\right){\tau }_1{\lambda }_{0,8}$
		$\left(\cos k_x+\cos k_y\right){\tau }_1{\lambda }_1$
	$i{\tau }_1\left(\sin k_x{\lambda }_{+}^{57}-\sin k_y{\lambda }_{-}^{57}\right)$	${\tau }_2\left(\sin k_x{\lambda }_{-}^{46}-\sin k_y{\lambda }_{+}^{46}\right)$
$B_2$		$\left(\cos k_x+\cos k_y\right){\tau }_1{\lambda }_3$
	$i{\tau }_1\left(\sin k_x{\lambda }_{-}^{57}-\sin k_y{\lambda }_{+}^{57}\right)$	$i{\tau }_2\left(\sin k_x{\lambda }_{+}^{46}-\sin k_y{\lambda }_{-}^{46}\right)$
$E$	$(\sin k_x,\sin k_y)i{\tau }_1{\lambda }_2$	$(\sin k_x,\sin k_y)i{\tau }_2{\lambda }_{0,8}$
		$(\sin k_x,\sin k_y)i{\tau }_2{\lambda }_{1,3}$
	$i{\tau }_2(\cos k_x{\lambda }_{-}^{57},\cos k_y{\lambda }_{+}^{57})$	${\tau }_2(\cos k_x{\lambda }_{-}^{46},\cos k_y{\lambda }_{+}^{46})$
	$i{\tau }_2(\cos k_x{\lambda }_{+}^{57},\cos k_y{\lambda }_{-}^{57})$	$i{\tau }_2(\cos k_x{\lambda }_{+}^{46},\cos k_y{\lambda }_{-}^{46})$

Equations

$$F_4=\beta |{\Delta }_1{|}^2|{\Delta }_2{|}^2+{\beta }^{\prime }({\Delta }_1^{\ast ,2}{\Delta }_2^2+c.c).$$

$${\mathcal{F}}_M=\alpha m_z^2+i\gamma m_z\left({\Delta }_1{\Delta }_2^{\ast }-{\Delta }_1^{\ast }{\Delta }_2\right),$$

$$\Delta =\sum _{\mathbf{k}}{\tau }_{ij}{f}^l(\mathbf{k})M_{ab}^l{\Delta }_{bj,ai}(\mathbf{k}),$$

$$\displaystyle\Delta_{1}:\cos k_{x}\cos k_{y}\tau_{0}\lambda_{1},\ \ \,\Delta_{2}:(\cos k_{x}+\cos k_{y})\tau_{1}\lambda_{3},\ \ \$$

$${\Delta }_1:(\cos k_x-\cos k_y){\tau }_1{\lambda }_3,{\Delta }_2:(\cos k_x-\cos k_y){\tau }_1{\lambda }_1.$$

$${\Delta }_1:\cos (k_x+k_y){\tau }_0{\lambda }_4,{\Delta }_2:\cos (k_x-k_y){\tau }_0{\lambda }_6,$$

$${\mathcal{H}}_0$$

$$\gamma$$

$$\gamma \approx -\frac{7\zeta (3)}{4{\pi }^3}\frac{{\lambda }_{soc}{N}_0}{(k_B{T}_c{)}^2},$$

$$T,C_4^1(z)T,C_2(z)T,C_4^3(z)T,$$

$${\sigma }_{x}T,{\sigma }_{y}T,{\sigma }_{x^{\prime }}T,{\sigma }_{y^{\prime }}T,$$

$$IT,S_4^{3}T,g({\sigma }_h,\tau )T,S_4^{1}T,$$

$$C_2(y)T,C_2(x)T,g(C_2(y^{\prime \prime }),{\vec{\tau }}_y)T,g(C_2(x^{\prime \prime }),{\vec{\tau }}_x)T,$$

$${\lambda }_1$$

$${\lambda }_3$$

$${\lambda }_5$$

$${\lambda }_7$$

$$H_0$$

$${\mathcal{H}}_0$$

$${\mathcal{H}}_{NNN}$$

$${\mathcal{H}}_{NN}$$

$$+2t_{xy}^{nn}(\cos k_x+\cos k_y)\left({\lambda }_0-\sqrt{3}{\lambda }_8\right)/3$$

$$+2t_2^{nn}(\cos k_x-\cos k_y){\lambda }_1$$

$$-2t_3^{nn}[(\sin k_x+\sin k_y){\lambda }_5+(\sin k_x-\sin k_y){\lambda }_7],$$

$${\mathcal{H}}_{soc}={\lambda }_{soc}{\tau }_0\left(\vec{L}\cdot \vec{\sigma }\right).$$

$${\mathcal{H}}_0={\mathcal{H}}_{NNN}+{\mathcal{H}}_{NN}.$$

$${\bar{\mathcal{H}}}_0$$

$${\bar{\mathcal{H}}}_{NNN}$$

$${\bar{\mathcal{H}}}_{NN}$$

$${\bar{\mathcal{H}}}_0={\bar{\mathcal{H}}}_{NNN}+{\bar{\mathcal{H}}}_{NN}+{\bar{\mathcal{H}}}_{soc},$$

$$G_e$$

$${\sigma }_y{G}_h{\sigma }_y$$

$$\displaystyle\begin{split}\mathcal{P}_{\pm}&=\frac{1}{2}\Big{\{}1\pm\big{(}2\bar{t}_{2}\sin k_{x}\sin k_{y}\bar{\lambda}_{3}\\ &+2\bar{t}_{7}(\cos k_{x}-\cos k_{y})\bar{\lambda}_{1}+\lambda_{soc}\bar{\lambda}_{2}\sigma_{z}\big{)}/E_{k}\Big{\}},\end{split}$$

$$E_k$$

$${ϵ}_{\pm }$$

$${\Delta }_1:(\cos k_x-\cos k_y){\tau }_1{\lambda }_1,{\Delta }_2:(\cos k_x-\cos k_y){\tau }_1{\lambda }_3.$$

Full text

The discovery of iron-based superconductors Kamihara et al. (2008) opened a new direction in the study on unconventional superconductivity Sigrist (1991); Van Harlingen (1995); Tsuei and Kirtley (2000). Significant progress has subsequently been made in searching for new superconductors Chen et al. (2008a, b); Ren et al. (2008); Rotter et al. (2008); Sasmal et al. (2008); Ni et al. (2010); Nakai et al. (2010); Wang et al. (2008); Borisenko et al. (2010); Mizuguchi et al. (2008); Hsu et al. (2008); Fang et al. (2008); Chen et al. (2009a); Guo et al. (2010); Qian et al. (2011), and their pairing mechanisms have attracted considerable attention Chubukov (2012); Wang and Lee (2011); Chen et al. (2014); Johnson et al. (2015a); Si et al. (2016). The parent compounds are metals with multiple Fermi surfaces around both $\Gamma$- and $M$- points. The possibility of the fully gapped extended $s_{\pm}$-wave superconducting gap function is supported by various experimental evidence and theoretical calculations Mazin et al. (2008); Kuroki et al. (2008); Parker et al. (2008); Chubukov et al. (2008); Bang et al. (2009); Seo et al. (2008); Chen et al. (2009b, b); Wang et al. (2009); Hu et al. (2015). On the other hand, several theoretical studies suggest that $s_{\pm}$ and $d_{x^{2}-y^{2}}$ pairings are nearly degenerate in the iron-pnictide superconductors Kuroki et al. (2008); Graser et al. (2009), leading to the possibility of a novel time-reversal (TR) symmetry breaking pairing $s_{\pm}+id_{x^{2}-y^{2}}$ Lee et al. (2009). It was proposed that a resonance mode carrying the $B_{1g}$-symmetry Lee et al. (2009); Scalapino and Devereaux (2009), which can be detected via Raman spectroscopy Böhm et al. (2014), exists if the $s_{\pm}$ and $d_{x^{2}-y^{2}}$ pairings are nearly degenerate. TR symmetry breaking pairing also naturally arises in mixed singlet and triplet pairing states Wu and Hirsch (2010); Hinojosa et al. (2014). A chiral $d+id$ pairing state is also found to spontaneously generate the gaps of the Haldane model Brydon et al. (2019).

Recently, the topological band structure of iron-based superconductors has aroused a great deal of attention. The FeSe1-xTex family with a wide range of composition $x$ is of particular interest Hsu et al. (2008); Fang et al. (2008); Yeh et al. (2008); Sales et al. (2009); Li (2009); LCC Ambolode et al. (2015); Wang et al. (2015); Wu et al. (2016); Johnson et al. (2015b); Rameau et al. (2019). Recent evidence shows that FeSe0.45Te0.55 is a strong topological insulator exhibiting a single Dirac cone on the (001) surface Zhang et al. (2018). With lowering the temperature below $T_{c}=14.5K$, both bulk and the surface become superconducting Miao et al. (2012); Yin et al. (2015). The surface superconductivity is predicted to be topologically non-trivial Wu et al. (2016); Wang et al. (2015); Xu et al. (2016). Excitingly, in the same system, Majorana zero modes in vortex cores have been observed Wang et al. (2018); Kong and Ding (2018) exhibiting the signature of spin-selective Andreev reflection Sun et al. (2016); Hu et al. (2016); He et al. (2014). Similar evidence to vortex core Majorana modes is also observed in (Li0.84Fe0.16)OHFeSe system Liu et al. (2018).

However, the non-trivial topology of the surface superconductivity in FeSe0.5Te0.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 Fu and Kane (2008). 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 Van Harlingen (1995); Tsuei and Kirtley (2000). In contrast, in a very recent laser-based angle-resolved photoemission spectroscopy (APERS) experiment on FeSe0.3Te0.7 Zaki et al. (2019), 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$ Lee et al. (2009); Platt et al. (2012); Khodas and Chubukov (2012); Fernandes and Millis (2013); Kang et al. (2018) and $s+is$ Stanev and Tešanović (2010); Marciani et al. (2013); Maiti and Chubukov (2013); Ahn et al. (2014) 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. Note that the degenerate Dirac cone in the surface state is protected by TR symmetry in the normal state due to the non-trivial band structure topology. Nevertheless, if the superconducting state breaks TR symmetry, the degeneracy of the Dirac cone is no longer protected in the superconducting state. 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)}\pm iA_{2g(u)}$, $B_{1g(u)}\pm iB_{2g(u)}$, and $E_{g(u)}\pm iE_{g(u)}$. In the normal state sufficiently close to $T_{c}$, the relative phase between two gap functions can still be locked at $\pm\frac{\pi}{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 FeSe1-xTex 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 $\Gamma$-point, as shown in Fig. 1 ($a$), described by an effective $k\cdot p$ Hamiltonian $H_{sf}=v(k_{x}\sigma_{y}-k_{y}\sigma_{x})-\mu$, where $\sigma$’s are Pauli matrices defined for the Kramers doublet at the $\Gamma$-point, and $\mu$ is the chemical potential. A superconducting gap $\Delta$ by itself, i.e., $H_{\Gamma}=-\mu+\Delta w_{-}+\Delta^{*}w_{+}$ where $w_{\pm}=w_{x}\pm iw_{y}$ are Pauli matrices in the Nambu space of the particle-particle channel, does not lift the degeneracy. To split the degeneracy, a mass term breaking TR symmetry is necessary, i.e., $\Delta H_{\Gamma}=-m_{z}\sigma_{z}$, and the associated Bogoliubov spectra become $\pm m_{z}+\sqrt{|\Delta|^{2}+\mu^{2}}$ around the $\Gamma$-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 $\sigma_{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.

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 Volovik and Gor’kov (1985); Sigrist (1991, 2005); Lee et al. (2009); Wu and Hirsch (2010), it has been shown that the mixing between two gap functions $\Delta_{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. $\Delta_{1,2}$ cannot form a symmetry invariant at the quadratic level, but they do at the quartic level via

[TABLE]

The $\beta^{\prime}$-term locks the relative phase between two gap functions, which equals $\beta^{\prime}|\Delta_{1}|^{2}|\Delta_{2}|^{2}\cos 2\Delta\varphi$, where $\Delta\varphi=\varphi_{1}-\varphi_{2}$, and $\varphi_{1,2}$ are the phases of two gap functions. When $\beta^{\prime}>0$, $\Delta\varphi$ is pinned at $\pm\frac{\pi}{2}$, giving rise to complex gap functions $\Delta_{1}\pm i\Delta_{2}$, which break TR symmetry spontaneously. This formalism also applies to the case that $\Delta_{1,2}$ form a two-dimensional (2D) irreducible representation. The complex gap functions $\Delta_{1}\pm i\Delta_{2}$ distribute more evenly over the Fermi surface than the real ones $\Delta_{1}\pm\Delta_{2}$, hence, they are energetically more preferable at the mean-field level Wu and Hirsch (2010). The corresponding Cooper pairs carry non-zero orbital moments, which could generate magnetic fields at boundaries as shown earlier Volovik and Gor’kov (1985); Sigrist (2005). 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 $\Delta_{1,2}$ through a cubic term as

[TABLE]

This satisfies both the $U(1)$ and TR symmetry. $\alpha>0$ is assumed in Eq. (2), and hence there is no spontaneous magnetic ordering by itself, rather, the magnetization is induced, $m_{z}=\frac{\gamma}{\alpha}|\Delta_{1}\Delta_{2}|\sin\Delta\varphi$, i.e, by coupling to the TR breaking superconducting orders. The sign of $m_{z}$ is determined by the relative phase $\Delta\varphi$ between $\Delta_{1,2}$.

The free energy density ${\cal 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 $\sqrt{2}\times\sqrt{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 $P4/nmm$ Fang et al. (2008); Li (2009); Miao et al. (2012); Yin et al. (2015), which is reviewed in the Supplemental Material (S.M.) I sup . 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\sim 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 FeSexTe1-x, the distribution of Se and Te breaks the $P4/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 $\sigma_{x},\sigma_{y}$, and are along the diagonal lines $x^{\prime}$ and $y^{\prime}$ denoted as $\sigma_{x^{\prime}}$ and $\sigma_{y^{\prime}}$, respectively.

We consider the order parameter properties under the crystalline symmetries. In multi-orbital systems, the superconducting gap function is expressed as

[TABLE]

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; $\tau$ is a $2\times 2$ matrix representing the sublattice channel; $f^{l}(\mathbf{k})$ is the angular form factor of momentum $\mathbf{k}$, and $M^{l}$ is the pairing matrix in the orbital channel; $l$ is the index for multiple combinations between $f^{l}(\mathbf{k})$ and $M^{l}$. The pairing matrix is defined as $\Delta_{bj,ai}(\mathbf{k})=\sum_{\mathbf{k}}i\sigma_{y,\alpha\beta}\langle|c^{\dagger}_{\mathbf{k},\alpha bj}c^{\dagger}_{-\mathbf{k},\beta ai}|\rangle$, where $i\sigma_{y}$ projects out the singlet pairing with Greek indices representing spin components, and $\langle||\rangle$ represents averaging over the thermal equilibrium state.

As required by Fermi statistics, for spin-singlet pairing, the product $\tau_{ij}f(\mathbf{k})M_{ab}$ in Eq. (3) needs to be even under the combined operations of $\mathbf{k}\to-\mathbf{k}$ and the transposes of $M$ and $\tau$. $f(\mathbf{k})$ can be an even function taking the forms of $1$, $\cos k_{x}\pm\cos k_{y}$, $\cos k_{x}\cos k_{y}$, $\sin k_{x}\sin k_{y}$, $\cos(k_{x}\pm k_{y})$. or, an odd function among $\sin k_{x}$, $\sin k_{y}$ and $\sin(k_{x}\pm k_{y})$. We choose the three $t_{2g}$-orbital bases, $d_{x^{\prime}z}$, $d_{y^{\prime}z}$ and $d_{xy}$, where $d_{x^{\prime}z}(d_{y^{\prime}z})$ extends along the diagonal $x^{\prime}(y^{\prime})$ direction as depicted in S. M. II sup . Hence, $M$ is a $3\times 3$ Hermitian matrix which is expanded in terms of the Gell-mann matrices $\lambda_{i}(i=1\sim 8)$ under the basis in the sequence of $(d_{x^{\prime}z},d_{y^{\prime}z},d_{xy})$, and the $3\times 3$ identity matrix $\lambda_{0}$, whose expressions are presented in S. M. II sup .

The representation of a gap function under the crystalline symmetry group is determined by the symmetry properties of $f(\mathbf{k})$, $M$, and $\tau$ as analyzed and presented in the S. M. II sup . 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 $\sigma_{x(y)}$ and $\sigma_{x^{\prime}(y^{\prime})}$ 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 sup . Different combinations of $f(\mathbf{k})$, $M$, and $\tau$ 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 Nica et al. (2017); Hu et al. (2018); Nakayama et al. (2018); König and Coleman (2019), and their importance have been analyzed in recent experiments Sprau et al. (2017); Kostin et al. (2018). 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(\mathbf{k})$. For example, consider the following two gap functions with even parity,

[TABLE]

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^{\prime}z}$ and $d_{y^{\prime}z}$-orbitals as shown in Fig. 3 ($a$) and ($b$). Since $d_{xz}\to d_{yz}$ and $d_{yz}\to-d_{xz}$ under the $90^{\circ}$ rotation, $\lambda_{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:

[TABLE]

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,

[TABLE]

The former (latter) describes the pairing between the $d_{x^{\prime}z}$-orbital ($d_{y^{\prime}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 $\Delta_{1}$ and $\Delta_{2}$ should contain that of $m_{z}$, i.e., $A_{2g}$. This yields the following possibilities of pairing symmetries: $B_{1g(u)}\pm iB_{2g(u)}$, $A_{1g(u)}\pm iA_{2g(u)}$ and $E_{g(u)}\pm iE_{g(u)}$. Examples of the above pairing symmetries with even parity are provided in Eqs. (4), (5), and (6). They carry the same symmetry under rotation, and hence, the experimental observables do not break the rotational symmetry. However, they mix different symmetries with respect to the vertical reflection plane, and such symmetries are also spontaneously breaking. For example, the magnetization $m_{z}$ is odd under such reflections.

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 $\Delta_{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 FeTe1-xSex around the $\Gamma$-point, which consists of the $t_{2g}$-orbitals $d_{x^{\prime}z}$, $d_{y^{\prime}z}$ and $d_{xy}$ Wang et al. (2015); Xu et al. (2016). Neglecting the small dispersion along the $z$-axis, the 3-band tight-binding Hamiltonian is expressed as Daghofer et al. (2010), $H_{0}=\sum_{k}\psi^{\dagger}(k){\cal H}_{0}(k)\psi(k)$ where $\psi=\left[d_{x^{\prime}z,\uparrow},d_{y^{\prime}z,\uparrow},d_{xy,\uparrow},d_{x^{\prime}z,\downarrow},d_{y^{\prime}z,\downarrow},d_{xy,\downarrow}\right]^{T}$ and the matrix kernel ${\cal H}_{0}(k)$ is given by,

[TABLE]

where $\mathcal{H}_{NNN}$ and $\mathcal{H}_{NN}$ represent the NNN and NN hoppings, respectively, with detailed forms presented in S. M. III sup . $\mathcal{H}_{soc}=\lambda_{soc}\tau_{0}\vec{L}\cdot\vec{\sigma}$ is the atomic spin-orbit coupling with $\lambda_{soc}$ the coupling strength and $\vec{L}$ representing the onsite orbital angular momentum projected to the $t_{2g}$-basis. Explicitly, $\vec{L}=\left((\lambda_{5}-\lambda_{7})/\sqrt{2},-(\lambda_{5}+\lambda_{7})/\sqrt{2},-\lambda_{2}\right)$.

Based on the band structure Eq. (7), the coupling coefficient $\gamma$ in Eq. (2) can be evaluated as

[TABLE]

where $\beta=1/(k_{B}T)$ is the inverse of temperature; $f^{l_{i}}_{i}$ and $M_{i}^{l_{i}}$ with $i=1,2$ are the angular form factors and orbital pairing matrix kernels of $\Delta_{1,2}$, respectively; $G_{e}(\mathbf{k},i\omega_{n})=(i\omega_{n}-\mathcal{H}_{0}(\mathbf{k}))^{-1}$ is the Matsubara Green’s function with $\omega_{n}=(2n+1)\pi/\beta$, and $G_{h}(\mathbf{k},i\omega_{n})=G_{e}^{*}(-\mathbf{k},-i\omega_{n})$. As shown in S. M. III sup , the hole-like Fermi pockets around the $\Gamma$-point are mainly from the bonding states between two Fe-sublattices, i.e., they are approximately eigenstates of $\tau_{1}$ with the eigenvalue of 1. Hence, only tracing over the spin and orbital channels are needed, and only these gap functions characterized by $\tau_{0,1}$ are considered. Gap functions with $\tau_{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 $\Delta_{1,2}$ leading to the spontaneous magnetization $m_{z}$. We begin with the cases of $B_{1g(u)}\pm iB_{2g(u)}$. For parity even, i.e., $B_{1g}\pm iB_{2g}$, we take $\Delta_{1,2}$ in the form of Eq. (4). As shown in S. M. IV sup , after further reducing the band Hamiltonian Eq. (7) to a two-band model only based on $d_{x^{\prime}(y^{\prime})z}$, Eq. (8) yields an analytic expression as

[TABLE]

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 $\gamma$ 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 sup ). The parity odd case, i.e., $B_{1u}\pm iB_{2u}$, is also numerically checked to yield a nonzero $\gamma$, for example, with $\Delta_{1}:\tau_{0}(\sin(k_{x}+k_{y})\lambda_{7}+\sin(k_{x}-k_{y})\lambda_{5})$ and $\Delta_{2}:\tau_{0}(\sin(k_{x}+k_{y})\lambda_{5}-\sin(k_{x}-k_{y})\lambda_{7})$. The above two cases also break the mirror symmetries of $\sigma_{x(y)}$ and $\sigma_{x^{\prime}(y^{\prime})}$ spontaneously. They are topologically non-trivial belonging to the C-class supporting the chiral Majorana edge modes Hasan and Kane (2010); Qi and Zhang (2011); Bansil et al. (2016); Chiu et al. (2016). We have also studied both cases of $A_{1g(u)}\pm iA_{2g(u)}$, which also yield non-zero $\gamma$’s as shown in S. M. V sup . The nodal pairing gap functions presented in Eq. (5) are used for the even parity case, and the nodeless pairing with $\Delta_{1}:\tau_{0}(\sin(k_{x}+k_{y})\lambda_{5}+\sin(k_{x}-k_{y})\lambda_{7})$ and $\Delta_{2}:\tau_{0}(\sin(k_{x}+k_{y})\lambda_{7}-\sin(k_{x}-k_{y})\lambda_{5})$ are used for the odd parity case. For the case of $E_{g(u)}\pm iE_{g(u)}$, we take the gap functions presented in Eq. (6) as an example. Since $E_{g(u)}\otimes E_{g(u)}=A_{1g}\oplus A_{2g}\oplus B_{1g}\oplus B_{2g}$, it also yields a nonzero $\gamma$ as calculated in S. M. IV sup , and consequently induces magnetization.

In strongly correlated superconductors, there exist strong superconducting phase fluctuations in the normal state close to $T_{c}$ Emery and Kivelson (1995). In this case, the phases $\varphi_{1}$ and $\varphi_{2}$ of gap functions are disordered such that $\langle\Delta_{1,2}\rangle=0$, but $|\Delta_{1,2}|$ remains finite. The $\beta^{\prime}$ term in Eq. (1) can still pin the relative phase $\Delta\varphi=\pm\frac{\pi}{2}$. This transition breaks TR symmetry and its critical temperature $T^{\prime}>T_{c}$. We still expect a weaker but still finite $m_{z}$ in the temperature window between $T_{c}$ and $T^{\prime}$.

The bulk spin magnetization induced by TR breaking pairing qualitatively explains the gap opening of the surface Dirac cone observed in the FeSe0.3Te0.7 superconductor Zaki et al. (2019). 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 sup ), which is still smaller than the lower critical magnetic field of the FeSe1-xTex superconductors ($H_{c1}\sim 30$ Gauss)Yadav and Paulose (2009); Abdel-Hafiez et al. (2013). 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 Shopova and Uzunov (2005); Chirolli et al. (2017). This is consistent with the observation that neutron spectroscopy does not detect a bulk magnetic field Xu et al. (2018). A very weak but finite internal magnetic field around $0.15$ Gauss is detected in the FeSe superconductor by the muon spin rotation measurement ($\mu$SR)Watashige et al. (2015); Matsuura et al. (2019), 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 $\pm\frac{\pi}{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 FeTe0.7Se0.3 Zaki et al. (2019).

Note added. During the review process, we become aware of a recent work by T. Kawakami Kawakami and Sato (2019), where they mainly discuss the topological odd-parity superconductivity in iron-based superconductors.

Acknowledgments–. L. H., P. D. J., and C. W. acknowledge helpful discussions with J. Tranquada, G. D. Gu. L. H. and C. W. are supported by AFOSR FA9550-14-1-0168. P. D. J. is supported by the Office of Science, U.S. Department of Energy under contract number DE-SC001274.

Appendix A P4/nmm non-symmorphic group

Here we review the $P4/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 $\left\{E,S_{4},C^{\prime}_{2}(z),S_{4}^{3},C_{2}(x),C_{2}(y),\sigma_{x^{\prime}},\sigma_{y^{\prime}}\right\}$, where $C^{\prime}_{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 $\sqrt{2}a$ 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}=\left\{E,C_{4}(z),C_{2}(z),C^{3}_{4}(z),\sigma_{x},\sigma_{y},\sigma_{x^{\prime}},\sigma_{y^{\prime}}\right\}$, which is used in this work to specify the irreducible representations of gap functions. The configurations of the reflection planes of $\sigma_{x},\sigma_{y},\sigma_{x^{\prime}}$ and $\sigma_{y^{\prime}}$ are given in Fig. 2 in the main text.

The $P4/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 $\sum_{i=1,2,3}l_{i}\vec{a}_{i}$ with integers $l_{1,2,3}$, and the lattice vectors are defined as $\vec{a}_{1}=(a,a,0)$, $\vec{a}_{2}=(a,-a,0)$ and $\vec{a}_{3}=(0,0,c)$ with $c$ the length of Fe-Fe bond along the $z$-direction. There are 8 cosets generated by $g\in C_{4v}$ point group:

[TABLE]

where $C_{4}^{1,3}(z)$ and $C_{2}(z)$ are rotations around the $z$-axis passing the Se atom, and $\sigma_{x(y)}$ and $\sigma_{x^{\prime}(y^{\prime})}$ 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

[TABLE]

where $S_{4}$ and $S_{4}^{3}$ 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(\sigma_{h},\vec{\tau})$ is the glide reflection with $\sigma_{h}$ the reflection with respect to the $xy$-plane and $\vec{\tau}=a(1,0,0)$; $g(C_{2}(x^{\prime\prime}),\vec{\tau}_{x})$ and $g(C_{2}(y^{\prime\prime}),\vec{\tau}_{y})$ are screw rotations with their axes $x^{\prime\prime}$ and $y^{\prime\prime}$ along the $45^{\circ}$ and $135^{\circ}$ degrees passing the Fe-Fe bond centers; $\vec{\tau}_{x}=\frac{1}{2}(-a,a,0)$ and $\vec{\tau}_{y}=\frac{1}{2}(a,a,0)$.

The character table of $C_{4v}$ is given in Table. 1, 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 $P4/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.

Appendix B The gap function symmetry

According to the crystal symmetry, the $d_{x^{\prime}z},d_{y^{\prime}z},d_{xy}$-orbitals on the A and B sublattices of Fe are constructed shown in Fig. 4, where the $d_{x^{\prime}z}$ and $d_{y^{\prime}z}$ are along the diagonal Fe-Fe lines. For the three dimensional orbital space, we define the $3\times 3$ Gell’mann matrices according to the basis of $(d_{x^{\prime}z},d_{y^{\prime}z},d_{xy})^{T}$ as

[TABLE]

In addition, $\lambda_{0}$ is the $3\times 3$ identity matrix. We also denote the Pauli matrices in the channel of the Fe A/B sublattice as $\tau_{1,2,3}$ and $\tau_{0}$ the $2\times 2$ identity matrix. Similarly, the Pauli matrices in the spin channel is denoted as $\sigma_{x,y,z}$, and $\sigma_{0}$ is also the $2\times 2$ identity matrix.

The classification for matrices in the sublattice channel ($\tau$), the orbital channel ($\lambda$), and the spin channel ($\sigma$) according to the $C_{4v}$ and inversion symmetries are summarized in the following Table. 2.

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.

Appendix C Three-orbital model

The corresponding tight-binding model can be constructed below,

[TABLE]

where $\psi(\mathbf{k})=[d_{A,x^{\prime}z,\uparrow},d_{A,y^{\prime}z,\uparrow},d_{A,xy,\uparrow},d_{B,x^{\prime}z,\uparrow},d_{B,y^{\prime}z,\uparrow},\\ d_{B,xy,\uparrow},d_{A,x^{\prime}z,\downarrow},d_{A,y^{\prime}z,\downarrow},d_{A,xy,\downarrow},d_{B,x^{\prime}z,\downarrow},d_{B,y^{\prime}z,\downarrow},d_{B,xy,\downarrow}\\ ]^{T}$, and the matrix kernel is defined as

[TABLE]

The $\mathcal{H}_{NNN}$ term describes the next-nearest neighbor (NNN) hopping Hamiltonian,

[TABLE]

where $\Delta_{xy}$ is the the on-site energy difference between the $d_{xy}$-orbital and the two degenerated $d_{x^{\prime}z/y^{\prime}z}$; the hopping integrals $t^{nnn}_{\sigma,\pi,\eta,xy}$ here describe the bonding along the NNN Fe-Fe bond (see Fig. 5(a-d)). $t^{nnn}_{\sigma}$ and $t^{nnn}_{\pi}$ are the $\sigma$ and $\pi$-bonding strengths of the $d_{x^{\prime}(y^{\prime}),z}$ orbitals; $t_{xy}^{nnn}$ is the bonding strength between two $d_{xy}$-orbitals; $t_{\eta}^{nnn}$ is the bonding between $d_{x^{\prime}(y^{\prime})z}$ and $d_{xy}$-orbitals.

The $\mathcal{H}_{NN}$ term represents the nearest neighbor (NN) hopping Hamiltonian,

[TABLE]

where $t^{nn}_{1,2,3,xy}$ describe the bonding along the NN Fe-Fe bond (see Fig. 5(e-h)). $t_{1}^{nn}$ and $t_{2}^{nn}$ describe the intra- and inter-orbital bonding between $d_{x^{\prime}z}$ and $d_{y^{\prime}z}$-orbitals; $t^{nn}_{xy}$ is the bonding between $d_{xy}$-orbitals; $t^{nn}_{3}$ is the bonding between $d_{xy}$ and $d_{x^{\prime}(y^{\prime})z}$ orbitals.

The last term $\mathcal{H}_{soc}$ is the spin-orbit coupling Hamiltonian,

[TABLE]

where $\vec{L}=\left((\lambda_{5}-\lambda_{7})/\sqrt{2},-(\lambda_{5}+\lambda_{7})/\sqrt{2},-\lambda_{2}\right)$.

Since $\tau_{1}$ is conserved in Eq.(41), we can label the eigenstates of ${\cal H}_{0}$ by $\tau_{1}$’s eigenvalues of $\pm 1$. There are mainly three hole-like pockets around the $\Gamma$-point, and all of them carry the eigenvalue of 1 of $\tau_{1}$. Projecting ${\cal H}_{0}$ into this sector, we arrive at a $6\times 6$ Hamiltonian matrix as

[TABLE]

The corresponding band structure is calculated along the high symmetry lines $\Gamma$-$X$-$M$-$\Gamma$ 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.

Appendix D The reduced two-orbital model

Due to the relatively large value of $\Delta_{xy}$, we can project out the $d_{xy}$-orbital and arrive a reduced two-orbital model with only $d_{x^{\prime}z}$ and $d_{y^{\prime}z}$-orbitals, whose matrix kernel read

[TABLE]

with

[TABLE]

where $\bar{t}^{nnn}_{\sigma,\pi}$ are hopping elements for NNN Fe-Fe bonding; $\bar{t}^{nn}_{1,2}$ are those for the NN Fe-Fe bonding; the $\bar{\lambda}_{1,2,3}$ are Pauli matrices defined in the orbital channel ($d_{x^{\prime}z}$ and $d_{y^{\prime}z}$), and $\bar{\lambda}_{0}$ is the $2\times 2$ identity matrix. Then there will be two hole-like pockets around the $\Gamma$-point and all the states on the Fermi surface are eigenstates of $\tau_{1}$ with the eigenvalues $+1$. Within this sector, the $\bar{\cal{H}}_{0}$ can be further simplified as

[TABLE]

where $\bar{\cal{H}}_{soc}=\lambda_{soc}\bar{\lambda}_{2}\sigma_{z}$. Then the Green’s function can be analytically solved as,

[TABLE]

where $\beta=1/T$ ($T$ is temperature) and the fermion Matsubara frequency $\omega_{n}=(2n+1)\pi/\beta$. The projection operators are defined as

[TABLE]

with

[TABLE]

where $\epsilon_{0}=2\bar{t}_{1}\cos k_{x}\cos k_{y}+2\bar{t}_{5}(\cos k_{x}+\cos k_{y})$. This simplification helps to calculate the Feymann diagram analytically discussed in the main text.

Appendix E Numerical calculations of $\gamma$

Based on the Hamiltonian in Eq. (45), 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 spin-orbit coupling $\lambda_{soc}$ region, $\gamma$ is linearly proportional to $\lambda_{soc}$. The five examples for $\Delta_{1,2}$ are studied with the following angular form factors and pairing matrices.

• 1.)

Example of $A_{1g}\pm iA_{2g}$ case:

[TABLE]

• 2.)

Example of $A_{1u}\pm iA_{2u}$ case:

[TABLE]

• 3.)

Example of $B_{1g}\pm iB_{2g}$ case:

[TABLE]

• 4.)

Example of $B_{1u}\pm iB_{2u}$ case:

[TABLE]

• 5.)

Example of $E_{g}\pm iE_{g}$ case:

[TABLE]

However, we find that for the $E_{u}\pm iE_{u}$ case, the trace operation in Eq. 4 in the main text already yields zero based on the above three-band model, yielding $\gamma=0$, although a non-zero $\gamma$ is allowed by symmetry.

Appendix F 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}$-pairingLee et al. (2009). 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 $\sigma_{x}$ and $\sigma_{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,\sigma\rangle,\}$ with $i=1\sim 4$ and $\sigma=\uparrow,\downarrow$ as follows,

[TABLE]

where $\phi^{A}$ and $\phi^{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\rangle_{\sigma}$ and $|4\rangle_{\sigma}$, and its Hamiltonian matrix is given by $H_{Neel}=\begin{pmatrix}H_{Neel,\uparrow}&0\\ 0&H_{Neel,\downarrow}\end{pmatrix}$, where

[TABLE]

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 two-fold degeneracy is protected by two reflection symmetries with respect to two perpendicular mirror planes, i.e., $\sigma_{x}$ and $\sigma_{y}$. When applying to half-integer spin fermions, actually $\sigma_{x}$ and $\sigma_{y}$ anti-commute with each other, giving rise to the protected double degeneracy.

Appendix G 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\mu_{B}$, where $\mu_{B}=e\hbar/2m_{e}c$ is the Bohr magneton, and $b$ is the relative magnetization, then the magnetization strength $M_{s}$ is given by,

[TABLE]

where $a_{0}=\hbar^{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,

[TABLE]

where $\Phi_{0}=hc/2e$ is the magnetic flux quantum and $\alpha$ is the fine-structure constant.

In the FeSe1-xTex 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_{s}^{3}$ is estimated as $680$. The upper bound of the relative magnetization $b$ can be estimated as the total hole density, which is $6\pi k_{F}^{2}/(2\pi)^{2}\approx 3\pi/200\approx 4.7\%$ with $k_{F}\approx\pi/10$ as the average Fermi momentum of the three hole pockets (see Fig. 6). Therefore,

[TABLE]

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.