#  Odd-frequency pairing inherent in Bogoliubov Fermi liquid

Tatsuya Miki, Shun-Ta Tamura, Shoma Iimura, and Shintaro Hoshino  Department
of Physics, Saitama University, Shimo-Okubo, Saitama 338-8570, Japan

###### Abstract

The disorder and interaction effects on Bogoliubov-Fermi surfaces with
preserved inversion symmetry are studied for a low-energy effective model
coupled to bosonic degrees of freedom. It is shown that the non-ideal
Bogoliubov quasiparticles (bogolons) generically induce the odd-frequency pair
amplitude which reflects a Cooper pairing at different time. The self-energy
of bogolons is mainly contributed by the disorder effects in the low frequency
limit as in the usual electron liquid. Depending on the choice of the
parameters, there are two kinds of solutions: one is frequency-independent
(but with sign function of frequency) and the other is proportional to the
inverse of the frequency, which exist in both the normal and anomalous parts
of the self-energy. These characteristic self-energy structures are clearly
reflected in the single-particle spectrum. Since the bogolons are originally
composed of electrons, the connection between the two is also sought using the
concrete $j=3/2$ fermion model, which reveals that the odd-frequency pairing
of bogolons is mainly made of the electrons’ odd-frequency pairing.

The superconductivity is induced by the interactions among electrons near the
Fermi surfaces. In the resultant ground state, the Cooper pair condensation
energy is gained by energy gap formation near the Fermi level. While the Fermi
surface usually disappears in the pairing state, they can remain in some
superconducting states [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], where the
elementary excitations near the Fermi surfaces are composed of Bogoliubov
quasiparticles (bogolons) [13]. For the time-reversal symmetry broken system
with preserved inversion symmetry, such Bogoliubov-Fermi surfaces are stable
as they are topologically protected [5]. Whereas the bogolon is not a simple
charged particle, it can carry energy. Hence the thermal properties such as
specific heat and thermal conductivity are expected to be similar to the
conventional Fermi liquid of electrons [14, 15] and are potentially observed
in the actual materials [16, 17, 18]. However, since the bogolons are
quasiparticles in the superconducting state, the physical properties should be
different from those of the electrons. Therefore it is desirable to clarify
the difference between the Fermi liquid and the Bogoliubov Fermi liquid, the
latter of which is realized for the non-ideal bogolons generically. Recently,
the effect of interactions are considered for bogolons and the possible
instabilities are studied [19, 20, 21, 22]. Here, we show that the impurity
and correlation effects on this Bogoliubov Fermi surface generate purely odd-
frequency pairing amplitude at low energies, which gives a clear distinction
from the normal Fermi liquid state of electrons.

The odd-frequency pairing is proposed as a possible exotic ground state of the
fermionic systems where the Cooper pair amplitude has an odd function in
relative time or frequency, meaning that the pair formation occurs only at
different time [23, 24, 25, 26, 27, 28, 29, 30, 31]. The actual realizations
as a bulk state have been proposed in correlated electron models [32, 33, 34,
35, 36, 37, 38, 39, 40] and Kondo lattices [41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52]. However, it has been argued that the spatially uniform odd-
frequency pairing cannot be a thermodynamically stable state in a conventional
framework [53]. Then, several possibilities to remedy this problem are
proposed such as spatially modulated superconducting states [53, 43, 48] and
non-Hermitian description [54, 55, 56]. As for the latter scenario, potential
problems have also been pointed out [57]. In any case, it has been recognized
that the spatially uniform and purely odd-frequency pairing state without
conventional (even-frequency) pairing cannot simply be realized. On the other
hand, such exotic pairing has been discussed also in edge or interface [58,
59, 60, 61], which are secondarily induced from the even-frequency pairing in
bulk.

Figure 1:  Schematic phase diagram for the Bogoliubov Fermi liquid. The
horizontal axis with $U$ indicates the deviation from the non-interacting
bogolons.

In this paper, we demonstrate the emergence of the spatially uniform odd-
frequency pair in the world not of electrons but of bogolons. This odd-
frequency pairing is induced by a non-ideality of bogolons, where the broken
gauge-symmetry plays a crucial role as different from the electronic
interactions. As shown in the following, the induced odd-frequency pair is
closely related to the energy structure of the usual Fermi-liquid self-energy
existing in the diagonal Green function for electrons, but for bogolons it is
reflected also in the off-diagonal part, i.e., pair amplitude and pair
potential.

In the previous works [19, 20, 21, 22], the possibility of the spontaneous
symmetry breaking induced by interactions among bogolons is considered. Here,
we concentrate on the superconducting state above such low-temperature phase
($T>T_{b}$) but below the superconducting transition temperature ($T<T_{c}$).
Namely, we deal with a “normal state” of bogolons. Our situation is sketched
in Fig. 1.

Below, we connect the two important concepts of superconductivity, i.e., the
Bogoliubov Fermi surface and odd-frequency pairing, based on the concrete
model. As in the theories of superconductivity and of Fermi liquid, the
degrees of freedom near the Fermi level are important for the low-temperature
and low-energy properties. Hence we also assume that the dominant contribution
enters through the degrees of freedom near the Bogoliubov Fermi surface. The
non-interacting part of the Hamiltonian for bogolons near the Fermi level is
written as

$\displaystyle\mathscr{H}_{0}$
$\displaystyle=\sum_{\bm{k}}\xi_{\bm{k}}\alpha_{\bm{k}}^{\dagger}\alpha_{\bm{k}}.$
(1)

In the following, we consider the inversion symmetric systems
($\xi_{-\bm{k}}=\xi_{\bm{k}}$) with a time-reversal symmetry breaking where
the Fermi surface is topologically protected [5]. In order to have an
intuition for the energy scales, we write the energy dispersion as
$\xi_{\bm{k}}=\frac{\bm{k}^{2}}{2m_{b}}-\varepsilon_{{\rm F}b}$ for
simplicity, where $m_{b}$ and $\varepsilon_{{\rm F}b}$ are effective mass and
Fermi energy for bogolons ($\hbar=k_{\rm B}=1$). As inferred from the original
electronic system [5] on which bogolons are based, each quantity is roughly
expressed as $m_{b}\sim\frac{\Delta}{\varepsilon_{{\rm F}e}}m_{e}$ and
$\varepsilon_{{\rm F}b}\sim\Delta$ where $m_{e}$ and $\varepsilon_{{\rm F}e}$
are mass and Fermi energy for electrons, and $\Delta$ is the energy scale for
the superconducting gap function. Here we have assumed the magnitude relation
$\varepsilon_{{\rm F}e}>\Delta$. The Fermi wavenumber is given by $k_{{\rm
F}b}=\sqrt{2m_{b}\varepsilon_{{\rm F}b}}$. Thus the Fermi velocity for
bogolons is similar to that of electrons: $v_{{\rm F}b}=\frac{k_{{\rm
F}b}}{m_{b}}\sim v_{{\rm F}e}$. Note that the potential term proportional to
$\alpha^{\dagger}_{\bm{k}}\alpha^{\dagger}_{-\bm{k}}$ cannot exist together
with the inversion symmetry. Such pair potential may be generated with
spontaneous symmetry breaking below the ordering temperature $T_{b}$ for
bogolons as shown in Fig. 1, although we work in the temperature regime
$T>T_{b}$ in the following.

For the superconducting state considered here, the time-reversal and gauge
symmetries are already broken. Nevertheless, in terms of bogolons as in Eq.
(1), the effects from these broken symmetries are not described apparently.
This paradox at first sight is rationalized by considering the self-energy
terms in which the symmetry of the system is reasonably reflected. Then, we
arrive at the important conclusion that the anomalous part can be finite in
the presence of correlations, and its symmetry must be even-parity pairing.
Since the bogolons have no internal degrees of freedom at the Fermi level, and
in order to satisfy the Pauli principle, the pairing state of bogolons must
have odd relative time dependence. Hence the odd-frequency pairing is
reasonably realized.

We proceed along with the conventional microscopic approach for the Fermi
liquid based on the weak-coupling limit [62]. The total Hamiltonian is written
as $\mathscr{H}=\mathscr{H}_{0}+\mathscr{H}_{\rm imp}+\mathscr{H}_{\rm int}$,
where the impurity scattering part is given by

$\displaystyle\mathscr{H}_{\rm imp}$
$\displaystyle=\frac{1}{V}\sum_{\bm{k},\bm{q}}\Big{[}\rho_{\bm{q}}u_{1}(\bm{k},\bm{q})\alpha^{\dagger}_{\bm{k}+\bm{q}}\alpha_{\bm{k}}$
$\displaystyle\hskip
42.67912pt+\rho_{\bm{q}}u_{2}(\bm{k},\bm{q})\alpha^{\dagger}_{\bm{k}+\bm{q}}\alpha^{\dagger}_{-\bm{k}}+{\rm
H.c.}\Big{]},$ (2)

where $\rho_{\bm{q}}=\sum_{i}\mathrm{e}^{\mathrm{i}\bm{q}\cdot\bm{R}_{i}}$ is
the structure factor for the impurity configuration $\\{\bm{R}_{i}\\}$, and
$V$ is a system volume. The second term with
$\alpha^{\dagger}\alpha^{\dagger}$ is characteristic of the bogolon systems.
We also consider the correlation effects in $\mathscr{H}_{\rm int}$. Among the
various interactions, we take the simple model where the bogolons are coupled
with bosons. The interaction term is given by the replacement

$\displaystyle
u_{j}(\bm{k},\bm{q})\rho_{\bm{q}}\to\mathrm{i}g_{j}(\bm{k},\bm{q})\sqrt{\frac{\omega_{0,\bm{q}}}{2}}(b_{\bm{q}}-b_{-\bm{q}}^{\dagger})$
(3)

in Eq. (2), where $g_{j=1,2}$ is the coupling constant, $\omega_{0,\bm{q}}$
the bare boson dispersion, and $b$ ($b^{\dagger}$) the annihilation (creation)
operator of boson. Note that, in the above model, we have assumed that the
Bogoliubov Fermi surfaces are stable as separated from the other bands located
at higher energies. Then, the disorders and interactions dominantly affect the
bogolons located near the Bogoliubov Fermi surfaces. Otherwise, the original
assumption of the stable Bogoliubov Fermi surfaces must be reconsidered.

Figure 2:  Diagrams for (a) normal and (b) anomalous self-energies which are
relevant to the Bogoliubov Fermi liquid.

In order to show the presence of the pure odd-frequency pairing, we use the
simple weak-coupling perturbation theory, which is a minimal description for
the Fermi liquid and is sufficient for generating anomalous Green functions.
The self-energy is contributed by the diagrams shown in Fig. 2, and their
calculation is parallel to the usual Fermi liquid theory [62] and is shown in
Supplementary Material (SM) A and B [63]. The resultant normal part of the
self-energy is given by

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\mathrm{i}\Gamma_{1\bm{k}}{\rm
sgn\,}\varepsilon_{n}+a_{\bm{k}}\mathrm{i}\varepsilon_{n}+\mathrm{i}b_{\bm{k}}(\pi^{2}T^{2}-\varepsilon_{n}^{2}){\rm
sgn\,}\varepsilon_{n},$ (4)

where $\varepsilon_{n}=(2n+1)\pi T$ is the fermionic Matsubara frequency. The
first term corresponds to a quasiparticle damping due to impurity scattering.
The second term represents a renormalization factor, and the third term is
responsible for the damping in the usual Fermi liquid theory. For bogolons,
the anomalous self-energy is also present, whose diagrammatic contribution is
very similar to the normal self-energy as shown in Fig. 2. The anomalous part
is obtained as

$\displaystyle S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\mathrm{i}\Gamma_{2\bm{k}}{\rm
sgn\,}\varepsilon_{n}+c_{\bm{k}}\mathrm{i}\varepsilon_{n}+\mathrm{i}d_{\bm{k}}(\pi^{2}T^{2}-\varepsilon_{n}^{2}){\rm
sgn\,}\varepsilon_{n},$ (5)

which is the same frequency dependence as the usual Fermi liquid. Namely, the
Fermi-liquid self-energy is originally odd in frequency and therefore it
matches well with the requirement of the bogolon anomalous self-energy under
the inversion symmetry. This is the reason why the odd-frequency pair
potential naturally appears for bogolons. We emphasize that the spatially
uniform and purely odd-frequency pairing is realized in the present setup,
which is necessarily accompanied by the normal self-energies.

From the Hermiticity relation, it can be shown that $a$, $b$ and $\Gamma_{1}$
are real, while $c$, $d$ and $\Gamma_{2}$ can be complex. As discussed in Ref.
[5], the Bogoliubov Fermi surfaces can be realized for a chiral $d$-wave
superconductivity with the gap function $\Delta_{\bm{k}}\sim
k_{z}(k_{x}+\mathrm{i}k_{y})$ of the original electrons. In this case, the
anomalous part $c$ (and also $d$) may include the contribution $\propto
k_{z}(k_{x}+\mathrm{i}k_{y})$ from the symmetry argument. Hence the broken
gauge symmetry and time-reversal symmetry are clearly reflected in the
anomalous self-energy of bogolons which is odd in frequency.

The effect of the presence of the odd-frequency pairing is best visualized in
the single-particle spectral functions. To see this, we consider the Green
functions given by

$\displaystyle\begin{pmatrix}G_{\bm{k}}&F_{\bm{k}}\\\
F^{\dagger}_{\bm{k}}&\bar{G}_{\bm{k}}\end{pmatrix}^{-1}$
$\displaystyle=\begin{pmatrix}\mathrm{i}\varepsilon_{n}-\xi_{\bm{k}}&\\\
&\mathrm{i}\varepsilon_{n}+\xi_{\bm{k}}\end{pmatrix}-\begin{pmatrix}\Sigma_{\bm{k}}&S_{\bm{k}}\\\
S^{\dagger}_{\bm{k}}&\bar{\Sigma}_{\bm{k}}\end{pmatrix},$ (6)

where the frequency dependence is omitted. The self-energies satisfy the
relations
$\bar{\Sigma}_{\bm{k}}(\mathrm{i}\varepsilon_{n})=-\Sigma_{\bm{k}}^{*}(\mathrm{i}\varepsilon_{n})$
and
$S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})=-S_{\bm{k}}^{*}(\mathrm{i}\varepsilon_{n})$
as derived from the Hermiticity and inversion symmetry. The dominant
contribution of the self-energies at low energies enters from the impurity
effect. Hence here we focus on the impurity self energies $\Gamma_{1,2}$ in
Eqs. (4) and (5). We note that $\Gamma_{1}>|\Gamma_{2}|$ is required for the
physical behavior, i.e., the positive weight of the spectrum, and it can
indeed be checked in the weak coupling limit [63]. In the following, we
neglect the $\bm{k}$-dependence in the self-energy for simplicity,
corresponding to the spatially local self-energies.

Figure 3:  (a) Density of states, (b) momentum distribution function, (c) pair
potential for the solution of the first kind. (d,e,f) are similar to (a,b,c),
respectively, but for the solution of the second kind. $D_{0}$ is a bare
density of states. The pair amplitude in (c) is normalized by $\Gamma_{1}$ and
in (d) by $|V|$ where the phases are chosen as zero.

The single particle spectrum $A_{\bm{k}}(\varepsilon)$ is obtained from the
imaginary part of the retarded Green function. The density of states is then
evaluated by performing the $\bm{k}$-integral as
$D(\varepsilon)=\sum_{\bm{k}}A_{\bm{k}}(\varepsilon)$, which is shown in Fig.
3(a). The larger $\Gamma_{2}$ makes the higher peak at the Fermi level.
Clearly, the presence of odd-frequency pair is responsible for this
characteristic feature near the Fermi level. Figure 3(b) shows the momentum
distribution function
$n_{\bm{k}}=\langle\alpha_{\bm{k}}^{\dagger}\alpha_{\bm{k}}\rangle$. The sharp
drop at the Fermi energy in the ideal limit is smeared by the damping
$\Gamma_{1}$, but is recovered with increasing $|\Gamma_{2}|$. The inset of
(b) shows the derivative of this function, where the change becomes more
abrupt for the larger pair potential $\Gamma_{2}$ and diverges when
$|\Gamma_{2}|\to\Gamma_{1}$.

We also show the spatially local anomalous pair amplitude
$F=\sum_{\bm{k}}F_{\bm{k}}$. Performing the Fourier transformation, we show in
Fig. 3(c) the imaginary-time $\tau$ dependence of the pair amplitude where we
have taken the zero temperature limit. The functional form clearly shows the
odd-frequency pairing. The value of the pair amplitude is largest at short
time and is discontinuous at $\tau=0$. The asymptotic behavior at long time is
$F^{\dagger}(\tau\to\infty)\sim\tau^{-1}$.

The information of quasiparticles are seen in the real-time evolution of the
retarded Green function $G^{\rm
R}_{\bm{k}}(\varepsilon)=G_{\bm{k}}(\varepsilon+\mathrm{i}0^{+})$ [63]. We can
explicitly evaluate it in the real-time domain and obtain $G_{\bm{k}}^{\rm
R}(t)\sim\mathrm{e}^{\mathrm{i}z_{\bm{k}}t}$. For the high-energy region with
$\xi_{\bm{k}}>|\Gamma_{2}|$, the complex energy is given by
$z_{\bm{k}}=\pm\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}+\mathrm{i}\Gamma_{1}$.
This is a standard form composed of the oscillating part with the
quasiparticle energy $\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}$ and the
damping $\Gamma_{1}$. On the contrary, for $\xi_{\bm{k}}<|\Gamma_{2}|$, we
have
$z_{\bm{k}}=\mathrm{i}(\pm\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}+\Gamma_{1})$
which is pure imaginary. Hence the low-energy part has no oscillating part and
has only damping with two relaxation rates.

For the impurity effect, we have also analyzed the equations with the self-
consistent treatment [63]. We have found that the frequency-independent
$\Gamma_{1,2}$ solution discussed above is not much modified at low energies.
This type is called the first-kind. On the other hand, owing to the non-
linearity of the equations, we have also found another type of solution
depending on the choice of the parameters. This solution of the second-kind
has the frequency dependence

$\displaystyle\Sigma(\mathrm{i}\varepsilon_{n})=\frac{|V|^{2}}{\mathrm{i}\varepsilon_{n}},\
\ \
S^{\dagger}(\mathrm{i}\varepsilon_{n})=\frac{V^{2}}{\mathrm{i}\varepsilon_{n}}$
(7)

at low energies, where $V$ is a complex constant. In this case, we have the
three single-particle excitation energies
$E_{\bm{k}}=0,\pm\sqrt{\xi_{\bm{k}}^{2}+2|V|^{2}}$ from the pole of Green
functions. This energy structure is due to the fact that the same absolute
value of $|V|$ is shared for both normal and anomalous part of self-energy.
The presence of the gap structure and the zero-energy peak at the Fermi level
is the characteristic feature for the density of states as shown in Fig. 3(d),
where the weight of the zero-energy peak is proportional to $|V|$. Figure
3(e,f) show the momentum distribution function and imaginary-time dependence
of local pair amplitude, respectively, for the case with Eq. (7). The pair
amplitude has odd function form but is now smooth at equal imaginary-time. The
asymptotic behaviors are identified as
$F^{\dagger}(\tau)\sim-|V|^{2}\tau\ln(\tau|V|)$ for $\tau\to 0$ and
$F^{\dagger}(\tau)\sim|V|$ for $\tau\to\infty$ [63].

The frequency dependence proportional to the inverse of $\varepsilon_{n}$ for
both the normal and anomalous parts are the same features as the pairing
states in the multichannel Kondo lattices where the conduction electrons
hybridized virtually with the localized fermions [43, 64, 65, 49, 66]. This
suggests that the $1/\varepsilon_{n}$ form is a ubiquitous form of the self-
energy for odd-frequency pairing in bulk, as it is found in the very different
two physical systems.

Thus, the Bogoliubov Fermi liquid shows characteristic single-particle
excitations at low energies, which can be much different from the usual
electron liquids. The above two different behaviors, the first- and second-
kind solutions, are observed depending on the choice of parameters [63], and
hence are dependent on the specific materials. For example, in Fe(Se,S), the
Fermi-liquid-like behavior is observed below the superconducting transition
temperature as probed by thermal measurement [16, 18]. The solution of the
second kind shows the presence of the localized level at the Fermi level,
which should not contribute to the transport phenomena. Hence, the solution of
the first-kind with the form (5) is likely realized in Fe(Se,S). While the
observed behaviors are similar to the Fermi liquid, which implies weak
disorder effects, the tuning of the system may make it clearer to detect the
fingerprints of Bogoliubov-Fermi liquids.

Figure 4:  Wavevector dependence of (a) $C_{M}(\bm{k})$ (multipoles), (b)
$C_{P}(\bm{k})$, and (c) $C_{P^{\dagger}}(\bm{k})$ (multiplet electron and
hole pairs) along the Fermi surface in the $j=3/2$ fermion model. The inset of
(a) shows the two Fermi surfaces near the $k_{x}$ axis located in the
$k_{y}=0$ plane. The horizontal axis of (a–c) is a path shown by arrows in the
inset of (a).

Finally we discuss the connection of the bogolon degrees of freedom to
original electrons. As a simple realization of the Bogoliubov Fermi surface,
we take $j=3/2$ fermion model proposed in Ref. [5], where the symmetric spin-
orbit interaction and time-reversal-symmetry-broken $d+\mathrm{i}d$ pairing
are considered. The pair amplitude of bogolons is connected to the order
parameters of the original electrons. At each $\bm{k}$ point on the Fermi
surface, we can define the electronic multipoles $M^{r}_{\bm{k}}$ and
multiplet pair amplitudes $P^{r}_{\bm{k}},P_{\bm{k}}^{\dagger,r}$ (see SM C
[63] for the concrete forms), which are classified by the rank $r$ defined up
to $2j$. For $M$ ($P,P^{\dagger}$), the rank $r=0,1,2,3$ respectively
correspond to monopole (singlet pair), dipole (triplet pair), quadrupole
(quintet pair), octupole (septet pair) [20]. The order parameters induced by
the odd-frequency pair amplitude of bogolons are given by

$\displaystyle X^{r}_{\bm{k}}(\tau)$
$\displaystyle=C^{r}_{X}(\bm{k})F_{\bm{k}}(\tau)$ (8)

for $X=M,P,P^{\dagger}$. The quantity $C$ is regarded as a kind of
‘susceptibility’, showing how much of multipoles and multiplet pair are
induced from the odd-frequency pair amplitude. For the multiplet pairs, only
the spin-triplet and spin septet pair can be finite due to the odd-function in
time in systems with inversion symmetry.

Figure 4 shows the value of $C_{X}$ along the Fermi surface. For exemplary
demonstration, we take the parameters as the symmetric spin-orbit coupling
$\beta/\varepsilon_{\rm Fe}=0.3$, $s$\- and $d$-wave pair potentials
$\Delta_{0}/\varepsilon_{\rm Fe}=\Delta_{1}/\varepsilon_{\rm Fe}=0.1$ in Refs.
[5, 20] (see also SM C [63]). The inset of Fig. 4(a) shows the Fermi surfaces
in the $k_{y}=0$ plane. We have two Fermi surfaces near the $xy$-plane, and it
has a donut-like shape because of the rotational symmetry around $z$-axis. The
figures (a–c) show that the dominantly induced components are odd-frequency
(spin-triplet and spin-septet) electron/hole pair amplitude ($P,P^{\dagger}$).
Namely, the bogolon odd-frequency pair is mainly composed of the electron odd-
frequency pair. We note that the diagonal multipole components [$M$ in (a)]
also mix with a similar-order of magnitudes.

To summarize, we have demonstrated that the non-ideal bogolons generate purely
odd-frequency pair potential and pair amplitude near the Fermi surface, as
different from the usual Fermi liquid of electrons. The odd-frequency pair
recovers the original symmetry of the electron system in terms of bogolons,
where the broken gauge and time-reversal symmetries are not reflected in the
level of the non-interacting Hamiltonian of bogolon. The effect of the odd-
frequency pairing is clearly seen in the single-particle spectral functions at
low energies where the disorder effect is dominantly present. The system with
Bogoliubov Fermi surfaces is a suitable playground for studying the properties
of odd-frequency pairing.

## Acknowledgment

This work was supported by JSPS KAKENHI Grants No. JP18K13490, No. JP18H01176,
No. JP18H04305 and No. JP19H01842.

## References

  * [1] G.E. Volovik, Phys. Lett. A 142, 282 (1989).
  * [2] G.E. Volovik, JETP Lett. 58, 469 (1993).
  * [3] W.V. Liu and F. Wilczek, Phys. Rev. Lett. 90, 047002 (2003).
  * [4] E. Gubankova, E.G. Mishchenko, and F. Wilczek, Phys. Rev. Lett. 94, 110402 (2005); Phys. Rev. B 74, 184516 (2006).
  * [5] D.F. Agterberg, P.M.R. Brydon, and C. Timm, Phys. Rev. Lett. 118, 127001 (2017).
  * [6] P.M.R. Brydon, D.F. Agterberg, H. Menke, and C. Timm, Phys. Rev. B 98, 224509 (2018).
  * [7] N.F.Q. Yuan and L. Fu, Phys. Rev. B 97, 115139 (2018).
  * [8] S. Sumita, T. Nomoto, K. Shiozaki, and Y. Yanase, Phys. Rev. B 99, 134513 (2019).
  * [9] H. Menke, C. Timm, and P.M.R. Brydon, Phys. Rev. B 100, 224505 (2019).
  * [10] J.M. Link, I. Boettcher, and I.F. Herbut, Phys. Rev. B 101, 184503 (2020).
  * [11] J.M. Link and I.F. Herbut, Phys. Rev. Lett. 125, 237004 (2020).
  * [12] S. Autti, J.T. Mäkinen, J. Rysti, G.E. Volovik, V.V. Zavjalov, and V.B. Eltsov, Phys. Rev. Res. 2, 033013 (2020).
  * [13] N.N. Bogoliubov, Sov. Phys. JETP 34, 41 (1958).
  * [14] C.J. Lapp, G. Börner, and C. Timm, Phys. Rev. B 101, 024505 (2020).
  * [15] C. Setty, Y. Cao, A. Kreisel, S. Bhattacharyya, and P.J. Hirschfeld, Phys. Rev. B 102, 064504 (2020).
  * [16] Y. Sato, S. Kasahara, T. Taniguchi, X. Xing, Y. Kasahara, Y. Tokiwa, Y. Yamakawa, H. Kontani, T. Shibauchi, and Y. Matsuda, Proc. Nat. Acad. Sci. 115 1227 (2018).
  * [17] C. Setty, S. Bhattacharyya, Y. Cao, A. Kreisel, and P.J. Hirschfeld, Nat. Commun. 11, 523 (2020).
  * [18] For a review, see T. Shibauchi, T. Hanaguri, and Y. Matsuda, J. Phys. Soc. Jpn. 89, 102002 (2020)
  * [19] H. Oh and E.-G. Moon, Phys. Rev. B 102, 020501(R) (2020).
  * [20] S.-T. Tamura, S. Iimura, and S. Hoshino, Phys. Rev. B 102, 024505 (2020).
  * [21] C. Timm, P.M.R. Brydon, and D.F. Agterberg, arXiv:2011.02798 (2020).
  * [22] I.F. Herbut and J.M. Link, arXiv:2012.10801 (2020).
  * [23] V.L. Berezinskii, JETP Lett. 20, 287 (1974).
  * [24] T.R. Kirkpatrick and D. Belitz, Phys. Rev. Lett. 66, 1533 (1991).
  * [25] D. Belitz and T.R. Kirkpatrick, Phys. Rev. B 46, 8393 (1992).
  * [26] A. Balatsky and E. Abrahams, Phys. Rev. B 45, 13125 (1992).
  * [27] V.J. Emery and S. Kivelson, Phys. Rev. B 46, 10812 (1992).
  * [28] J.R. Schrieffer, A.V. Balatsky, E. Abrahams, and D.J. Scalapino, J. Supercond. 7, 501 (1994).
  * [29] For a review, F.S. Bergeret, A.F. Volkov, and K.B. Efetov, Rev. Mod. Phys 77, 1321 (2005).
  * [30] For a review, Y. Tanaka, M. Sato, and N. Nagaosa, J. Phys. Soc. Jpn. 81, 011013 (2012).
  * [31] For a review, J. Linder and A.V. Balatsky, Rev. Mod. Phys. 91, 045005 (2019).
  * [32] A.V. Balatsky and J. Bon$\check{\rm c}$a, Phys. Rev. B 48, 7445 (1993).
  * [33] M. Vojta and E. Dagotto, Phys. Rev. B 59, 713(R) (1999).
  * [34] Y. Fuseya, H. Kohno, and K. Miyake, J. Phys. Soc. Jpn. 72, 2914 (2003).
  * [35] K. Yada, S. Onari, Y. Tanaka, and K. Miyake, arXiv:0806.4241 (2008).
  * [36] T. Hotta, J. Phys. Soc. Jpn. 78, 123710 (2009).
  * [37] K. Shigeta, S. Onari, K. Yada, and Y. Tanaka, Phys. Rev. B 79, 174507 (2009).
  * [38] K. Shigeta, Y. Tanaka, K. Kuroki, S. Onari, and H. Aizawa, ibid. 83, 140509 (2011).
  * [39] H. Kusunose, Y. Fuseya, and K. Miyake, J. Phys. Soc. Jpn. 80, 044711 (2011).
  * [40] Y. Yanagi, Y. Yamashita, and K. Ueda, J. Phys. Soc. Jpn. 81, 123701 (2012).
  * [41] V.J. Emery and S.A. Kivelson, Phys. Rev. Lett. 71, 3701 (1993).
  * [42] P. Coleman, E. Miranda, and A. Tsvelik, Phys. Rev. Lett. 70, 2960 (1993).
  * [43] P. Coleman, E. Miranda, and A. Tsvelik, Phys. Rev. B 49, 8955 (1994).
  * [44] O. Zachar, S.A. Kivelson, and V.J. Emery, Phys. Rev. Lett. 77, 1342 (1996).
  * [45] M. Jarrell, H. Pang, D.L. Cox, Phys. Rev. Lett. 78, 1996 (1997).
  * [46] F.B. Anders, Phys. Rev. B 66 020504(R) (2002).
  * [47] F.B. Anders, Eur. Phys. J. B 28, 9 (2002).
  * [48] S. Hoshino and Y. Kuramoto, Phys. Rev. Lett. 112, 167204 (2014).
  * [49] S. Hoshino, Phys. Rev. B 90, 115154 (2014).
  * [50] J. Otsuki, Phys. Rev. Lett. 115, 036404 (2015).
  * [51] H. Kusunose, J. Phys. Soc. Jpn. 85, 113701 (2016).
  * [52] A.M. Tsvelik, Proc. Nat. Acad. Sci. 116, 12729 (2019).
  * [53] R. Heid, Z. Phys. B 99, 15 (1995).
  * [54] D. Belitz and T. R. Kirkpatrick, Phys. Rev. B 60, 3485 (1999).
  * [55] D. Solenov, I. Martin, and D. Mozyrsky, Phys. Rev. B 79, 132502 (2009).
  * [56] H. Kusunose, Y. Fuseya, and K. Miyake, J. Phys. Soc. Jpn. 80, 054702 (2011).
  * [57] Ya.V. Fominov, Y. Tanaka, Y. Asano, and M. Eschrig, Phys. Rev. B 91, 144514 (2015).
  * [58] F.S. Bergeret, A.F. Volkov, and K.B. Efetov, Phys. Rev. Lett. 86, 4096 (2001).
  * [59] M. Eschrig, J. Kopu, J.C. Cuevas, and G. Schön, Phys. Rev. Lett. 90, 137003 (2003).
  * [60] Y. Tanaka and A. A. Golubov, Phys. Rev. Lett. 98, 037003 (2007).
  * [61] Y. Tanaka, Y. Tanuma, and A.A. Golubov, Phys. Rev. B 76, 054522 (2007).
  * [62] A.A. Abrikosov, L.P. Gorkov, and I.E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Dover, 1963).
  * [63] See Supplementary Materials.
  * [64] P. Coleman, A.M. Tsvelik, N. Andrei, and H.Y. Kee, Phys. Rev. B 60, 3608 (1999).
  * [65] R. Flint, M. Dzero, and P. Coleman, Nat. Phys. 4, 643 (2008).
  * [66] S. Iimura, M. Hirayama, and S. Hoshino, Phys. Rev. B 100, 094532 (2019).
  * [67] K. Miyake, T. Matsuura, and C.M. Varma, Solid State Commun.71, 1149 (1989).
  * [68] A.C. Jacko, J.O. Fjærestad, and B.J. Powell, Nat. Phys. 5, 422 (2009).

SUPPLEMENTARY MATERIAL FOR
“Odd-frequency pairing inherent in Bogoliubov Fermi liquid”
T. Miki, S.-T. Tamura, S. Iimura, and S. Hoshino
(Dated: )

## SM A. Impurity effects

We consider the impurity potential with the form

$\displaystyle\mathscr{H}_{\rm imp}$
$\displaystyle=\frac{1}{V}\sum_{\bm{k},\bm{q}}\rho_{\bm{q}}u_{1}(\bm{k},\bm{q})\alpha^{\dagger}_{\bm{k}+\bm{q}}\alpha_{\bm{k}}$
$\displaystyle+\frac{1}{V}\sum_{\bm{k},\bm{q}}\rho_{\bm{q}}u_{2}(\bm{k},\bm{q})\alpha^{\dagger}_{\bm{k}+\bm{q}}\alpha^{\dagger}_{-\bm{k}}+{\rm
H.c.},$ (S1)

which breaks the translational symmetry or momentum conservation. The
information for the disorder is included in the structure factor

$\displaystyle\rho_{\bm{q}}=\sum_{i}\mathrm{e}^{-\mathrm{i}\bm{q}\cdot\bm{R}_{i}},$
(S2)

where the summation with respect to $i$ is performed for the impurity
positions $\bm{R}_{i}$, and $n_{\mathrm{imp}}=V^{-1}\sum_{i}1$. This quantity
satisfies

$\displaystyle\overline{\rho_{\bm{q}}\rho_{\bm{q}^{\prime}}}=Vn_{\mathrm{imp}}\delta_{\bm{q},-\bm{q}^{\prime}},$
(S3)

where the overline indicates the average over the impurity configurations. The
connection to the original electrons is shown in SM C. We have the relations

$\displaystyle u_{1}(\bm{k},\bm{q})=u_{1}^{*}(\bm{k}+\bm{q},-\bm{q}),$ (S4)
$\displaystyle u_{2}(\bm{k},\bm{q})=-u_{2}(-\bm{k}-\bm{q},\bm{q}),$ (S5)
$\displaystyle u_{1,2}(\bm{k},\bm{q})=u_{1,2}(-\bm{k},-\bm{q}),$ (S6)

each of which originates from the Hermiticity, Pauli principle
(anticommutation relation), and the inversion symmetry, respectively.

The single-particle Green function is given by the $2\times 2$ matrix as

$\displaystyle\hat{G}_{\bm{k}}(\tau)$
$\displaystyle=-\langle\mathcal{T}\begin{pmatrix}\alpha_{\bm{k}}(\tau)\\\
\alpha^{\dagger}_{-\bm{k}}(\tau)\end{pmatrix}\begin{pmatrix}\alpha_{\bm{k}}^{\dagger}&\alpha_{-\bm{k}}\end{pmatrix}\rangle,$
(S7)

where $\mathcal{T}$ indicates a imaginary time ordering, and
$A(\tau)=\mathrm{e}^{\tau\mathscr{H}}A\mathrm{e}^{-\tau\mathscr{H}}$ is the
Heisenberg picture. The Fourier transformation is defined by

$\displaystyle\hat{G}_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=\int_{0}^{1/T}\mathrm{d}\tau\hat{G}_{\bm{k}}(\tau)\mathrm{e}^{\mathrm{i}\varepsilon_{n}\tau}.$
(S8)

As given in the main text, the normal and anomalous self-energies are
introduced by

$\displaystyle\hat{G}_{\bm{k}}^{-1}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=\begin{pmatrix}\mathrm{i}\varepsilon_{n}-\xi_{\bm{k}}-\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})&-S_{\bm{k}}(\mathrm{i}\varepsilon_{n})\\\
-S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})&\mathrm{i}\varepsilon_{n}+\xi_{\bm{k}}-\bar{\Sigma}_{\bm{k}}(\mathrm{i}\varepsilon_{n})\end{pmatrix}.$
(S9)

In the following of this section, we will derive explicit form of self-
energies.

### 1\. Born approximation

We consider the lowest-order contribution of the impurity potential. Second
order self-energies (Born approximation) are given as follows:

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})+2\Sigma_{2,\bm{k}}(\mathrm{i}\varepsilon_{n}),$
(S10) $\displaystyle\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=n_{\mathrm{imp}}\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}|u_{1}(\bm{q},\bm{k}-\bm{q})|^{2}G^{0}_{\bm{q}}(\mathrm{i}\varepsilon_{n}),$
(S11) $\displaystyle\Sigma_{2,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-2n_{\mathrm{imp}}\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}|u_{2}(\bm{q},\bm{k}-\bm{q})|^{2}G^{0}_{\bm{q}}(-\mathrm{i}\varepsilon_{n}),$
(S12) $\displaystyle S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-2n_{\mathrm{imp}}\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}u_{1}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})$
$\displaystyle\quad\times[G^{0}_{\bm{q}}(-\mathrm{i}\varepsilon_{n})-G^{0}_{\bm{q}}(\mathrm{i}\varepsilon_{n})],$
(S13)

which corresponds to the diagrams shown in Fig. 2 of the main text. For the
evaluation of the wavevector integral, we replace the coupling constant by its
averaged value with respect to $\bm{q}$. Then we only have to evaluate the
integration as follows:

$\displaystyle\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}G^{0}_{\bm{q}}(\pm\mathrm{i}\varepsilon_{n})$
$\displaystyle=\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}\frac{1}{\pm\mathrm{i}\varepsilon_{n}-\xi_{\bm{q}}}$
$\displaystyle=\mp\mathrm{i}\pi D_{0}\mathrm{sgn}\,\varepsilon_{n},$ (S14)

where $D_{0}$ is a density of states of bogolon at Fermi level. Then,

$\displaystyle\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\mathrm{i}\pi
n_{\mathrm{imp}}D_{0}\mathrm{sgn}\,\varepsilon_{n}\langle|u_{1}(\bm{q},\bm{k}-\bm{q})|^{2}\rangle_{\bm{q}},$
(S15) $\displaystyle\Sigma_{2,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-2\mathrm{i}\pi
n_{\mathrm{imp}}D_{0}\mathrm{sgn}\,\varepsilon_{n}\langle|u_{2}(\bm{q},\bm{k}-\bm{q})|^{2}\rangle_{\bm{q}},$
(S16) $\displaystyle S_{\bm{k}}^{\dagger}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-4\mathrm{i}\pi
n_{\mathrm{imp}}D_{0}\mathrm{sgn}\,\varepsilon_{n}\langle
u_{1}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})\rangle_{\bm{q}},$
(S17)

where
$\langle\cdots\rangle_{\bm{q}}=\int\mathrm{d}\bm{q}\cdots/\int\mathrm{d}\bm{q}1$.
Hence, we obtain the self-energies

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\mathrm{i}\Gamma_{1,\bm{k}}\mathrm{sgn}\,\varepsilon_{n},$
(S18) $\displaystyle S_{\bm{k}}^{\dagger}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\mathrm{i}\Gamma_{2,\bm{k}}\mathrm{sgn}\,\varepsilon_{n},$
(S19)

where

$\displaystyle\Gamma_{1,\bm{k}}$ $\displaystyle=\pi
n_{\mathrm{imp}}D_{0}\big{(}\langle|u_{1}(\bm{q},\bm{k}-\bm{q})|^{2}\rangle_{\bm{q}}+4\langle|u_{2}(\bm{q},\bm{k}-\bm{q})|^{2}\rangle_{\bm{q}}\big{)},$
(S20) $\displaystyle\Gamma_{2,\bm{k}}$ $\displaystyle=4\pi
n_{\mathrm{imp}}D_{0}\langle
u_{1}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})\rangle_{\bm{q}}.$
(S21)

One sees that the usual Born approximation result is reproduced if $u_{2}$,
which is specific to bogolons, is set to zero.

As pointed out in the main text, the spectral function shows unphysical
behavior for $|\Gamma_{2,\bm{k}}|>\Gamma_{1,\bm{k}}$. However, we can show
that the condition $|\Gamma_{2,\bm{k}}|<\Gamma_{1,\bm{k}}$ is satisfied based
on the expressions obtained above together with the magnitude relation between
arithmetic and geometric means.

### 2\. Self-consistent treatment

Figure S1:  Feynman diagrams relevant to the self-consistent Born
approximation.

We also consider the self-energy with self-consistent treatment. Namely, we
consider the diagrams for the self-energy shown in Fig. S1. The self-
consistent equations are derived as

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})=n_{\mathrm{imp}}\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}\Big{[}|u_{1}(\bm{q},\bm{k}-\bm{q})|^{2}G_{\bm{q}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle-4|u_{2}(\bm{q},\bm{k}-\bm{q})|^{2}G_{\bm{q}}(-\mathrm{i}\varepsilon_{n})$
$\displaystyle-2u_{1}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})F_{\bm{q}}(-\mathrm{i}\varepsilon_{n})$
$\displaystyle-2u_{1}^{\ast}(\bm{q},\bm{k}-\bm{q})u_{2}(\bm{q},\bm{k}-\bm{q})F^{\dagger}_{\bm{q}}(-\mathrm{i}\varepsilon_{n})\Big{]},$
(S22) $\displaystyle
S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})=n_{\mathrm{imp}}\int\frac{\mathrm{d}\bm{q}}{(2\pi)^{3}}\Big{[}2u_{1}^{\ast}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})$
$\displaystyle\times\big{(}G_{\bm{q}}(\mathrm{i}\varepsilon_{n})-G_{\bm{q}}(-\mathrm{i}\varepsilon_{n})\big{)}$
$\displaystyle+\frac{1}{2}u_{1}^{\ast}(\bm{q},\bm{k}-\bm{q})u_{1}^{\ast}(\bm{q},\bm{k}-\bm{q})\big{(}F^{\dagger}_{\bm{q}}(\mathrm{i}\varepsilon_{n})-F^{\dagger}_{\bm{q}}(-\mathrm{i}\varepsilon_{n})\big{)}$
$\displaystyle+2u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})u_{2}^{\ast}(\bm{q},\bm{k}-\bm{q})\big{(}F_{\bm{q}}(\mathrm{i}\varepsilon_{n})-F_{\bm{q}}(-\mathrm{i}\varepsilon_{n})\big{)}\Big{]}.$

Obviously, these equations reduce to those in the Born approximation in the
last subsection if one drops the self-energies in the right-hand side. In
order to search for the concrete solutions, we simplify the equation by
replacing the coefficients by the wavevector independent ones, and then the
self-energies are also $\bm{k}$-independent. We have considered all the
parameter space within this approximation, and always have found physical
solutions. Defining the energy-dependent functions
$\Gamma_{1,2}(\mathrm{i}\varepsilon_{n})$ by
$\Sigma(\mathrm{i}\varepsilon_{n})=-\mathrm{i}\Gamma_{1}(\mathrm{i}\varepsilon_{n})$
and
$S^{\dagger}(\mathrm{i}\varepsilon_{n})=-\mathrm{i}\Gamma_{2}(\mathrm{i}\varepsilon_{n})$,
we obtain the equations

$\displaystyle\Gamma_{1}(\mathrm{i}\varepsilon_{n})=\frac{\mathrm{i}\pi
n_{\mathrm{imp}}D_{0}}{\sqrt{-(\varepsilon_{n}+\mathrm{Re}\Gamma_{1}(\mathrm{i}\varepsilon_{n}))^{2}+|\Gamma_{2}(\mathrm{i}\varepsilon_{n})|^{2}}}$
$\displaystyle\times\Big{[}\big{(}|u_{1}|^{2}+4|u_{2}|^{2}\big{)}(\varepsilon_{n}+\mathrm{Re}\Gamma_{1}(\mathrm{i}\varepsilon_{n}))$
$\displaystyle-2u_{1}u_{2}^{\ast}\Gamma_{2}^{\ast}(i\varepsilon_{n})-2u_{1}^{\ast}u_{2}\Gamma_{2}(\mathrm{i}\varepsilon_{n})\Big{]}$
(S24) $\displaystyle\Gamma_{2}(\mathrm{i}\varepsilon_{n})=\frac{\mathrm{i}\pi
n_{\mathrm{imp}}D_{0}}{\sqrt{-(\varepsilon_{n}+\mathrm{Re}\Gamma_{1}(\mathrm{i}\varepsilon_{n}))^{2}+|\Gamma_{2}(\mathrm{i}\varepsilon_{n})|^{2}}}$
$\displaystyle\times\Big{[}4u_{1}^{\ast}u_{2}^{\ast}(\varepsilon_{n}+\mathrm{Re}\Gamma_{1}(\mathrm{i}\varepsilon_{n}))$
$\displaystyle-
u_{1}^{\ast}u_{1}^{\ast}\Gamma_{2}(i\varepsilon_{n})-4u_{2}^{\ast}u_{2}^{\ast}\Gamma_{2}^{\ast}(\mathrm{i}\varepsilon_{n})\Big{]}.$
(S25)

Performing a suitable transformation for the dimensionless expressions of the
equations, we find that the equations are controlled by the complex parameter
$u_{1}/|u_{2}|$. We note that the solution is identical to the Born
approximation discussed in the last subsection at high frequencies.

Figure S2:  Examples of the solutions of Eqs. (S24) and (S25). We define
dimensionless self-energies and frequency as $\tilde{\Sigma}=\Sigma/(\pi
n_{\mathrm{imp}}D_{0}|u_{2}|^{2}),\tilde{S}^{\dagger}=S^{\dagger}/(\pi
n_{\mathrm{imp}}D_{0}|u_{2}|^{2}),\tilde{\varepsilon}_{n}=\varepsilon_{n}/(\pi
n_{\mathrm{imp}}D_{0}|u_{2}|^{2})$. The parameters are chosen as (a)
$u_{1}/|u_{2}|=3\exp(\mathrm{i}\pi/6),\arg u_{2}=\pi/12$ and (b)
$u_{1}/|u_{2}|=3\exp(\mathrm{i}\pi/3),\arg u_{2}=\pi/12$. The inset shows the
self-energy with $\varepsilon_{n}$ multiplied, indicating
$-\varepsilon_{n}\mathrm{Im}\,\Sigma(\mathrm{i}\varepsilon_{n})=\varepsilon_{n}|S^{\dagger}(\mathrm{i}\varepsilon_{n})|$
for $\varepsilon_{n}\to 0$.

We can solve the above simultaneous equations at each frequency. The
calculation examples are shown in Fig. S2, where we have the two kinds of
solutions depending on the parameters both of which are odd function in
frequency or imaginary time. The characteristic behaviors are seen in the low-
energy limit $\varepsilon_{n}\to 0$. The first kind has the form

$\displaystyle\Sigma(\mathrm{i}\varepsilon_{n}\to 0)$
$\displaystyle=-\mathrm{i}\Gamma_{1}\mathrm{sgn}\,\varepsilon_{n},$ (S26)
$\displaystyle S^{\dagger}(\mathrm{i}\varepsilon_{n}\to 0)$
$\displaystyle=-\mathrm{i}\Gamma_{2}\mathrm{sgn}\,\varepsilon_{n},$ (S27)

which is similar to the results in the simple Born approximation. We can
explicitly obtain the solution as

$\displaystyle\frac{\Gamma_{2}}{\Gamma_{1}}=\begin{cases}\displaystyle\frac{2u_{2}^{\ast}}{\mathrm{Re}u_{1}}&(2<|\mathrm{Re}u_{1}|/|u_{2}|),\vspace{5pt}\\\
\displaystyle\frac{\mathrm{Re}u_{1}}{2u_{2}}&(\mathrm{Im}u_{1}/|u_{2}|=0\text{
and }0<|\mathrm{Re}u_{1}|/|u_{2}|<2)\end{cases}$ (S28)

for $|\mathrm{Re}\,u_{1}|/|u_{2}|>2$ or $\mathrm{Im}\,u_{1}/|u_{2}|=0$. The
calculation example is shown in Fig. S2(a).

On the other hand, the solution of the second kind is obtained for
$|\mathrm{Re}\,u_{1}|/|u_{2}|<2$ and $\mathrm{Im}\,u_{1}/|u_{2}|\neq 0$, which
has the frequency dependence

$\displaystyle\Sigma(\mathrm{i}\varepsilon_{n}\to 0)$
$\displaystyle=\frac{|V|^{2}}{\mathrm{i}\varepsilon_{n}},$ (S29)
$\displaystyle S^{\dagger}(\mathrm{i}\varepsilon_{n}\to 0)$
$\displaystyle=\frac{V^{2}}{\mathrm{i}\varepsilon_{n}},$ (S30)

where $V$ is a complex constant. The calculation example is shown in Fig.
S2(b). The absolute value $|V|$ cannot be written in the simple analytic form,
but is determined easily from the numerical calculations.

### 3\. Single-particle spectral functions

From the normal and anomalous Green functions, we obtain several physical
quantities of interest. The first one is the single-particle spectral function
and density of states defined by

$\displaystyle A_{\bm{k}}(\varepsilon)$
$\displaystyle=-\frac{1}{\pi}\mathrm{Im}\,G_{\bm{k}11}(\varepsilon+\mathrm{i}\eta),$
(S31) $\displaystyle D(\varepsilon)$
$\displaystyle=\sum_{\bm{k}}A_{\bm{k}}(\varepsilon)\simeq
D_{0}\int\mathrm{d}\xi_{\bm{k}}A_{\bm{k}}(\varepsilon).$ (S32)

The momentum distribution function is given by

$\displaystyle n_{\bm{k}}$
$\displaystyle=\langle\alpha^{\dagger}_{\bm{k}}\alpha_{\bm{k}}\rangle.$ (S33)

In addition, the time-dependent pair amplitude can also be calculated. Whereas
the static quantity such as $\langle\alpha_{\bm{k}}\alpha_{-\bm{k}}\rangle$ is
zero due to the inversion symmetry, the pair amplitude is finite at different
times. Specifically we consider the spatially local quantity

$\displaystyle F(\tau)$
$\displaystyle=T\sum_{n,\bm{k}}F_{\bm{k}}(\mathrm{i}\varepsilon_{n})\mathrm{e}^{-\mathrm{i}\varepsilon_{n}\tau}.$
(S34)

The above quantities are graphically shown in the main text.

The real time dynamics is also obtained by performing the analytic
continuation. The retarded Green function is given by

$\displaystyle G^{\rm
R}(t)=-\mathrm{i}\theta(t)\langle\\{\alpha_{\bm{k}}(t),\alpha_{\bm{k}}^{\dagger}\\}\rangle,$
(S35)

where $\theta(t)$ is the step function. Here we have considered the real-time
Heisenberg picture
$A(t)=\mathrm{e}^{\mathrm{i}t\mathscr{H}}A\mathrm{e}^{-\mathrm{i}t\mathscr{H}}$.
In the following, we summarize the above physical quantities for the solutions
of both the first and second kinds.

#### (i) Solution of the first kind

We consider the self-energies

$\displaystyle\Sigma(\mathrm{i}\varepsilon_{n})=-\mathrm{i}\Gamma_{1}\mathrm{sgn}\,\varepsilon_{n},\
\
S^{\dagger}(\mathrm{i}\varepsilon_{n})=-\mathrm{i}\Gamma_{2}\mathrm{sgn}\,\varepsilon_{n},$
(S36)

with $\Gamma_{1}>|\Gamma_{2}|$. The density of states is

$\displaystyle\frac{D(\varepsilon)}{D_{0}}$
$\displaystyle=\frac{|\varepsilon|+\mathrm{i}\Gamma_{1}}{2\sqrt{(|\varepsilon|+\mathrm{i}\Gamma_{1})^{2}+|\Gamma_{2}|^{2}}}+{\rm
c.c.}$ (S37)

The momentum distribution function is

$\displaystyle n_{\bm{k}}$
$\displaystyle=\frac{\xi_{\bm{k}}}{2\pi\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}}\ln\left(\frac{\Gamma_{1}-\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}}{\Gamma_{1}+\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}}\right)+\frac{1}{2}$
(S38)

at zero temperature. The spatially local pair amplitude at $T=0$ is

$\displaystyle\frac{F^{\dagger}(\tau)}{\Gamma_{2}}$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}\varepsilon\left[\frac{\mathrm{i}}{2\sqrt{(|\varepsilon|+\mathrm{i}\Gamma_{1})^{2}+|\Gamma_{2}|^{2}}}+{\rm
c.c.}\right]\mathrm{e}^{-\varepsilon\tau}$ (S39)

for $\tau>0$, which is numerically integrated.

The real time dependence of the retarded Green function ($t>0$) is given by

$\displaystyle G^{\rm R}_{\bm{k}}(t)$
$\displaystyle=-\frac{\mathrm{i}}{2}\left(1+\frac{\xi_{\bm{k}}}{\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}}\right)\mathrm{e}^{-\mathrm{i}\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}\,t-\Gamma_{1}t}$
$\displaystyle-\frac{\mathrm{i}}{2}\left(1-\frac{\xi_{\bm{k}}}{\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}}\right)\mathrm{e}^{\mathrm{i}\sqrt{\xi_{\bm{k}}^{2}-|\Gamma_{2}|^{2}}\,t-\Gamma_{1}t}$
(S40)

for $|\xi_{\bm{k}}|>|\Gamma_{2}|$ and

$\displaystyle G^{\rm R}_{\bm{k}}(t)$
$\displaystyle=-\frac{1}{2}\left(\mathrm{i}+\frac{\xi_{\bm{k}}}{\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}}\right)\mathrm{e}^{-(\Gamma_{1}-\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}})t}$
$\displaystyle-\frac{1}{2}\left(\mathrm{i}-\frac{\xi_{\bm{k}}}{\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}}}\right)\mathrm{e}^{-(\Gamma_{1}+\sqrt{|\Gamma_{2}|^{2}-\xi_{\bm{k}}^{2}})t}$
(S41)

for $|\xi_{\bm{k}}|<|\Gamma_{2}|$.

#### (ii) Solution of the second kind

We summarize the results for the self-energies

$\displaystyle\Sigma(\mathrm{i}\varepsilon_{n})=|V|^{2}/\mathrm{i}\varepsilon_{n},\
\ S^{\dagger}(\mathrm{i}\varepsilon_{n})=V^{2}/\mathrm{i}\varepsilon_{n}.$
(S42)

The density of states is

$\displaystyle\frac{D(\varepsilon)}{D_{0}}$
$\displaystyle=\frac{\pi}{\sqrt{2}}|V|\delta(\varepsilon)+\frac{\varepsilon^{2}-|V|^{2}}{\varepsilon\sqrt{\varepsilon^{2}-2|V|^{2}}}\theta(\varepsilon-\sqrt{2}|V|).$
(S43)

The momentum distribution function is

$\displaystyle n_{\bm{k}}$
$\displaystyle=\frac{E_{\bm{k}}-\xi_{\bm{k}}[1-2f(E_{\bm{k}})]}{2E_{\bm{k}}},$
(S44)

where $E_{\bm{k}}=\sqrt{\xi_{\bm{k}}^{2}+2|V|^{2}}$ and
$f(x)=1/(\mathrm{e}^{x/T}+1)$. The spatially local pair amplitude is

$\displaystyle\frac{F^{\dagger}(\tau)}{V}$ $\displaystyle=\frac{\pi
V}{2\sqrt{2}|V|}-V\int_{\sqrt{2}|V|}^{\infty}\mathrm{d}\varepsilon\frac{\mathrm{e}^{-\varepsilon\tau}}{\varepsilon\sqrt{\varepsilon^{2}-2|V|^{2}}}$
(S45)

for $\tau>0$, which is numerically integrated. The real time dependence of the
Green function is given by

$\displaystyle G^{\rm R}_{\bm{k}}(t)$
$\displaystyle=-\frac{\mathrm{i}}{2E_{\bm{k}}^{2}}\Big{[}2|V|^{2}+(\xi_{\bm{k}}^{2}+|V|^{2}+\xi_{\bm{k}}E_{\bm{k}})\mathrm{e}^{-\mathrm{i}E_{\bm{k}}t}$
$\displaystyle+(\xi_{\bm{k}}^{2}+|V|^{2}-\xi_{\bm{k}}E_{\bm{k}})\mathrm{e}^{\mathrm{i}E_{\bm{k}}t}\Big{]}$
(S46)

for $t>0$.

## SM B. Interaction with bosons

The interaction Hamiltonian is defined as Eq (2) and Eq (3) in the main text.
With the Hermiticity and inversion symmetry, we obtain the relations similar
to Eqs. (S4–S6) for $g_{1,2}(\bm{k},\bm{q})$. The derivation of the self-
energies is similar to the standard procedure used in the Fermi liquid theory
[62].

Figure S3:  (a) Self-energy for bosons. (b) Wavevector dependence of $h_{1,2}$
defined in Eqs. (S55) and (S59) at several temperatures. The cutoff wavevector
is chosen as $\omega_{\mathrm{C}}/\varepsilon_{\mathrm{F}b}=0.3$

### 1\. Boson self-energy

We first derive the boson self-energy. The Green function of boson is defined
by

$\displaystyle
D_{\bm{q}}(\tau)=-\langle\mathcal{T}\phi_{\bm{q}}(\tau)\phi_{-\bm{q}}\rangle,$
(S47)
$\displaystyle\phi_{\bm{q}}=\mathrm{i}\sqrt{\frac{\omega_{0,\bm{q}}}{2}}(b_{\bm{q}}-b_{-\bm{q}}^{\dagger}).$
(S48)

Fourier transformation is given by

$\displaystyle D_{\bm{q}}(\mathrm{i}\omega_{n})$
$\displaystyle=\int_{0}^{1/T}\mathrm{d}\tau
D_{\bm{q}}(\tau)\mathrm{e}^{\mathrm{i}\omega_{n}\tau},$ (S49)

where $\omega_{n}=2n\pi T$ is the bosonic Matsubara frequency. The self-energy
$\Pi_{\bm{q}}(\mathrm{i}\omega_{n})$ is defined by

$\displaystyle D_{\bm{q}}(\mathrm{i}\omega_{n})^{-1}$
$\displaystyle=D^{0}_{\bm{q}}(\mathrm{i}\omega_{n})^{-1}-\Pi_{\bm{q}}(\mathrm{i}\omega_{n}),$
(S50)

with the zeroth-order Green function
$D^{0}_{\bm{q}}(\mathrm{i}\omega_{n})=\omega_{0,\bm{q}}^{2}/[(\mathrm{i}\omega_{n})^{2}-\omega_{0,\bm{q}}^{2}]$.
The diagrammatic contributions are illustrated in Fig. S3(a), and are then
given as follows:

$\displaystyle\Pi_{\bm{q}}(\mathrm{i}\omega_{n})=\Pi_{1,\bm{q}}(\mathrm{i}\omega_{n})+2\Pi_{2,\bm{q}}(\mathrm{i}\omega_{n}),$
(S51)
$\displaystyle\Pi_{1,\bm{q}}(\mathrm{i}\omega_{n})=\frac{T}{(2\pi)^{3}}\sum_{m}\int\mathrm{d}\bm{k}\,|g_{1}(\bm{k},-\bm{q})|^{2}$
$\displaystyle\quad\times
G^{0}_{\bm{k}}(\mathrm{i}\varepsilon_{m})G^{0}_{\bm{k}-\bm{q}}(\mathrm{i}\varepsilon_{m}-\mathrm{i}\omega_{n}),$
(S52)
$\displaystyle\Pi_{2,\bm{q}}(\mathrm{i}\omega_{n})=-\frac{T}{(2\pi)^{3}}\sum_{m}\int\mathrm{d}\bm{k}|g_{2}(\bm{k},-\bm{q})|^{2}G^{0}_{\bm{k}}(\mathrm{i}\varepsilon_{m})$
$\displaystyle\quad\times[G^{0}_{-\bm{k}+\bm{q}}(-\mathrm{i}\varepsilon_{m}+\mathrm{i}\omega_{n})+G^{0}_{-\bm{k}+\bm{q}}(-\mathrm{i}\varepsilon_{m}-\mathrm{i}\omega_{n})],$
(S53)

where $G^{0}_{\bm{k}}(\mathrm{i}\varepsilon_{n})$ is a free particle Green
function of bogolon. First, we find an explicit form of
$\Pi_{1,\bm{q}}(\mathrm{i}\omega_{n})$. For the concrete calculation, we
introduce the wavevector-averaged coupling constant
$\overline{\Gamma}_{1,\bm{q}}\equiv\langle|g_{1}(\bm{k},-\bm{q})|^{2}\rangle_{\bm{k}}$.
Then we obtain the real part of retarded self-energy at low $\omega$ as
follows:

$\displaystyle\mathrm{Re}\Pi^{\mathrm{R}}_{1,\bm{q}}(\omega)=-\frac{\overline{\Gamma}_{1,\bm{q}}m_{b}k_{\mathrm{F}b}}{(2\pi)^{2}}h_{1}\left(\frac{q}{2k_{\mathrm{F}b}}\right),$
(S54)

where

$\displaystyle h_{1}(x)\equiv
1+\frac{1-x^{2}}{2x}\ln\left|\frac{1+x}{1-x}\right|$ (S55)

is the Lindhard function shown in Fig. S3(b). We also get the imaginary part
as

$\displaystyle\mathrm{Im}\Pi^{\mathrm{R}}_{1,\bm{q}}(\omega)=-\frac{\overline{\Gamma}_{1,\bm{q}}m_{b}^{2}}{2\pi
q}\omega\theta(2k_{\mathrm{F}b}-q).$ (S56)

These are the same form as those for the usual electrons. On the other hand,
in order to calculate $\Pi^{\mathrm{R}}_{2,\bm{q}}(\omega)$, we introduce the
wavevector-averaged coupling constant
$\overline{\Gamma}_{2,\bm{q}}\equiv\langle|g_{2}(\bm{k},-\bm{q})|^{2}\rangle_{\bm{k}}$.
The result of the calculation of $\Pi^{\mathrm{R}}_{2,\bm{q}}(\omega)$ at low
$\omega$ is as follows:

$\displaystyle\mathrm{Re}\Pi^{\mathrm{R}}_{2,\bm{q}}(\omega)$
$\displaystyle=-\frac{\overline{\Gamma}_{2,\bm{q}}m_{b}k_{\mathrm{F}b}}{(2\pi)^{2}}h_{2,T}\left(\frac{q}{2k_{\mathrm{F}b}}\right),$
(S57) $\displaystyle\mathrm{Im}\Pi^{\mathrm{R}}_{2,\bm{q}}(\omega)$
$\displaystyle=-\frac{\overline{\Gamma}_{2,\bm{q}}m_{b}^{2}}{2\pi
q}\omega\theta\big{(}(1+\sqrt{3})k_{\mathrm{F}b}-q\big{)},$ (S58)

where

$\displaystyle h_{2,T}(x)\equiv$
$\displaystyle\int_{\sqrt{1-\omega_{\mathrm{C}}/\varepsilon_{\mathrm{F}b}}}^{\sqrt{1+\omega_{\mathrm{C}}/\varepsilon_{\mathrm{F}b}}}\mathrm{d}y\frac{y}{x}\tanh\left(\frac{y^{2}-1}{2T/\varepsilon_{\mathrm{F}b}}\right)\ln\left|\frac{y^{2}-1+2x^{2}+2xy}{y^{2}-1+2x^{2}-2xy}\right|,$
(S59)

where the energy integration is performed within the cutoff frequency of
bosons $\omega_{\mathrm{C}}$. We have assumed that the dominant contributions
enter at small $q$. Note that $h_{2,T}(q/2k_{\mathrm{F}b})$ diverges
logarithmically in the infrared regime at zero temperature limit. This
corresponds to the Cooper instability intrinsic to Fermi surfaces. However,
this divergence is suppressed at finite $T$ as shown in Fig. S3(b), in which
we are interested.

For later discussion, we derive the explicit form of boson Green function. The
Dyson equation is given as follows:

$\displaystyle[D^{\mathrm{R}}_{\bm{q}}(\omega)]^{-1}=\frac{1}{\omega_{0,\bm{q}}^{2}}\left(\omega^{2}-\omega_{\bm{q}}^{2}-2\mathrm{i}\omega\gamma_{\bm{q}}\right),$
(S60)

where

$\displaystyle\omega_{\bm{q}}$
$\displaystyle=\omega_{0,\bm{q}}\sqrt{1-\eta_{\bm{q}}},$ (S61)
$\displaystyle\eta_{\bm{q}}$
$\displaystyle=\frac{\overline{\Gamma}_{1,\bm{q}}m_{b}k_{\mathrm{F}b}}{(2\pi)^{2}}h_{1}\left(\frac{q}{2k_{\mathrm{F}b}}\right)+2\frac{\overline{\Gamma}_{2,\bm{q}}m_{b}k_{\mathrm{F}b}}{(2\pi)^{2}}h_{2,T}\left(\frac{q}{2k_{\mathrm{F}b}}\right),$
(S62) $\displaystyle\gamma_{\bm{q}}$
$\displaystyle=\frac{\overline{\Gamma}_{1,\bm{q}}m_{b}^{2}\omega_{0,\bm{q}}^{2}}{4\pi
q}\theta(2k_{\mathrm{F}b}-q)$
$\displaystyle\quad+\frac{\overline{\Gamma}_{2,\bm{q}}m_{b}^{2}\omega_{0,\bm{q}}^{2}}{2\pi
q}\theta\big{(}(1+\sqrt{3})k_{\mathrm{F}b}-q\big{)}.$ (S63)

Therefore, we get the explicit form of the Green function as

$\displaystyle
D_{\bm{q}}^{\mathrm{R}}(\omega)=\frac{\omega_{0,\bm{q}}^{2}}{2\omega_{\bm{q}}}\left(\frac{1}{\omega-\omega_{\bm{q}}+\mathrm{i}\gamma_{\bm{q}}}-\frac{1}{\omega+\omega_{\bm{q}}+\mathrm{i}\gamma_{\bm{q}}}\right).$
(S64)

which will be used to obtain the bogolon self-energies.

### 2\. Bogolon normal self-energy

Next, we derive the bogolon normal self-energy given in Eq. (S9). The diagram
of the second-order self-energy is shown in Fig. 2(b) of the main text. The
corresponding self-energy $\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$ is
given as follows:

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})+2\Sigma_{2,\bm{k}}(\mathrm{i}\varepsilon_{n}),$
(S65) $\displaystyle\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=-\frac{T}{(2\pi)^{3}}\sum_{m}\int\mathrm{d}\bm{k}^{\prime}\,|g_{1}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})|^{2}$
$\displaystyle\quad\times
G^{0}_{\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{m})D_{\bm{k}-\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{n}-\mathrm{i}\varepsilon_{m}),$
(S66) $\displaystyle\Sigma_{2,\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=\frac{2T}{(2\pi)^{3}}\sum_{m}\int\mathrm{d}\bm{k}^{\prime}\,|g_{2}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})|^{2}$
$\displaystyle\quad\times
G^{0}_{\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{m})D_{\bm{k}-\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{n}+\mathrm{i}\varepsilon_{m}).$
(S67)

We calculate $\Sigma_{1,\bm{k}}(\mathrm{i}\varepsilon_{n})$ at first. For the
evaluation of $\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$, we replace the
Matsubara sums by the energy integral. By using the spectral representation
for the non-interacting Green function

$\displaystyle
G^{0,\mathrm{R}}(\bm{k},\varepsilon)=\int_{-\infty}^{\infty}\mathrm{d}\varepsilon^{\prime}\frac{A_{\bm{k}}^{0}(\varepsilon^{\prime})}{\varepsilon+\mathrm{i}\eta-\varepsilon^{\prime}},$
(S68)

the retarded self-energy has the form

$\displaystyle\Sigma^{\mathrm{R}}_{1,\bm{k}}(\varepsilon)=\frac{1}{(2\pi)^{4}}\int\mathrm{d}\bm{k}^{\prime}\int_{-\infty}^{\infty}\mathrm{d}\varepsilon^{\prime}\int_{-\infty}^{\infty}\mathrm{d}\omega$
$\displaystyle\quad\times|g_{1}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})|^{2}\frac{A_{\bm{k}^{\prime}}^{0}(\varepsilon^{\prime})\mathrm{Im}D^{\mathrm{R}}_{\bm{k}-\bm{k}^{\prime}}(\omega)}{\omega+\varepsilon^{\prime}-\varepsilon-\mathrm{i}\eta}$
$\displaystyle\quad\times\left[\tanh\left(\frac{\varepsilon^{\prime}}{2T}\right)+\coth\left(\frac{\omega}{2T}\right)\right].$
(S69)

To evaluate the integral, we again introduce the wavevector-averaged coupling
constant
$\overline{\Gamma}_{3,\bm{k}}\equiv\langle|g_{1}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})|^{2}\rangle_{\bm{k}^{\prime}}\in\mathbb{R}$.
Using the variable transformation
$\bm{k}^{\prime}\to(q,\xi_{\bm{k}^{\prime}})$ with
$\bm{q}=\bm{k}-\bm{k}^{\prime}$, we obtain

$\displaystyle\Sigma^{\mathrm{R}}_{1,\bm{k}}(\varepsilon)=\frac{\overline{\Gamma}_{3,\bm{k}}m_{b}}{(2\pi)^{3}k}\int_{0}^{k_{1}}\mathrm{d}q\,q\int_{-\infty}^{\infty}\mathrm{d}\omega\int_{\xi_{|\bm{k}|-q}}^{\xi_{|\bm{k}|+q}}\mathrm{d}\varepsilon^{\prime}$
$\displaystyle\quad\times\frac{\mathrm{Im}D^{\mathrm{R}}_{\bm{q}}(\omega)}{\omega+\varepsilon^{\prime}-\varepsilon-\mathrm{i}\eta}\left[\tanh\left(\frac{\varepsilon^{\prime}}{2T}\right)+\coth\left(\frac{\omega}{2T}\right)\right],$
(S70)

where $k_{1}=\mathrm{min}\\{q_{\mathrm{C}},2k_{\mathrm{F}b}\\}$, and
$q_{\mathrm{C}}$ is a wavevector cutoff of bosons. We can approximate
$\xi_{|\bm{k}|-q}$ and $\xi_{|\bm{k}|+q}$ by $-\infty$ and $+\infty$,
respectively [62]. This procedure is checked for
$\varepsilon_{\mathrm{Fb}}/\omega_{\rm C}\gg 1$ by comparing the result with
the numerical calculation without using the extrapolation, regardless of the
cutoff wavevector $k_{1}$. Moreover, we use the relation
$\mathrm{Re}D^{\mathrm{R}}_{\bm{q}}(\omega)=\omega_{0,\bm{q}}^{2}/(\omega^{2}-\omega_{\bm{q}}^{2})$.
Then, the real part of $\Sigma_{1,\bm{k}}^{\mathrm{R}}(\varepsilon)$ is given
as follows:

$\displaystyle\mathrm{Re}\Sigma^{\mathrm{R}}_{1,\bm{k}}(\varepsilon)$
$\displaystyle=-\frac{\overline{\Gamma}_{3,\bm{k}}\zeta}{\omega^{*}}\varepsilon$
(S71)

for $\varepsilon\ll\omega_{\mathrm{C}}$, where we have introduced the
quantities

$\displaystyle\omega^{*}$
$\displaystyle=\frac{\displaystyle\int_{0}^{k_{1}}\mathrm{d}q\,q\frac{\omega_{0,\bm{q}}^{2}}{\omega_{\bm{q}}^{2}}}{\displaystyle\int_{0}^{k_{1}}\mathrm{d}q\,q\frac{\omega_{0,\bm{q}}^{2}\gamma_{\bm{q}}}{\omega_{\bm{q}}^{4}}},$
(S72) $\displaystyle\frac{\zeta}{\omega^{*}}$
$\displaystyle=\frac{m_{b}}{4\pi^{2}k_{\mathrm{F}b}}\int_{0}^{k_{1}}\mathrm{d}q\,q\frac{\omega_{0,\bm{q}}^{2}}{\omega_{\bm{q}}^{2}}$
(S73)

to make the notations simple. With
$\mathrm{Im}D^{\mathrm{R}}_{\bm{q}}(\omega)\simeq-2\omega_{0,\bm{q}}^{2}\gamma_{\bm{q}}\omega/\omega_{\bm{q}}^{4}$
for $\omega\ll\gamma_{\bm{q}}$, the imaginary part of
$\Sigma_{1,\bm{k}}^{\mathrm{R}}(\varepsilon)$ is obtained as

$\displaystyle\mathrm{Im}\Sigma^{\mathrm{R}}_{1,\bm{k}}(\varepsilon)$
$\displaystyle=-\frac{\overline{\Gamma}_{3,\bm{k}}\zeta}{\omega^{*2}}(\pi^{2}T^{2}+\varepsilon^{2}).$
(S74)

The other contribution can also be evaluated with a similar procedure. We list
the result below:

$\displaystyle\Sigma^{\mathrm{R}}_{2,\bm{k}}(\varepsilon)=-\frac{2\overline{\Gamma}_{4,\bm{k}}\zeta}{\omega^{*}}\varepsilon-\mathrm{i}\frac{2\overline{\Gamma}_{4,\bm{k}}\zeta}{\omega^{*2}}(\pi^{2}T^{2}+\varepsilon^{2}),$
(S75)

where
$\overline{\Gamma}_{4,\bm{k}}=\langle|g_{2}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})|^{2}\rangle_{\bm{k}^{\prime}}\in\mathbb{R}$.

### 3\. Bogolon anomalous self-energy

In this subsection, we derive the bogolon anomalous self-energy. The diagram
of the second-order self-energy is shown in Fig.2 (c) of the main text. With
inversion symmetry, the corresponding self-energy
$S_{\bm{k}}^{\dagger}(\mathrm{i}\varepsilon_{n})$ is given as follows:

$\displaystyle
S^{\dagger}_{\bm{k}}(\mathrm{i}\varepsilon_{n})=-\frac{2T}{(2\pi)^{3}}\sum_{m}\int\mathrm{d}\bm{k}^{\prime}\,g_{1}^{\ast}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})g_{2}^{\ast}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})$
$\displaystyle\quad\times
G^{0}_{\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{m})[D_{\bm{k}-\bm{k}^{\prime}}(-\mathrm{i}\varepsilon_{n}+\mathrm{i}\varepsilon_{m})-D_{\bm{k}-\bm{k}^{\prime}}(\mathrm{i}\varepsilon_{n}+\mathrm{i}\varepsilon_{m})].$
(S76)

The explicit form of the retarded self-energy
$S_{\bm{k}}^{\dagger\mathrm{R}}(\varepsilon)$ can be derived by a similar
manner to the previous subsection. The result is written as

$\displaystyle
S_{\bm{k}}^{\dagger\mathrm{R}}(\varepsilon)=-\frac{4\overline{\Gamma}_{5,\bm{k}}\zeta}{\omega^{*}}\varepsilon-\mathrm{i}\frac{4\overline{\Gamma}_{5,\bm{k}}\zeta}{\omega^{*2}}(\pi^{2}T^{2}+\varepsilon^{2}),$
(S77)

where $\overline{\Gamma}_{5,\bm{k}}=\langle
g_{1}^{\ast}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})g_{2}^{\ast}(\bm{k}^{\prime},\bm{k}-\bm{k}^{\prime})\rangle_{\bm{k}^{\prime}}$.

### 4\. Analytic continuation

For the characterization of the frequency-dependent functional forms, it is
suitable to see the physical quantities on the imaginary (or Matsubara) axis.
The results on the retarded bogolon self-energies obtained in the previous
subsections can be analytically continued to the imaginary axis as

$\displaystyle\Sigma_{\bm{k}}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=a_{\bm{k}}\mathrm{i}\varepsilon_{n}+\mathrm{i}b_{\bm{k}}\left[\pi^{2}T^{2}+(\mathrm{i}\varepsilon_{n})^{2}\right]\mathrm{sgn}\,\varepsilon_{n},$
(S78) $\displaystyle S_{\bm{k}}^{\dagger}(\mathrm{i}\varepsilon_{n})$
$\displaystyle=c_{\bm{k}}\mathrm{i}\varepsilon_{n}+\mathrm{i}d_{\bm{k}}\left[\pi^{2}T^{2}+(\mathrm{i}\varepsilon_{n})^{2}\right]\mathrm{sgn}\,\varepsilon_{n},$
(S79)

where

$\displaystyle a_{\bm{k}}$
$\displaystyle=-(\overline{\Gamma}_{3,\bm{k}}+4\overline{\Gamma}_{4,\bm{k}})\frac{\zeta}{\omega^{*}},$
(S80) $\displaystyle b_{\bm{k}}$
$\displaystyle=\frac{a_{\bm{k}}}{\omega^{\ast}},$ (S81) $\displaystyle
c_{\bm{k}}$
$\displaystyle=-4\overline{\Gamma}_{5,\bm{k}}\frac{\zeta}{\omega^{*}},$ (S82)
$\displaystyle d_{\bm{k}}$ $\displaystyle=\frac{c_{\bm{k}}}{\omega^{\ast}},$
(S83)

from which we can construct both the retarded and advanced Green functions.
Note that $a_{\bm{k}}$ and $b_{\bm{k}}$ are real, while $c_{\bm{k}}$ and
$d_{\bm{k}}$ are complex.

The electron-electron interaction is also expected to lead to a similar effect
discussed above. The imaginary part should be calculated in a manner similar
to Ref. [62], and then the real part may be reconstructed through the Kramers-
Kronig relation [67, 68].

## SM C. Connection to $j=3/2$ fermion model

### 1\. Electronic Hamiltonian and bogolons

We consider the $j=3/2$ fermion model as the simplest description for the
inversion symmetric Bogoliubov Fermi surfaces [5, 20]. The Hamiltonian is
given by

$\displaystyle\mathscr{H}=\sum_{\bm{k}}\vec{c}_{\bm{k}}^{\dagger}\left[\Big{(}\frac{\bm{k}^{2}}{2m_{e}}-\varepsilon_{{\rm
F}e}\Big{)}\hat{1}+\beta(\bm{k}\cdot\hat{\bm{J}})^{2}\right]\vec{c}_{\bm{k}}$
$\displaystyle+\sum_{\bm{k}}\vec{c}_{\bm{k}}^{\dagger}\left[\Delta_{1}k_{z}(k_{x}+\mathrm{i}k_{y})\hat{E}+\frac{2\Delta_{0}}{\sqrt{3}}\lceil\hat{J}_{z}(\hat{J}_{x}+\mathrm{i}\hat{J}_{y})\rfloor\hat{E}\right]\vec{c}_{-\bm{k}}^{\dagger\rm
T}$ $\displaystyle+{\rm H.c.},$ (S84)

where
$\vec{c}_{\bm{k}}=(c_{\bm{k},3/2},c_{\bm{k},1/2},c_{\bm{k},-1/2},c_{\bm{k},-3/2})^{\rm
T}$ is the spin $3/2$ spinor of electrons. The vector $\hat{\bm{J}}$, which is
$4\times 4$ matrix, represents a spin-3/2 operator (or dipole), and $\hat{E}$
is the antisymmetric matrix defined in Ref. [20]. The square bracket
symmetrize the product of matrices as
$\lceil\hat{A}\hat{B}\rfloor=(\hat{A}\hat{B}+\hat{B}\hat{A})/2$.

The time-reversal-symmetry broken superconducting state with inversion
symmetry generally realizes the Bogoliubov-Fermi surfaces which are
topologically protected [5]. Then the resultant effective low-energy
Hamiltonian in Eq. (1) of the main text is derived and the bogolon operators
are given by

$\displaystyle\alpha_{\bm{k}}$
$\displaystyle=\sum_{m}\big{(}u_{\bm{k},m}c_{\bm{k},m}+v_{\bm{k},m}c^{\dagger}_{-\bm{k},m}\big{)},$
(S85)

where the coefficients $u$ and $v$ are obtained by diagonalizing the
Hamiltonian (S84). The odd-frequency pair amplitude of bogolons is now defined
by

$\displaystyle
F_{\bm{k}}(\tau)=-\langle\mathcal{T}\alpha_{-\bm{k}}(\tau)\alpha_{\bm{k}}\rangle,$
(S86)

which is an odd function with respect to time and is induced from the disorder
and correlation effects.

Rank | $\eta$
---|---
$r=0$ (monopole/singlet) | 1
$r=1$ (dipole/triplet) | $x,y,z$
$r=2$ (quadrupole/quintet) | $xy,yz,zx,z^{2},x^{2}-y^{2}$
$r=3$ (octupole/septet) |  | $xyz,xz^{2},yz^{2},z^{3},z(x^{2}-y^{2}),$
---
$x(x^{2}-3y^{2}),y(3x^{2}-y^{2})$
Table S1:  List of the types of diagonal and offdiagonal physical quantities
classified by the rank $r$. See Ref. [20] for the detailed forms of the
matrices $\hat{O}^{\eta}$.

Now we consider the time-dependent order parameters in terms of the original
electrons as

$\displaystyle M^{\eta}_{\bm{k}}(\tau)$
$\displaystyle=\langle\mathcal{T}\vec{c}_{\bm{k}}^{\dagger}\hat{O}^{\eta}\vec{c}_{\bm{k}}(\tau)\rangle,$
(S87) $\displaystyle P^{\eta}_{\bm{k}}(\tau)$
$\displaystyle=\langle\mathcal{T}\vec{c}_{\bm{k}}^{\dagger}\hat{O}^{\eta}\hat{E}\vec{c}^{\dagger\rm
T}_{-\bm{k}}(\tau)\rangle,$ (S88) $\displaystyle
P^{\dagger\eta}_{\bm{k}}(\tau)$
$\displaystyle=\langle\mathcal{T}\vec{c}_{-\bm{k}}^{\rm
T}(\hat{O}^{\eta}\hat{E})^{\dagger}\vec{c}_{\bm{k}}(\tau)\rangle.$ (S89)

The $4\times 4$ matrices $\hat{O}^{\eta}$ are defined in Ref. [20], where all
the order parameters are exhausted by this expression. The pair amplitude
corresponds to the multiplet pairs, i.e., electron ($P$) and hole
($P^{\dagger}$) pair amplitudes. The index $\eta$ represents an identifier for
the multipoles and multiplet pairs, and is classified by the rank $r$ which is
defined up to $2j$. (Note that the symbol $\eta$ in this paper corresponds to
$\eta^{\prime}$ in Ref. [20]). For the diagonal quantity, $r=0,1,2,3$
corresponds to one monopole, three dipoles, five quadrupole, and seven
octupole. For the off-diagonal quantity, $r=0,1,2,3$ correspond to singlet,
triplet, quintet, and septet pairs. We list possible indices in Tab. S1.

The time-dependent multipole functions are represented in terms of bogolon’s
physical quantities. In order to see the contributions from the odd-frequency
pairing of bogolons near the Fermi level, we define the odd-frequency
multipoles $\tilde{M}^{\eta}$, $\tilde{P}^{\eta}$, $\tilde{P}^{\dagger\eta}$
which are induced solely by the degrees of freedom near the Bogoliubov Fermi
surfaces and are written in the forms

$\displaystyle\tilde{M}^{\eta}_{\bm{k}}(\tau)$
$\displaystyle=C^{\eta}_{M}(\bm{k})F_{\bm{k}}(\tau),$ (S90)
$\displaystyle\tilde{P}^{\eta}_{\bm{k}}(\tau)$
$\displaystyle=C^{\eta}_{P}(\bm{k})F_{\bm{k}}(\tau),$ (S91)
$\displaystyle\tilde{P}^{\dagger\eta}_{\bm{k}}(\tau)$
$\displaystyle=C^{\eta}_{P^{\dagger}}(\bm{k})F_{\bm{k}}(\tau).$ (S92)

The quantity $C$ is regarded as a kind of ‘susceptibility’, showing how much
of odd-frequency multipoles and multiplet pair of electrons are induced from
the odd-frequency pair amplitude of bogolons. We then define the quantities
depending only on the rank

$\displaystyle C^{r}_{X}(\bm{k})$
$\displaystyle=\sqrt{\displaystyle\sum_{\eta\in r}|C_{X}^{\eta}(\bm{k})|^{2}}$
(S93)

for $X=M,P,P^{\dagger}$. The result is shown in Fig. 4 of the main text with a
specific choice of the parameters.

### 2\. Impurity potential

The impurity scattering term is given in the language of the original
electrons by

$\displaystyle\mathscr{H}_{\rm imp}$
$\displaystyle=\sum_{i}\int\mathrm{d}\bm{r}\sum_{\eta}U_{\eta}(\bm{r}-\bm{R}_{i})\vec{c}^{\dagger}(\bm{r})\hat{O}^{\eta}\vec{c}(\bm{r}),$
(S94) $\displaystyle\vec{c}(\bm{r})$
$\displaystyle=\frac{1}{\sqrt{V}}\sum_{\bm{k}}\vec{c}_{\bm{k}}\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}},$
(S95)

where we consider the isotropic ($\eta=1$) and anisotropic
($\eta=xy,yz,zx,z^{2},x^{2}-y^{2}$) scattering centers located at
$\bm{R}_{i}$, both of which are electric degrees of freedom compatible with
nonmagnetic impurities. The full list of $4\times 4$ matrices are defined by
using the $\hat{\bm{J}}$ matrix in Ref. [20]. The impurity position
$\bm{R}_{i}$ is to be averaged.

The impurity potential can be rewritten in terms of bogolon’s operators by
using the relation

$\displaystyle c_{\bm{k}m}$
$\displaystyle=u_{\bm{k}m}^{*}\alpha_{\bm{k}}+v_{-\bm{k},m}\alpha_{-\bm{k}}^{\dagger},$
(S96)

where only the fermions that have Fermi surfaces are considered in the right-
hand side. We then obtain the impurity potential term for bogolons as

$\displaystyle\mathscr{H}_{\rm imp}$
$\displaystyle=\frac{1}{V}\sum_{\bm{k},\bm{q}}\rho_{\bm{q}}u_{1}(\bm{k},\bm{q})\alpha_{\bm{k}+\bm{q}}^{\dagger}\alpha_{\bm{k}}$
$\displaystyle\quad+\frac{1}{V}\sum_{\bm{k},\bm{q}}\rho_{\bm{q}}u_{2}(\bm{k},\bm{q})\alpha_{\bm{k}+\bm{q}}^{\dagger}\alpha_{-\bm{k}}^{\dagger}+\mathrm{H.c.}$
$\displaystyle\quad+\mathrm{Const.},$ (S97)

where

$\displaystyle u_{1}(\bm{k},\bm{q})$
$\displaystyle=\sum_{\eta}\sum_{m,m^{\prime}}U_{\eta}(\bm{q})[u_{\bm{k}+\bm{q},m}O^{\eta}_{mm^{\prime}}u_{\bm{k},m^{\prime}}^{\ast}$
$\displaystyle\qquad-
v_{\bm{k},m}^{\ast}O_{mm^{\prime}}^{\eta}v_{\bm{k}+\bm{q},m^{\prime}}],$ (S98)
$\displaystyle u_{2}(\bm{k},\bm{q})$
$\displaystyle=\sum_{\eta}\sum_{m,m^{\prime}}U_{\eta}(\bm{q})u_{\bm{k}+\bm{q},m}O^{\eta}_{mm^{\prime}}v_{-\bm{k},m^{\prime}}.$
(S99)

### 3\. Single-particle spectral function

We connect the single particle spectrum of the electrons to that of the
bogolons. We define the single-particle spectrum for $j=3/2$ electrons by

$\displaystyle A^{\rm elec}_{\bm{k}}(\varepsilon)$
$\displaystyle=\frac{1}{2\pi
f(\varepsilon)}\int_{-\infty}^{\infty}\mathrm{d}t\langle\vec{\psi}_{\bm{k}}^{\dagger}\vec{\psi}_{\bm{k}}(t)\rangle\mathrm{e}^{\mathrm{i}\varepsilon
t},$ (S100)

where we have introduced the Nambu spinor
$\vec{\psi}_{\bm{k}}=(\vec{c}_{\bm{k}}^{\rm
T},\vec{c}_{-\bm{k}}^{\dagger})^{\rm T}$. This is written in terms of
bogolons, and near the Fermi level it has the form

$\displaystyle A^{\rm elec}_{\bm{k}}(\varepsilon)$
$\displaystyle=\frac{1}{2\pi
f(\varepsilon)}\int_{-\infty}^{\infty}\mathrm{d}t\big{[}\langle\alpha_{\bm{k}}^{\dagger}\alpha_{\bm{k}}(t)\rangle+\langle\alpha_{-\bm{k}}\alpha_{-\bm{k}}^{\dagger}(t)\rangle\big{]}\mathrm{e}^{\mathrm{i}\varepsilon
t}$
$\displaystyle=-\frac{1}{\pi}\mathrm{Im}\,\mathrm{Tr}\,\hat{G}_{\bm{k}}^{\rm
R}(\varepsilon),$ (S101)

where the $2\times 2$ matrix $\hat{G}_{\bm{k}}^{\rm R}(\varepsilon)$ is the
retarded version of the Green function defined in Eq. (S9).

References
See the list of references in the main text.