Introduction

Optical studies have been very important in superconductivity research since the superconducting gap was first observed by far-infrared optical measurements1. Not only does the optical absorption gap directly reveal the superconducting gap size, but also the loss of spectral weight of the optical conductivity in the superconducting transition shows the superfluid density2,3,4. Optical responses in superconductors are well understood by the dirty-limit theory of Mattis and Bardeen5 and its extensions to arbitrary purity6,7. Impurity is essential in the Mattis–Bardeen theory because Bogoliubov quasiparticles cannot be excited by uniform light when momentum is conserved in the Bardeen–Cooper–Schrieffer (BCS) model8,9. In this paradigm, optical responses are due to impurity scattering and correspond to the Drude responses remaining in the superconducting state. They are thus completely described within a single-band model such as the BCS model and approaches the Drude formula as the photon energy increases above the gap [Fig. 1].

Fig. 1: Dirty and clean optical responses.
figure 1

a Band structure in the Bogoliubov de-Gennes (BdG) formalism. Superconducting pairing opens the bandgap at the Fermi level EF by mixing electron (red) and hole (blue) bands as shown in (b). ck and \({c}_{-{\bf{k}}}^{\dagger }\) indicate the electron and hole annihilation operators, respectively. b Optical excitations in dirty and clean limits by spatially uniform light (q = 0) across the superconducting gap 2Δ. The momentum transfer p − k in the left figure is supplied from the impurity potential. The crossover from dirty to clean optical responses occurs when the mean free path l exceeds the superconducting coherence length ξ0 by a factor \({x}_{c}={({k}_{F}{\xi }_{0})}^{2}{\alpha }^{-2}\). Here, kF is the Fermi wave number, and α is a degree of multiband pairing (see the text above Eq. (6) for its more precise definition). The red and blue colors schematically represent the electron–hole band mixing in the superconducting state. c Real part of the optical conductivity in clean systems (x > 1). σs = σdis + σint in the superconducting state is the sum of disorder-mediated and intrinsic parts. σint ∝ ω−1 is insensitive to x while σdis ∝ x−1ω−2 depends significantly on x. The Drude conductivity in the normal state σn is also shown as dashed lines. d, e Ratio of superconducting and normal state conductivities. d Disorder-mediated part. e Intrinsic part.

Full size image

On the other hand, there have been cumulative studies revealing the relevance of multiband effects in superconductivity. Strong gap anisotropy and multiple gap signatures due to orbital-dependent pairing have been observed in various superconductors, including elemental metals Nb, Ta, V, and Pb10,11, compound MgB212, strontium titanates13, iron pnictides and chalcogenides14,15,16, and heavy fermion compounds17,18. Multiband effects are also considered to be important in the superconductivity of strontium ruthenates19,20, some half-Heusler compounds with J = 3/2 degrees of freedom21,22,23, and twisted bilayer graphene24,25,26. This raises the question of whether multi-band effects can modify the optical responses. However, there has been no observation of a significant deviation from the Mattis–Bardeen theory in any materials.

In this work, we challenge the Mattis–Bardeen paradigm by showing that a significant portion of the far-infrared optical response in a clean multiband superconductor FeSe is due to intrinsic momentum-conserving excitations. We show that such intrinsic optical excitations are allowed by multiband effects. To establish the criteria for nonzero intrinsic responses systematically, we present a tenfold way classification of optical excitations as well as selection rules due to unitary symmetries. Here, the tenfold way classification is by three symmetry operations \({\mathfrak{T}}\), \({\mathfrak{C}}\), and S that leave momentum invariant27, whereas the original tenfold classification by Altland and Zirnbauer28,29 is by three spatially local symmetries, including time reversal T, particle–hole conjugation C, and chiral S symmetries [see Table 1]. Since T and C reverses the momentum, the combination of them with spatial inversion P (or any other momentum-reversing unitary operation) defines \({\mathfrak{T}}\) and \({\mathfrak{C}}\). In this classification, \({\mathfrak{C}}\) symmetry is the key player that imposes a new selection rule. We find that the absence of intrinsic optical excitations in single-band models can be attributed to \({\mathfrak{C}}\) symmetry in the superconducting state.

Table 1 Tenfold way classification of the lowest optical excitations in superconductors.
Full size table

As real materials are always accompanied by disorder, intrinsic responses coexist with disorder-mediated responses. A superconductor is considered to be clean when the mean free path l is larger than the superconducting coherence length ξ0 and to be dirty when l is smaller than ξ0. We show that the crossover from disorder-mediated to intrinsic optical responses occurs at a very clean regime \(l \sim {({k}_{F}{\xi }_{0})}^{2}{\alpha }^{-2}{\xi }_{0}\gg {\xi }_{0}\), as illustrated in Fig. 1b, where kF is the Fermi wave number, and 0 ≤ α ≤ 1 is a degree of multi-band pairing explained below. When l goes above this value, the optical conductivity follows the ω−1 behavior of the intrinsic response, deviating from the Drude-like ω−2 behavior [Fig. 1d, e]. We discuss the optical response of superconducting FeSe, which is closest to this crossover regime.

Results

Setting

Our theory is based on the mean-field theory of superconductors. We assume uniform illumination of light at zero temperature and the conservation of momentum. In momentum space, the single-particle mean-field Hamiltonian has the Bogoliubov-de Gennes (BdG) form

$$H({\bf{k}})=\left(\begin{array}{ll}h({\bf{k}})&{{\Delta }}({\bf{k}})\\ -{{{\Delta }}}^{* }(-{\bf{k}})&-{h}^{T}(-{\bf{k}})\end{array}\right)$$
(1)

in the basis of the Nambu spinor defined by \(\hat{{{\Psi }}}={({\hat{c}}_{\rho s{\bf{k}}},{\hat{c}}_{\rho s-{\bf{k}}}^{\dagger })}^{T}\), where \({\hat{c}}_{\rho s{\bf{k}}}\) is the electronic quasiparticle annihilation operator with orbital ρ and spin s = ↑,↓ indices. Here, h(k) is the normal-state Hamiltonian, and the pairing function \({{\Delta }}({\bf{k}})\propto \left\langle {c}_{{\bf{k}}}{c}_{-{\bf{k}}}\right\rangle\) satisfies Δ(k) = −ΔT(−k) due to Fermi statistics of electrons. The BdG Hamiltonian always has particle–hole symmetry CH(k)C−1 = −H(−k) under C = τxK, where τx is a Pauli matrix for the particle–hole indices and K is the complex conjugation operator.

The electromagnetic field couples to the normal-state Hamiltonian through the minimal coupling \({\bf{k}}\to {\bf{k}}+\frac{e}{\hslash }{\bf{A}}\), where q = −e (+e) for the electron (hole) sector. It follows that the velocity operator is

$${V}^{a}({\bf{k}})=\frac{1}{e}{\left.\frac{\partial H}{\partial {A}_{a}}\right|}_{{\bf{A}} = 0}=\frac{1}{\hslash }\left(\begin{array}{ll}{\partial }_{{k}_{a}}h({\bf{k}})&0\\ 0&{\partial }_{{k}_{a}}[{h}^{T}(-{\bf{k}})]\end{array}\right).$$
(2)

Matrix elements of this operator are important in our analysis because they describe the transition amplitudes. In the clean limit, the real part of the optical conductivity tensor is given by

$${\sigma }^{ca}(\omega )=\frac{\pi {e}^{2}}{2\hslash \omega }\int _{{\bf{k}}}\sum _{n,m}{f}_{nm}({\bf{k}}){V}_{nm}^{c}({\bf{k}}){V}_{mn}^{a}({\bf{k}})\delta (\omega -{\omega }_{mn}({\bf{k}})),$$
(3)

where ω is the frequency of light, fnm = fn − fm is the difference between the Fermi distribution of the nth band fn, \({V}_{mn}^{a}=\left\langle m| {V}^{a}| n\right\rangle\), and ωmn = ωm − ωn, where \(H\left|n\right\rangle =\hslash {\omega }_{n}\left|n\right\rangle\)9,30. The delta function is replaced by the Lorentzian distribution when the mean free path is finite.

Selection rules

Equation (3) is positive-semidefinite when c = a. Therefore, interband transitions are completely forbidden only when symmetries impose \({V}_{mn}^{a}({\bf{k}})=0\) at every k30. The relevant symmetry operators should be k-local (k → k). A unitary symmetry imposes selection rules by λm(k) = λV(k)λn(k), where λm,n and λV are symmetry eigenvalues of m,n states, and the velocity operator, respectively. We always have λV = 1 because k-local symmetry operations leave Va invariant, as one might expect because the velocity operator should transform like k. The selection rules thus simply become

$${V}_{mn}^{a}({\bf{k}})=0\quad \,\text{when}\,{\lambda }_{m}({\bf{k}})\;\ne \;{\lambda }_{n}({\bf{k}}),$$
(4)

meaning that optical excitations are forbidden between two different eigenspaces [Fig. 2a].

Fig. 2: Selection rules in clean superconductors.
figure 2

a Selection rule by unitary symmetry. λ1,2 are eigenvalues of a k-local unitary symmetry operator. No optical excitation occurs between two states with different eigenvalues. b Selection rule by \({\mathfrak{C}}\) symmetry. The case with \({\mathfrak{C}}={\rm{PC}}\) is shown. No optical transition occurs between PC-related states when (PC)2 = −1. c Optical excitation channels in \({\mathfrak{C}}\)-symmetric superconductors in the clean limit. At low photon energies comparable to the superconducting gap 2Δ, the relevant excitations are spectrum-inversion-symmetric (SIS) ones, i.e., from energy −E to E. For nondegenerate bands, they are transitions between \({\mathfrak{C}}\)-related pairs. d, e Optical excitations in spin-degenerate systems with and without spin–orbit coupling, respectively. (here, the order of d and e has been changed in order to match the label in the figure.) In d, \({\mathfrak{T}}={\rm{PT}}\) symmetry with (PT)2 = −1 imposes Kramers degeneracy. As a state \(\left|n\right\rangle\) can be excited to one of two SIS states, \({\rm{PC}}\left|n\right\rangle\) and \({\rm{PT}}({\rm{PC}}\left|n\right\rangle )={\rm{TC}}\left|n\right\rangle\), the excitation from \(\left|n\right\rangle\) is possible even when one transition channel, from \(\left|n\right\rangle\) to \({\rm{PC}}\left|n\right\rangle\), is blocked by (PC)2 = −1. The same applies to the excitation from \({\rm{PT}}\left|n\right\rangle\). In e, the boxes labeled up and down indicate the spin up and down eigenspace (sz = ℏ/2 and − ℏ/2, respectively). Since C reverses the spin (the anti-particle of a spin-up electron carries the down spin) while P does not change the spin, PC reverses the spin. Its combination with the spin rotation around the y axis, which is iσy for spin-singlet pairing, acts within a sz sector. For spin-triplet pairing, the spin rotation around the y-axis acts on the particle and hole sector with an opposite sign due to the spin carried by the Cooper pair, so the additional τz is introduced (see section 3 in the “Methods”). Optical excitations are forbidden when \({\mathfrak{C}}\) defined within a spin sector, which is −iPCσy for singlet pairing (−iPCτzσy for triplet pairing), satisfies \({{\mathfrak{C}}}^{2}=-1\). EF is the Fermi level in all figures.

Full size image

Let us consider the transition between two states in the same eigenspace of the unitary symmetry group. The remaining k-local symmetries come in three types: anti-unitary \({\mathfrak{T}}\), anti-unitary anti-symmetry \({\mathfrak{C}}\), and unitary anti-symmetry S, where anti-symmetry means that the operator anti-commutes with the Hamiltonian. They form ten EAZ symmetry classes27,28,29,31 shown in Table 1. \({\mathfrak{T}}\) (or \({\mathfrak{C}}\)) comes as a combination of T (or C) with a k-reversing unitary operator such as spatial inversion P in any dimensions or twofold rotation C2z in two dimensions. S is the combination \({\mathfrak{T}}{\mathfrak{C}}\) up to a phase factor. We find that only \({\mathfrak{C}}\)-type symmetry can additionally exclude transition channels within an eigenspace of the unitary symmetry group. By using that \({\mathfrak{C}}\) is anti-unitary and that the velocity operator is invariant under \({\mathfrak{C}}\) as shown in the “Methods” section 1, we have

$$\left\langle {\mathfrak{C}}\cdot n{\bf{k}}| {V}^{a}({\bf{k}})| n{\bf{k}}\right\rangle =0\quad \,\text{when}\,\ {{\mathfrak{C}}}^{2}=-1.$$
(5)

This constrains, in particular, the lowest-energy excitations, as illustrated in Fig. 2b, c. If bands are nondegenerate in each eigenspace, Eq. (5) indicates that the excitations across the gap are forbidden when \({{\mathfrak{C}}}^{2}=-1\) [Fig. 2b]. See class C and CI in Table 1.

We find that the absence of optical excitations in single-band metal models, described by a two-band BdG Hamiltonian, can be attributed to the existence of \({\mathfrak{C}}=i{\tau }_{y}K\) symmetry. A single-band metal has the \({\mathfrak{C}}\) symmetry in the superconducting state, independent of the pairing symmetry, when it has symmetry ξ(k) = ξ(k) in the normal state, where h(k) = ξ(k) is the 1 × 1 Hamiltonian. Since the formation of Cooper pairs at the Fermi level requires such a symmetry relating k and −k, it means that typically no optical excitations can occur in superconductors originating from single-band metals. One can extend this result to show the absence of optical excitations in multiband systems satisfying a generalized single-band pairing condition, the so-called zero superconducting fitness condition (see section 2 in the “Methods”).

We have three ways of generating nontrivial optical excitations in an eigenspace. When bands are nondegenerate within an eigenspace, one can (i) break \({\mathfrak{C}}\) symmetry (EAZ class A, AI, AII, and AIII) or (ii) realize \({{\mathfrak{C}}}^{2}=+1\) (class D and BDI). (iii) Or, when bands are Kramers degenerate due to \({\mathfrak{T}}\) symmetry satisfying \({{\mathfrak{T}}}^{2}=-1\), lowest-energy excitations are generally allowed irrespective of the sign of \({{\mathfrak{C}}}^{2}\) (class DIII and CII). The first condition (i) just means breaking inversion symmetry when other unitary symmetries do not exist, which was demonstrated in ref. 30. The second (ii) implies that the superconductor may host stable BFSs27,32. Since \({{\mathfrak{C}}}^{2}=+1\) protects 0D \({{\mathbb{Z}}}_{2}\) topological charges, \({{\mathbb{Z}}}_{2}\)-stable nodal surfaces/lines/points in 3D/2D/1D can appear after superconducting pairing on the Fermi surfaces, respectively, which we call as BFSs without distinguishing their dimension. Let us note that these are twofold degenerate BFSs. On the other hand, a stable nondegenerate BFS can appear in the EAZ classes A, AI, and AII. For instance, a superconductor with broken inversion and time-reversal symmetries can host stable BFSs33,34,35,36. Their stability is guaranteed by the change of the number of occupied BdG bands across the BFS, which is a \({\mathbb{Z}}\) topological charge. Since these classes correspond to the case (i), the symmetry protection of the stable BFSs, whether it is twofold degenerate or not, indicates that the lowest-energy optical excitations are possible. The last possibility (iii) is realized in T- and P-symmetric systems with spin–orbit coupling. Because of the twofold band degeneracy imposed by \({\mathfrak{T}}={\rm{PT}}\) symmetry, there are two excitation channels where \(\left|n\right\rangle\) at energy −∣E∣ can be excited to +∣E∣ as shown in Fig. 2d. Even when one channel from \(\left|n\right\rangle\) to \({\rm{PC}}\left|n\right\rangle\) is excluded by (PC)2 = −1, there exits another channel from \(\left|n\right\rangle\) to \({\rm{PT}}({\rm{PC}}\left|n\right\rangle )={\rm{TC}}\left|n\right\rangle\). In the absence of spin–orbit coupling, each spin sector forms nondegenerate states so that conditions (i) and (ii) apply [see Fig. 2e and section 3 in the “Methods”].

Crossover to clean-limit optical responses

In real materials, the disorder is present even in very clean samples. We thus need to compare the magnitude of the disorder-mediated and intrinsic responses to characterize the detectability of the latter. We estimate the magnitude of the intrinsic response by counting the dimension of the conductivity tensor in Eq. (3), which gives \({\sigma }_{{\rm{int}}}(\omega )\simeq \frac{{e}^{2}}{h}\frac{1}{\hslash \omega }\frac{{(2{{\Delta }})}^{2}}{{E}_{F}}{k}_{F}^{d-2}{\alpha }^{2}\) above the gap, where kF and EF are the Fermi wave number and Fermi energy. Here, 0 ≤ α ≤ 1 is the ratio between the dominant pairing Δ and the pairing that are responsible for the optical conductivity. Comparing this with the disorder-mediated response, we obtain

$$\frac{{\sigma }_{{\rm{int}}}(\omega )}{{\sigma }_{{\rm{dis}}}(\omega )}\simeq \frac{\omega }{2{{\Delta }}}\frac{l}{{\xi }_{0}}{\left(\frac{2{{\Delta }}}{{E}_{F}}\right)}^{2}{\alpha }^{2}$$
(6)

above the superconducting gap, where we use that σdis(ω) ≈ σn(ω)5,7 as shown in Fig. 1d, where σn is the Drude conductivity in the normal state (see section 5 in the Methods). We thus find that σint ≳ σdis at ℏω ~ 2Δ when

$$x\equiv l/{\xi }_{0}\,> \, {x}_{c}={({k}_{F}{\xi }_{0})}^{2}{\alpha }^{-2}.$$
(7)

Since xc ≫ 1 in general because kFξ0 ~ EF/Δ ≫ 1 (EF/Δ is about 104 for pure metals and 102 for most unconventional superconductors), the clean limit for optical responses is realized in samples much cleaner than that are usually thought to be clean, just satisfying x > 1. This explains how the Mattis–Bardeen-type theories have successfully calculated optical conductivity even in clean superconductors.

Application to FeSe

In FeSe, however, intrinsic optical responses can make up a significant portion of the observed signal in far-infrared optical measurements. FeSe is a clean quasi-two-dimensional material that has a remarkably large ratio Δ/EF ≳ 0.137 with significant spin–orbit coupling comparable to the Fermi energy38 and strongly orbital-dependent pairing14, such that α ~1 is expected. It, therefore, satisfies all the requirements for significant intrinsic optical responses. Here, we use the low-energy model of FeSe in ref. 39 to demonstrate our theory, focusing on the Fermi surface near Γ = (0, 0) for simplicity (see section 6 in the “Methods”).

We consider six constant pairing functions Δ1, Δ2, Δ3, Δ4a, Δ4b, and Δ5 that preserve time-reversal symmetry whose matrix forms and symmetries are given in the “Methods” and Table 2. All of them have even parity, but spin-triplet pairing can occur due to their multi-orbital nature (Δ3 and Δ4a,4b). As we show in the “Methods” section 7, for even parity pairing, optical transitions are not forbidden within each Mz eigenspace in spin–orbit coupled systems. For Δi = 1,2,3,5 pairing, \({\mathfrak{C}}\) symmetry does not exist within a mirror sector. On the other hand, for Δi = 4a,4b, \({\mathfrak{C}}\) symmetry exists but satisfies \({{\mathfrak{C}}}^{2}=1\) in each mirror sector. In accordance with our analysis, all multi-band pairing Δi = 2,3,4a,4b,5 allow for non-zero optical responses [Fig. 3a–d]. In the case of Δi = 4a,4b, and Δ5 pairing, optical conductivity tensors are non-zero down to zero frequency because of their gapless spectrum due to the BFS and Dirac points, respectively [Fig. 3d, e].

Table 2 Properties of time-reversal-symmetric constant pairing functions in a 2D model of FeSe at Γ.
Full size table
Fig. 3: Optical conductivity in a model of superconducting FeSe near Γ at zero temperature.
figure 3

Model parameters in the normal state are adapted from ref. 39 (see section 6 in the Methods). The xx and yy components of the conductivity tensor is shown in red and blue, respectively, in (a–d, f, g). a–d Nonzero intrinsic optical conductivity tensors for each constant pairing function. See Methods and Table 2 for the matrix form and symmetries of the six constant pairing functions Δ1–Δ4a, Δ4b, and Δ5. The case of the Δ1 pairing is not shown as the conductivity is identically zero. e Superconducting gap similar to the experimentally observed gap. FS and the ellipse enclosing it represent the Fermi surface. Red and blue curves correspond to the choice of pairing functions (i) Δ1 = 4.09 meV, Δ2 = 4.82 meV, and Δ3 = 1.93 meV or (ii) Δ1 = 8.98 meV, Δ2 = 9.39 meV, and Δ3 = 0 meV, respectively. They are least-square fits with and without Δ3 to \({{\Delta }}(\theta )=| 2.06+1.42\cos (2\theta )-0.44\cos (4\theta )|\) (shown as a black curve) that was obtained in ref. 15 from experimental data. f, g Conductivity with pairing functions used in (e). σint is the internal optical conductivity in the superconducting state (solid lines), and σn is the Drude conductivity in the normal state (dashed lines). The disorder-mediated conductivity in the superconducting state is expected to be comparable to σn. h Ratio of σint and σn. Red and blue curves are for parameters in (f), and magenta and cyan are for parameters in (g). xx and yy indicate the component of the conductivity tensor. σxy and σyx are not shown in all plots because they vanish due to Mx symmetry.

Full size image

In experiments, a highly anisotropic pairing gap was observed14,15, having a sinusoidal shape with 2–3 meV peak at ky = 0 and almost zero deep at kx = 0. Supposing that the gap function belongs to the trivial representation of the symmetry group, we can obtain a similar anisotropic gap with various combinations of Δ1, Δ2, and Δ3. For example, we obtain (i) Δ1 = 4.09 meV, Δ2 = 4.82 meV, and Δ3 = 1.93 meV and (ii) Δ1 = 8.98 meV, Δ2 = 9.39 meV, and Δ3 = 0 meV, respectively, by least-square-fitting with and without spin–orbit coupled pairing Δ3 to the function \({{\Delta }}(\theta )=| 2.06+1.42\cos (2\theta )-0.44\cos (4\theta )|\) derived from experimental data15. See Fig. 3e. Despite the huge difference in pairing functions for (i) and (ii), the obtained conductivity is quite similar, as shown in Fig. 3f–h. When ℏτ−1 = 0.1 meV (l/ξ0 ~ 102), the xx component of the intrinsic optical conductivity in the superconducting state \({\sigma }_{{\rm{int}}}^{xx}\) exceeds the Drude conductivity in the normal state \({\sigma }_{n}^{xx}\) at around ℏω ~2 meV in both cases. Since the disorder-mediated response in the superconducting state is comparable to the normal-state response σn, the intrinsic response dominates above ℏω ~2 meV along the x direction.

Discussion

Our theory establishes the existence of the true clean-limit optical responses beyond the Mattis–Bardeen theory in multiband superconductors. While we focus on linear responses in the current work, our classification of optical transitions applies to nonlinear optical responses also30. Since nonlinear optical conductivity tensors have more components than the linear counterpart, they give richer information on the symmetry of the system. For instance, it is hard to detect inversion symmetry breaking from linear optical responses. On the other hand, since second-order optical responses are allowed only when inversion symmetry is broken, they directly reveal the presence of inversion symmetry30. As such, various optical measurements can be used in the study of clean multiband superconductors. We anticipate an immediate impact of our work on the optical study of the exotic superconductivity in FeSe. Furthermore, our theory may be relevant to the recently discovered 2D superconductivities reaching Δ/EF > 0.1 in twisted trilayer graphene40,41 and ZrNCl42. As the synthesis of extremely clean superconductors advances further, our results will become relevant to more materials.

Methods

Selection rule by \({\mathfrak{C}}\) symmetry

Equation (5) can be simply derived as follows.

$$\left\langle {\mathfrak{C}}\cdot n{\bf{k}}| {V}^{a}({\bf{k}})n{\bf{k}}\right\rangle =\left\langle {\mathfrak{C}}{V}^{a}({\bf{k}})\cdot n{\bf{k}}| {{\mathfrak{C}}}^{2}\cdot n{\bf{k}}\right\rangle \\ ={{\mathfrak{C}}}^{2}\left\langle {\mathfrak{C}}\cdot n{\bf{k}}| {\mathfrak{C}}{V}^{a}({\bf{k}}){{\mathfrak{C}}}^{-1}| n{\bf{k}}\right\rangle \\ ={\epsilon }_{{\mathfrak{C}},V}{{\mathfrak{C}}}^{2}\left\langle {\mathfrak{C}}\cdot n{\bf{k}}| {V}^{a}({\bf{k}})| n{\bf{k}}\right\rangle .$$
(8)

We use that \({\mathfrak{C}}\) is a anti-unitary operator in the first line, use that \({{\mathfrak{C}}}^{2}=\pm 1\) is a number, and Va is Hermitian in the second line, and define \({\epsilon }_{{\mathfrak{C}},V}=\pm1\) by \({\mathfrak{C}}{V}^{a}({\bf{k}}){{\mathfrak{C}}}^{-1}={\epsilon }_{{\mathfrak{C}},V}{V}^{a}({\bf{k}})\) in the third line. Let us recall that \({V}^{a}({\bf{k}})={\tau }_{z}{\partial }_{{k}_{a}}H({\bf{k}})\) for a BdG Hamiltonian H(k). \({\mathfrak{C}}\) anti-commutes with τz because C anti-commutes with τz while a physical unitary operator Ug that combine to define \({\mathfrak{C}}={U}_{g}C\) commutes with τz. Since the \({\mathfrak{C}}\) symmetry condition imposes \({\mathfrak{C}}H({\bf{k}}){{\mathfrak{C}}}^{-1}=-H({\bf{k}})\), we obtain \({\epsilon }_{{\mathfrak{C}},V}=1\), i.e., \({\mathfrak{C}}\) commutes with Va(k). Equation (5) then follows.

We note that normal-state systems with an emergent \({\mathfrak{C}}\) symmetry follow a different selection rule because they satisfy \({\epsilon }_{{\mathfrak{C}},v}=-1\). The difference comes from the fact that, in the normal state, all quasi-particles are electronic quasi-particles that couple to the gauge field with the equal charge −e (i.e., the emergent \({\mathfrak{C}}\) does not reverse the gauging charge). In this case, the velocity operator is \({v}^{a}({\bf{k}})={\partial }_{{k}_{a}}h({\bf{k}})\), where h is the \({\mathfrak{C}}\)-symmetric normal-state Hamiltonian, so that \({\mathfrak{C}}{v}^{a}({\bf{k}}){{\mathfrak{C}}}^{-1}={\partial }_{{k}_{a}}{\mathfrak{C}}h({\bf{k}}){{\mathfrak{C}}}^{-1}=-{v}^{a}({\bf{k}})\). The selection rule is then \(\left\langle {\mathfrak{C}}\cdot n{\bf{k}}| {v}^{a}({\bf{k}})| n{\bf{k}}\right\rangle =0\) when \({{\mathfrak{C}}}^{2}=+1\), which is opposite to the superconducting case.

Optical excitations with a single-band condition

Let us suppose that the normal state is described by a single band, i.e., h(k) = ξ(k) is a 1 × 1 matrix. We consider the normal state having time-reversal symmetry or inversion symmetry (or other symmetries whose action is equivalent to them43) because only then pairing between two electrons ck and c−k effectively occurs at the Fermi level. Then, ξ(k) = ξ(−k) such that Va(k) = ℏ−1∂aξ(k)τ0 has vanishing inter-band matrix components for all k, where τ0 is the 2 × 2 identity matrix with the particle–hole indices. It follows that superconductivity in a single-band metal cannot exhibit nontrivial optical conductivity in the clean limit. The same is true when the band has spin degeneracy at every k, because h(k) = ξ(k)σ0 is again proportional to the identity matrix such that the velocity operator is diagonal.

These constraints can be understood from the \({\mathfrak{C}}\) symmetry. Let us note that no identity term appears in the BdG Hamiltonian because we consider ξ(k) = ξ(−k). Thus, general two-band BdG Hamiltonian takes the form H = gxτx + gyτy + gzτz. It always satisfies \({\mathfrak{C}}H({\bf{k}}){{\mathfrak{C}}}^{-1}=-H({\bf{k}})\) for \({\mathfrak{C}}=i{\tau }_{y}K\), which satisfies \({{\mathfrak{C}}}^{2}=-1\). This symmetry blocks optical transitions by Eq. (5). Let us consider the case where bands are twofold degenerate due to \({{\mathfrak{T}}}^{2}=-1\) symmetry of the BdG Hamiltonian, where σy is an effective-spin Pauli matrix. As we assume time reversal or inversion symmetry, we also have \({\mathfrak{C}}\) and S symmetries. We consider two cases with \({{\mathfrak{C}}}^{2}=1\) and \({{\mathfrak{C}}}^{2}=-1\) by taking {\({\mathfrak{T}}=i{\sigma }_{y}K\), \({\mathfrak{C}}={\tau }_{x}K\), S = τxσy} and {\({\mathfrak{T}}=i{\tau }_{z}{\sigma }_{y}K\), \({\mathfrak{C}}=i{\tau }_{y}K\), S = τxσy}, respectively, such that H = ξτzσ0 + Δsτyσy, and \(H=\xi {\tau }_{z}{\sigma }_{0}+({{\boldsymbol{\Delta }}}_{t}\cdot {{\boldsymbol{\sigma }}}_{y}i{\sigma }_{y}{\tau }_{+}+h.c.)\), where τ+ = (τx + iτy)/2. Since these correspond to the spin-singlet and spin-triplet pairing, there is a continuous spin rotation symmetry around an axis, such that the spin around the axis is a good quantum number. Each spin sector is thus described by a two-band BdG Hamiltonian having a \({\mathfrak{C}}\) symmetry \({{\mathfrak{C}}}^{2}=-1\).

We can extend the above results to show that the lowest-energy excitations, from − E to E, are forbidden when a multi-band system satisfies the zero superconducting fitness20,44 condition, i.e., [h(k), Δ*(−k)] = 0, which always holds when Δ(k) or h(k) is proportional to the identity matrix. After we take simultaneous eigenstates of h(k) and Δ*(−k), the BdG Hamiltonian decomposes into a set of 2 × 2 blocks (4 × 4 blocks, in the presence of spin degeneracy), each of which correspond to a single-band superconductivity. It follows that transitions between two states with energies −E and E are forbidden.

Let us explain it in more detail. When the zero superconducting fitness condition is satisfied, one can take eigenstates \(\left|\alpha {\bf{k}}\right\rangle\) that satisfy \(h({\bf{k}})\left|\alpha {\bf{k}}\right\rangle ={\xi }_{\alpha }({\bf{k}})\left|\alpha {\bf{k}}\right\rangle\), \({{{\Delta }}}^{* }(-{\bf{k}})\left|\alpha {\bf{k}}\right\rangle ={{{\Delta }}}_{\alpha }^{* }(-{\bf{k}})\left|\alpha {\bf{k}}\right\rangle\). In this basis, the BdG Hamiltonian is diagonalized into blocks labeled by α:

$${H}_{\alpha }({\bf{k}})=\left(\begin{array}{ll}{\xi }_{\alpha }({\bf{k}})&{{{\Delta }}}_{\alpha }({\bf{k}})\\ -{{{\Delta }}}_{\alpha }^{* }(-{\bf{k}})&-{\xi }_{\alpha }(-{\bf{k}}),\end{array}\right)$$
(9)

where the basis states (1 0)T and (0 1)T correspond to \(\left|\alpha {\bf{k}}\right\rangle\) and \({\left|\alpha -{\bf{k}}\right\rangle }^{* }\), respectively. We assume either time reversal symmetry or inversion symmetry of the normal states, such that ξα(k) = ξα(−k).

The energy eigenstates of the BdG Hamiltonian are then

$$\begin{array}{lll}\left|\alpha ,+,{\bf{k}}\right\rangle &=&\left(\begin{array}{l}\cos {\theta }_{\alpha }({\bf{k}})\left|\alpha {\bf{k}}\right\rangle \\ \sin {\theta }_{\alpha }({\bf{k}}){\left|\alpha -{\bf{k}}\right\rangle }^{* }\end{array}\right),\\ \left|\alpha ,-,{\bf{k}}\right\rangle &=&\left(\begin{array}{l}-\sin {\theta }_{\alpha }({\bf{k}})\left|\alpha {\bf{k}}\right\rangle \\ \cos {\theta }_{\alpha }({\bf{k}}){\left|\alpha -{\bf{k}}\right\rangle }^{* }\end{array}\right),\end{array}$$
(10)

where \(\cos {\theta }_{\alpha }({\bf{k}})={\xi }_{\alpha }({\bf{k}})/{E}_{\alpha ,+}({\bf{k}})\), \(\sin {\theta }_{\alpha }({\bf{k}})={{{\Delta }}}_{\alpha }({\bf{k}})/{E}_{\alpha ,+}({\bf{k}})\), and \({E}_{\alpha ,\pm }({\bf{k}})=\pm \sqrt{{\xi }_{\alpha }^{2}({\bf{k}})+{{{\Delta }}}_{\alpha }^{2}({\bf{k}})}\). We have

$$\left\langle \alpha ,+,{\bf{k}}| {v}^{a}| \beta ,-,{\bf{k}}\right\rangle =\cos \theta ({\bf{k}})\sin \theta ({\bf{k}})\left[{v}_{\alpha \beta }^{a}({\bf{k}})+{v}_{\beta \alpha }^{a}(-{\bf{k}})\right]\\ =0\quad (\,\text{for}\,\alpha =\beta ),$$
(11)

where the normal-state velocity operator

$${v}_{\alpha \beta }^{a}({\bf{k}})={\hslash }^{-1}\left\langle \alpha {\bf{k}}| {\partial }_{a}h({\bf{k}})| \beta {\bf{k}}\right\rangle$$
(12)

satisfy \({v}_{\alpha \alpha }^{a}({\bf{k}})=-{v}_{\alpha \alpha }^{a}(-{\bf{k}})\) due to either time-reversal or inversion symmetry. Thus, assuming nondegenerate states, we immediately see that all transitions from Eα,−(k) to Eα,+(k) = −Eα,−(k) are forbidden because of the vanishing velocity matrix elements.

The transition between two states with energies −E and E is forbidden even when spin (or pseudo-spin in spin–orbit coupled systems) degeneracy exists. Spin degeneracy can be imposed either by PT symmetry or spin rotation symmetry. In either case, the normal-state velocity matrix element between two spin-degenerate states \(\left|\alpha {\bf{k}}\right\rangle\) and \(\left|\beta {\bf{k}}\right\rangle\) vanishes by symmetry, i.e., \({v}_{\alpha \beta }^{a}({\bf{k}})=0\). Let us first consider the case where the degeneracy is due to PT symmetry satisfying (PT)2 = − 1. Then we have \({v}_{\alpha ,\beta \,=\, {\rm{PT}}\alpha }^{a}({\bf{k}})=-{v}_{\alpha ,{\rm{PT}}\alpha }^{a}({\bf{k}}).\) When spin rotation symmetry exists, the constraint

$${v}_{\uparrow \downarrow }^{a}({\bf{k}})=0$$
(13)

simply follows from the spin-rotation invariance of the velocity operator. Equation (11) then shows \(\left\langle \alpha ,+,{\bf{k}}| {v}^{a}| \beta ,-,{\bf{k}}\right\rangle =0\) for α and β related by either PT or a spin rotation. Combining this with \(\left\langle \alpha ,+,{\bf{k}}| {v}^{a}| \alpha ,-,{\bf{k}}\right\rangle =0\), we see that all transition channels from Eα,− = Eβ,− to −Eα,− = −Eβ,− are forbidden.

Optical excitations with spin rotation symmetries

Let us consider the EAZ classes within spin sectors for example. We first assume that time reversal, spatial inversion, a spin U(1) rotation and a spin π rotation (around an axis perpendicular to the U(1) axis) symmetries are all present. Let us take energy eigenstates such that they carry a definite spin along the z direction. In the case of triplet pairing, this means that we take the z-direction as the triplet spin direction because continuous spin rotation symmetries around other directions are broken. Within each spin sector, bands are nondegenerate at generic momenta. As we show in Fig. 2d, PC flips the spin because of particles-hole conjugation. However, the combination of spin π rotation, which is −iσy for singlet pairing (−iτzσy for triplet pairing), and PC acts within a spin sector, so symmetry under this \({\mathfrak{C}}\)-type operation constraints optical excitations through Eq. (5). Optical excitations between \({\mathfrak{C}}\)-related pairs are allowed when \({{\mathfrak{C}}}^{2}={(-i{\sigma }_{y}{\rm{PC}})}^{2}=-{(\rm{PC})}^{2}=1\) and \({{\mathfrak{C}}}^{2}={(-i{\tau }_{z}{\sigma }_{y}{\rm{PC}})}^{2}={(\rm{PC})}^{2}=1\), respectively, for singlet and triplet pairing, where BFSs are stable. Since (PC)2 = +1 (−1) for even- and odd-parity pairing (see section 4 in the “Methods” below), this requires exotic odd-parity singlet pairing or even-parity triplet pairing, which is possible only when multi-orbital pairing, e.g., an orbital triplet, is realized. Alternatively, if inversion symmetry is broken or spin–orbit coupling is not negligible, optical excitations are allowed with fully gapped superconductivity. Let us note that sz-preserving spin-orbit coupling is enough to allow optical excitations, because it breaks spin rotation symmetries around other axes such that \({\mathfrak{C}}\propto i{\sigma }_{y}PC\) is broken in each spin sector. In general, spin–orbit coupling breaks all spin rotation symmetries, so two excitation channels are allowed [Fig. 2e].

Symmetry operator and pairing symmetry

Let ug be a unitary operator that acts on space as g:k → gk. Suppose that it is a symmetry operator of the normal state, i.e.

$${u}_{g}h({\bf{k}}){u}_{g}^{-1}=h(g{\bf{k}}),$$
(14)

and the pairing function has eigenvalues \({e}^{i{\theta }_{g}}\) under Ug, i.e.

$${u}_{g}{{\Delta }}({\bf{k}}){u}_{g}^{T}={e}^{i{\theta }_{g}}{{\Delta }}(g{\bf{k}}).$$
(15)

Due to the non-trivial symmetry transformation of the pairing function, the BdG Hamiltonian is symmetric under

$${U}_{g}=\left(\begin{array}{ll}{u}_{g}&0\\ 0&{e}^{i{\theta }_{g}}{u}_{g}^{* }\end{array}\right),$$
(16)

which rotates the hole sector by \({e}^{i{\theta }_{g}}\) more

$${U}_{g}H({\bf{k}}){U}_{g}^{-1}=H(g{\bf{k}}).$$
(17)

Ug satisfies the following commutation relation with the particle–hole conjugation operator C

$${U}_{g}C={e}^{i{\theta }_{g}}C{U}_{g}.$$
(18)

Let us take two examples.

  1. 1.

    Ug = P is spatial inversion: \({e}^{i{\theta }_{g}}=+1\) and −1 indicates even-parity and odd-parity pairing. Thus, PC = +CP (PC = −CP) for even-parity (odd-parity) pairing.

  2. 2.

    Ug is a spin rotation around the y-axis by π: \({e}^{i{\theta }_{g}}=+1\) always for a spin-singlet pairing, and \({e}^{i{\theta }_{g}}=+1\) (−1) when the pairing function is a spin-triplet with its spin parallel (perpendicular) to the y-axis.

Estimates of disorder-mediated and intrinsic responses in the clean regime

The disorder-mediated response in the superconducting state is comparable to the Drude response in the normal state. When the light frequency ω is much larger than the inverse relaxation time Γ, which is the case in the clean regime ℏΓ ≪ Δ ≲ ℏω, the Drude conductivity is \({\sigma }_{n}(\omega )={\sigma }_{0}\frac{{{{\Gamma }}}^{2}}{{\omega }^{2}\,+\,{{{\Gamma }}}^{2}}\approx {\sigma }_{0}\frac{{{{\Gamma }}}^{2}}{{\omega }^{2}}\). Here, \({\sigma }_{0}=\frac{n{e}^{2}\tau }{m}\approx \frac{{e}^{2}}{h}{k}_{F}^{d-2}({\hslash }^{-1}{E}_{F}\tau )\) is the DC conductivity, where \(m={\hslash }^{2}{k}_{F}^{2}/(2{E}_{F})\), and \(n \sim {k}_{F}^{d}\). Since σdis ~σn, we have

$${\sigma }_{{\rm{dis}}}(\omega ) \sim \frac{{e}^{2}}{h}{k}_{F}^{d-2}\frac{{E}_{F}}{2{{\Delta }}}\left(\frac{{\xi }_{0}}{l}\right){\left(\frac{2{{\Delta }}}{\hslash \omega }\right)}^{2},$$
(19)

where we use EF ~ ℏvFkF, \({{\Delta }} \sim \hslash {v}_{F}{\xi }_{0}^{-1}\), and l = vFΓ−1.

To estimate the intrinsic response. let us note that the inter-band velocity operator is linear in the leading order of \({{\Delta }}^{\prime} /{E}_{F}\), where \({{\Delta }}^{\prime}\) is the largest multiband pairing that is allowed to generate interband transitions by the selection rules in Eqs. (4) and (5). The conductivity tensor in Eq. (3) is thus proportional to \({{\Delta }}{^{\prime} }^{2}\). Using that the delta function has the dimension ω−1, and using the Fermi energy and wave number as other scales, we obtain

$${\sigma }_{{\rm{int}}}(\omega ) \sim \frac{{e}^{2}}{h}\frac{1}{\hslash \omega }\frac{{(2{{\Delta }})}^{2}}{{E}_{F}}{k}_{F}^{d-2}{\alpha }^{2},$$
(20)

where \(\alpha =\left|{{\Delta }}^{\prime} /{{\Delta }}\right|\) is equal to or smaller than one since Δ is the dominant pairing strength by definition. It follows that

$$\frac{{\sigma }_{{\rm{int}}}(\omega )}{{\sigma }_{{\rm{dis}}}(\omega )} \sim \frac{\omega }{2{{\Delta }}}\frac{l}{{\xi }_{0}}{\left(\frac{2{{\Delta }}}{{E}_{F}}\right)}^{2}{\alpha }^{2}.$$
(21)

above the superconducting gap frequency.

FeSe model at the Γ point

If we regard FeSe as a 2D system, it has three Fermi surfaces around Γ = (0, 0), X = (π, 0), and Y = (0, π), respectively, in the 1-Fe Brillouin zone16. Here, we consider the Fermi surface near Γ. The Fermi surface of FeSe at Γ consists mainly of two orbital degrees of freedom dyz and dxz, so we take \(\psi ({\bf{k}})={\left({d}_{yz}({\bf{k}}),{d}_{xz}({\bf{k}})\right)}^{T}\) as the basis state. At zero temperature, the normal-state Hamiltonian has the form

$${h}_{{{\Gamma }}}={h}_{0}+{h}_{{\rm{nem}}}+{h}_{{\rm{SOC}}},$$
(22)

where \({h}_{0}={\epsilon }_{{{\Gamma }}}-A({k}_{x}^{2}+{k}_{y}^{2})+B({k}_{x}^{2}-{k}_{y}^{2}){\rho }_{z}-2C{k}_{x}{k}_{y}{\rho }_{x}\) is the most general spinless Hamiltonian in the tetragonal phase up to second order in k, hnem = −Dρz is the constant part of the nematic terms that develops below Tnem ~90 K45,46,47, and hSOC = λρyσz is the constant spin–orbit coupling. Here, ρi=x,y,z and σi=x,y,z are the Pauli matrices for orbital and spin degrees of freedom, respectively. h0 and hSOC have tetragonal D4h symmetry under mirror mx = −iρzσx, my = iρzσy, and mz = −iσz, and fourfold rotation \({c}_{4z}=i{\rho }_{y}{e}^{-i\frac{\pi }{4}{\sigma }_{z}}\), while hnem breaks c4z symmetry down to c2z symmetry. All have time-reversal t = iσyK symmetry. In our numerical calculations, we take parameters used in ref. 39, which are ϵΓ = −9 meV, A = 700 meVÅ2, B = C = 484 meVÅ2, D = 15 meV, and λ = 10 meV.

The bulk FeSe shows superconductivity below 8 K without a sign of time-reversal symmetry breaking16,45. We consider constant pairing functions invariant under time reversal. Since there are six 4 × 4 matrices that are invariant under t = iσyK, the pairing function has the form

$${{\Delta }}({\bf{k}}) =\left({{{\Delta }}}_{1}+{{{\Delta }}}_{2}{\rho }_{z}+{{{\Delta }}}_{3}{\rho }_{y}{\sigma }_{z}\right.\\ \left.\quad+\;{{{\Delta }}}_{4a}{\rho }_{y}{\sigma }_{x}+{{{\Delta }}}_{4b}{\rho }_{y}{\sigma }_{y}+{{{\Delta }}}_{5}{\rho }_{x}\right)i{\sigma }_{y},$$
(23)

where Δ1, Δ2, Δ3, Δ4a, Δ4b, and Δ5 are all independent of k. The pairing symmetry of each term under the D4h point group is shown in Table 2. Let us note that all are even-parity pairing, i.e., invariant under p = mxmymz = 1. The orbital-singlet nature of Δ3 and Δ4a,4b allows them to be spin triplet even though they have even parity.

Optical excitations in M z, P, and T-symmetric 2D superconductors

In two-dimensional systems perpendicular to the z axis, Mz symmetry divides eigenstates into two distinct eigenspaces with Mz eigenvalues λ = ±i. It imposes a selection rule.

Let MzC = ηCMz, where η = ±1, and \({M}_{z}\left|n\right\rangle ={\lambda }_{n}\left|n\right\rangle\). Then, \({M}_{z}{\rm{PC}}\left|n\right\rangle =\eta {\rm{PC}}{M}_{z}\left|n\right\rangle =\eta {\rm{PC}}{\lambda }_{n}\left|n\right\rangle =\eta {\lambda }_{n}^{* }{\rm{PC}}\left|n\right\rangle =-\eta {\lambda }_{n}{\rm{PC}}\left|n\right\rangle\). Thus, \(\left|n\right\rangle\) and \({\rm{PC}}\left|n\right\rangle\) has the different eigenvalues when η = 1 and has same eigenvalues when η = −1. In the former case, the optical transition from \(\left|n\right\rangle\) to \({\rm{PC}}\left|n\right\rangle\) is forbidden by Mz symmetry because they are indifferent eigenspaces (and the velocity operator does not change the eigenspace), but the transition from \(\left|n\right\rangle\) to \(S\left|n\right\rangle\) is allowed. On the other hand, in the latter case, the transition from \(\left|n\right\rangle\) to \({\rm{PC}}\left|n\right\rangle\) within the same mirror sector is forbidden when the pairing is odd-parity such that (PC)2 = −1.

In summary, optical transitions between particle-hole- and chiral-related states are forbidden in P and T-symmetric systems when MzC = −CMz and PC = −CP. Similar constraints can appear in one-dimensional systems also due to mirror symmetries.