Efficient slave-boson approach for multiorbital two-particle response functions and superconductivity

We develop an efficient approach for computing two-particle response functions and interaction vertices for multiorbital strongly correlated systems based on fluctuation around rotationally-invariant slave-boson saddle-point. The method is applied to the degenerate three-orbital Hubbard-Kanamori model for investigating the origin of the s-wave orbital antisymmetric spin-triplet superconductivity in the Hund's metal regime, previously found in the dynamical mean-field theory studies. By computing the pairing interaction considering the particle-particle and the particle-hole scattering channels, we identify the mechanism leading to the pairing instability around Hund's metal crossover arises from the particle-particle channel, containing the local electron pair fluctuation between different particle-number sectors of the atomic Hilbert space. On the other hand, the particle-hole spin fluctuations induce the s-wave pairing instability before entering the Hund's regime. Our approach paves the way for investigating the pairing mechanism in realistic correlated materials.

Recently, RISB has been reformulated as a quantum embedding theory, where the interacting lattice problem is mapped to an impurity problem coupled to a self-consistently determined environment [21], similar to DMFT and density matrix embedding theory (DMET) [8,23,24]. In particular, the RISB saddle-point equations are equivalent to the "non-interacting bath" DMET (NIB-DMET) self-consistent equations when setting the quasiparticle renormalization matrix to unity and enforcing an additional constraint on the structure of the physical density matrix [25,26]. In addition, the two methods, originally proposed for describing the ground state or low-temperature properties, have been extended to study the finite-temperature effects, the nonequilibrium dynamics, the excited states, and the singleparticle spectral functions in correlated systems [27][28][29][30][31][32][33].
So far, RISB is mostly used for investigating the singleparticle spectral functions and the static local observ-ables. However, the two-particle response functions and the corresponding interaction vertices are also important for explaining the emergent phenomena in correlated materials, e.g., the spin-fluctuation mediated pairing in unconventional superconductors [34]. It is, therefore, important to extend RISB to study these quantities. Indeed, it is possible to compute the two-particle response functions with the Gaussian fluctuation approach around the slave-boson saddle-point [3,4,[35][36][37][38][39][40][41][42][43][44]. However, the technique is so far restricted to the single-orbital Hubbard model. On the other hand, the development of the time-dependent Gutzwiller approximation has been extended to multiorbital systems and applied to the twoorbital Hubbard model for spin susceptibilities [45][46][47][48][49][50][51]. To the best of our knowledge, the theories have not been generalized to compute arbitrary two-particle response functions and quasiparticle interaction vertices for multiorbital systems.
In this work, we develop an efficient approach to compute general susceptibilities and quasiparticle interaction vertices based on fluctuation around the RISB saddlepoint, allowing a diagrammatic analysis for the pairing mechanism. We apply our method to the degenerate three-orbital Hubbard-Kanamori model to investigate the origin of the s-wave orbital-antisymmetric spintriplet pairing instability in the Hund's metal regime, previously found in the DMFT and GA studies [52][53][54][55]. We show that, in agreement with DMFT [53], our approach captures the s-wave spin-triplet pairing instability around the Hund's metal crossover. By investigating the pairing interaction considering the particleparticle and the particle-hole scattering channel, we identify that the mechanism leading to the local s-wave orbital-antisymmetric spin-triplet pairing arises from the particle-particle channel, containing the local electron pair fluctuation between different particle-number sectors of the local Hilbert space. Interestingly, the particlehole spin-fluctuation mechanism for the s-wave pairing, considered also in previous works [53,56,57], induces the s-wave pairing instability slightly before entering the Hund's regime. Possible applications of our formalism to NIB-DMET are also discussed.

II. MODEL
We consider the following generic multi-orbital Hubbard-Kanamori model: where α is the orbital index, σ is the spin index, i is the unit-cell label, and k is the momentum conjugate to i. As a proof of principle and for pedagogical reason, we will assume a three-orbital degenerate model with the energy dispersion of a two-dimensional square lattice with the nearest neighbor hopping: where α ∈ {1, 2, 3}, and we will set t = 1 as the energy unit. However, we note that our formalism applies to multiorbital Hubbard models with general hopping matrix and arbitrary number of orbitals. The term H loc represents the following operator: which contains the Kanamori interaction [58] in the cubic-harmonic basis. The first term is the intra-orbital Coulomb interaction, the second term and the third term is the inter-orbital Coulomb interaction, and the last term contains the spin-flip and the pairing hopping interaction. Throughout our paper, we assume the rotationally invariant condition U = U − 2J and set J = U/4. Note that, with this choice of parameters, the bare orbital-antisymmetric spin-triplet pairing interaction is repulsive, i.e., U − J > 0. The electron occupancy is controlled by the chemical potential µ 0 . Due to the O(3) ⊗ SU (2) symmetry in the degenerate three-orbital model, the orbital-antisymmetric spintriplet pairing channels [59][60][61] are related to each other by a rotation in the orbital and the spin space. Consequently, we focus on the pairing fluctuation in one of the orbital-antisymmetric spin-triplet pairing channels: Similarly, we have the following independent operators for the charge, spin, orbital, and spin-orbital fluctuation channels: where we label the fluctuation channels by s ∈ {ch, sp, orb, so, orb*, so*, P} throughout the paper.
Here, λ 0 is the 3 × 3 identity matrix and λ i are the Gell-Mann matrices (see Appx. A), while σ 0 is the 2 × 2 identity matrix and σ i (i = x, y, z) are the Pauli matrices.

III. METHOD
Our fluctuation approach around the RISB normalstate saddle-point is entirely encoded in the following Lagrange function [62] (see Appx. (B)): (6) where: Equation (7) encodes the contribution of the so-called "quasiparticle fermionic" degrees of freedom. Specifically, the matrix: with the hopping term in the Nambu basis characterizes the "quasiparticle Hamiltonian": where is a Nambu spinor, f kaσ are the fermionic quasiparticle modes, and M is the total number of orbitals. The matrix R is the so-called "quasiparticle renormalization matrix" and Λ is a matrix of Lagrange multipliers enforcing the RISB constraints [9,62]: where ∆ i corresponds to the local quasiparticle density matrices [14], and the symbol ... T denotes the thermal average of the non-interacting quasiparticle Hamiltonian H qp at temperature T . The second term L emb (Eq. (8)) encodes the contribution of the slave-boson amplitudes, that here we expressed directly in terms of the corresponding "quantum embedding" states |Φ i and interacting embedding Hamiltonians [21] (see Appx. B 1): is the Nambu spinor for the bath orbitals. The matrix is the sign exchange matrix generated from the embedding mapping (see Appx. B 1), where 1 is the 2M × 2M identity matrix. The variable E c i is a Lagrange multiplier enforcing the normalization of |Φ i : The matrix Λ c i , describing the embedding Hamiltonian bath potential, is a matrix of Lagrange multipliers enforcing the RISB constraints: The matrix D i , describing the hybridization between the impurity and the bath orbitals, is a matrix of Lagrange multipliers, enforcing the definition of the renormalization matrix [9,14,62] The third term L mix (Eq. (9)) contains the Lagrange multipliers from both L qp and L emb . All physical observables can be obtained from the above variational variables at the saddle-point solution of Eq. (6). The total energy is equal to the Lagrange function (Eq. (6)) evaluated at the saddle-point. The expectation value of generic local operatorsÔ In particular, the local (physical) single-particle density matrix is obtained from: The quasiparticle weight is determined from the R matrix through Z i = R † i R i . Note that within the context of NIB-DMET, Eq. (12) corresponds to the so-called "low-level mean-field" Hamiltonian when setting R = I, and Λ is termed "correlation potential". Equation (14) corresponds to the so-called "high-level many-body Hamiltonian" in NIB-DMET, where the two-particle interaction on the bath orbitals is set to zero [24].
A. Parameterization of the single-particle matrices To enforce the symmetry conditions of the Lagrange function, we introduce the following parameterization of the renormalization matrix, R i , and the Lagrange multipliers, Λ i , ∆ i , D i , and Λ c i [14]: where 1 is the 4M × 4M identity matrix, and h s and h s are the symmetry-adapted matrix basis of the above single-particle matrices. The structure of the matrix basis h s andh s is determined from the group symmetry analysis of the model in the presence of the fluctuating operators (e.g., Eqs. (4)-(5)) [14]. This parameterization allows us to classify the fluctuations of the variational parameters (r s , l s , etc.) to a specific symmetry channel s, associated to h s andh s . For example, in the degenerate three-orbital Hubbard-Kanamori model, the h s andh s (see Appx. C 2) are associated to the fluctuation channels s ∈ {ch, sp, orb, so, orb*, so*, P} in Eqs. (4)- (5). In addition, for computing the susceptibility of a given channel s, the embedding wavefunction |Φ i has to break the corresponding symmetry, e.g., the particle number conservation of |Φ i has to be broken for the pairing susceptibility calculations.
For later convenience, we introduce the following vector of parameters: and assume that all of its entries are real, which is sufficient for static quantities (e.g., static susceptibilities and Landau parameters [39,42]). Note that our assumption of real variables is applicable for our model without spinorbit coupling. The generalization to spin-orbit coupled systems can be straightforwardly obtained using the same procedure proposed in this work, by including in the Lagrangian also the imaginary part of R and D.

B. Saddle-point approximation
The first step of our fluctuation approach is to determine the normal-state saddle-point solution without any ordering. We assume a spatially homogeneous saddlepoint solution, where x i does not depend on i.
Performing the partial derivatives of Eq. (48) with respect to x, we arrive the following saddle-point equations: where f T is the Fermi function and H qp k = R˜ k R † + Λ is the saddle-point quasiparticle Hamiltonian. Equations (27)-(32) can be solved numerically utilizing quasi-Newton methods [14,21]. Note that our saddle-point equations yield consistent results compared to the formalism in Ref. [62].
It is also interesting to point out that Eqs. (27)-(32) are equivalent to the NIB-DMET self-consistent equations when setting the renormalization matrix to unity R = I and enforcing the so-called "quasiparticle constraint" that we will introduce later in Sec. IV [25].
Given the saddle-point solution in the normal phase, we want to compute the corresponding susceptibilities. This will be accomplished using the approach described below.

C. Calculation of susceptibilities
Here we describe the formalism for calculating the susceptibilities in multi-orbital systems within the RISB framework. For concreteness, we focus on uniform susceptibilities in this section, where x i is independent of i and we suppress the i index in the following derivation. The generalization to susceptibilities with finite momentum transfer is described in Sec. IV.
Let us consider the RISB Lagrange function (Eq. (6)) in the presence of a local perturbation, proportional to a generic operatorÔ: where we have modified the embedding part of the Lagrangian to which was obtained by adding a field ξ coupled toÔ in the embedding Hamiltonian of Eq. (14) and expressing the variational parameters in terms of the vector x, see Eq. (26).
To calculate the linear response of the system to the perturbationÔ, we need to evaluate how the saddle-point variational parameters x of Eq. (33) evolves as a function of ξ. For this purpose, it is convenient to introduce the following functional: where |Φ(ξ, x) and E c (ξ, x) are the ground state ofĤ emb and its eigenvalue, respectively, see Eq. (30). Within these definitions, the saddle-point solution of x for a given ξ, that we call x(ξ), is defined by: and the linear response for the operatorÔ is given by the following equation (see Appx. D for derivation): where we introduced the susceptibilities: The so-called "fluctuation matrix" is: Here, the indices µ and ν run through all the variational variables in Eq. (26), i.e., r s , l s , d s , D s , l c s . To keep track of the structure of the fluctuation matrix (where different second order derivatives are computed through different equations, see Appx. E), from now on we will often use these variational variables as matrix subscripts. For example, M Ds,l c s corresponds to the second order derivatives with respect to D s and l c s (see Eq. (E21)). It is important to note that M is not invertible. The reason is that the functional Ω is invariant with respect to the gauge transformation (Eq. (H6)), so M is not unique because of the would-be Goldstone modes. As explained in Appx. I, this redundancy can be systematically resolved by operating a gauge fixing process that removes from the onset of the would-be Goldstone modes [47]. A simpler alternative is to solve the overdetermined linear system (Eq. (40)) by introducing the Moore-Penrose pseudo-inverse of the fluctuation matrix, which we are going to indicate asM −1 . In terms of the pseudo-inverse, the susceptibility can be formally expressed as follows: Note that Eq. (41) applies for general multiorbital Hubbard models, and the procedure for evaluating each element, Eqs. (38), (39), (40), is described in Appx. E. We now discuss the application of our formalism to the degenerate three-orbital Hubbard-Kanamori model. For the considered model, the fluctuation matrix M reduces to a block-diagonal matrix, constructed by seven 5 × 5 matrices shown schematically in Fig. 1 (one for each fluctuation channel s), because of the orthonormality of the fluctuation basis Tr[h s h † s ] = δ ss . Furthermore, for a given channel s, χ emb µO (see Eq. (39)) is nonzero only for the components µ = D s and l s . Therefore, Eq. (41), for a given channel s, can be further simplified to: . We only need to evaluate the 5 × 5 fluctuation matrix and its pseudo-inversion within each s block to compute the corresponding susceptibility. Note that the blockdiagonal structure is not directly applicable to generic systems, because of effects such as orbital differentiation or spin-orbit coupling. In these cases, one has to compute the full fluctuation matrix for calculating response functions.

IV. FERMI-LIQUID APPROXIMATION AND DIAGRAMMATIC APPROACH
The Landau Fermi-liquid theory allows one to describe the thermodynamic properties of metals in terms of an effective non-interacting picture. Importantly, this framework applies only to conserved quantities. In particular, since the superconducting order parameterÔ P does not commute with Eq. (1), the corresponding susceptibility is not rigorously expressible in terms of quasiparticle parameters. Nevertheless, as we are going to show below, within the RISB framework, it is possible to derive an approximate (but accurate) expression for the superconducting susceptibility in terms of the quasiparticle Green's function and interaction vertices. Moreover, the susceptibility can be formulated in terms of the Bethe-Salpeter equation, allowing further diagrammatic analysis for the pairing mechanism.
From the point of view of the RISB methodology, the reason why the superconducting susceptibility cannot be calculated in terms of quasiparticle parameters is that: for s = P, i.e., the physical density matrix is, in general, not the same as the quasiparticle density matrix.
Here we propose to modify the spatially inhomogeneous RISB Lagrange function (Eq. (6)) by imposing the constraint: which is accomplished by introducing additional Lagrange multipliers ζ i,s into Eq. (26) so the x vector be-comes: We also introduce x q , which is the momentum conjugate to x i . The Lagrange function now has the following form: where where we have introduced the physical Green's function: and the quasiparticle Green's function: Similar to the previous section, we also introduced a field ξ k1−k2 coupled to a generic quasiparticle operator This modification will allow us to derive momentum dependent susceptibilities, for investigating the finite momentum (commensurate or incommensurate) instabilities. From now on, we refer to Eq. (44) as the "quasiparticle constraint." Since utilizing the Lagrange equation Eqs. (47)-(49) amounts to solve the RISB equations Eqs. (27)-(32) within a reduced variational space, the corresponding solution is an approximation to the original one. In principle, enforcing the constraint (Eq. (44)) does not affect the results for the conserving channels, where the fluctuating operator commutes with the Hamiltonian, e.g., the charge and spin channels. However, it reduces slightly the variational freedom when the constraint is imposed on the non-conserving channel, e.g., the pairing channel. Nevertheless, as we are going to show, it is always possible to verify a-posteriori the accuracy of the approximation, by comparison to the formalism without the constraint (see also Appx. J).
It is also interesting to point out that Eq. (44) corresponds to the density matrix mapping constraint in DMET [24]. Therefore, the formalism presented in this section is also applicable to the NIB-DMET, by removing the r s sector of the fluctuation basis (Eq. (45)) and setting R = I [25]. This application is discussed in Appx. N.

A. Susceptibility: diagrammatic expression
Here we show how the susceptibility evaluated with the quasiparticle constraint can be expressed in terms of the Feynman diagram in perturbation theory.
Following the procedure in Sec. III C, we introduce the following functional: where now L qp depends on the field ξ q . The linear response for a generic operator is given by the following equation: where the bare susceptibilities are Note again that µ runs through all the elements in Eq. (45), and we use the variational parameters as subscripts. We also introduced the saddle-point Green's function The fluctuation matrix M now depends on momentum q and has an additional component ζ s (see Eq. (45)). The specific form of M is given in Appx. E. Furthermore, M is now an invertible matrix because the quasiparticle constraint breaks the gauge symmetry. Note that Eq. (53) applies for generic multiorbital Hubbard models.
We now discuss the application of our approach to the degenerate three-orbital Hubbard-Kanamori model. As described in the previous section, for the degenerate model considered here, M is a block-diagonal matrix shown schematically in Fig. 1. Also, from Eqs. (4)- (5) and for each fluctuation channel s. Therefore, the susceptibility can be simplified to: where The M −1 rsrs (q) denotes the µ = r s and ν = r s component of M −1 µν (q) and similarly applies to M −1 rsls (q) and M −1 lsls (q). We only need to evaluate the 6 × 6 fluctuation matrix and its inversion within each s block to compute the corresponding susceptibility.
To make a connection to perturbation theory, we compare Eq. (56) with the Bethe-Salpeter representation of the susceptibility: whereΓ s αβγδ (k, k , q) is the (reducible) interaction vertex. To extract theΓ s αβγδ (k, k , q) from Eq. (56), we introduced the following three-leg vertices: such that the susceptibilities can be written as: Substituting Eqs. (62) and (63) into Eq. (56), we obtain the interaction vertex (see Eq. (59)):  64)). The double wavy line corresponds to the dressed bosonic propagator containing the infinite summation of the particle-particle or the particle-hole fermionic bubbles. The black circles denotes the three-leg verticesΛ αβµ (see main text for detail).
describing the effective interaction between quasiparticles mediated by the bosonic propagator M −1 µν in the corresponding channel.
The diagrammatic representation of Eq. (59) is shown in Fig. 2 (a), where the solid line corresponds to the Nambu propagator, the grey circle corresponds toh s , and the grey rectangle corresponds to the interaction vertex Γ s αβγδ (k, k , q). The diagrammatic representation for the interaction vertexΓ s αβγδ (k, k , q) is shown in Fig. 2(b), where the solid circles correspond to the three-leg vertices Λ αβµ . The double wavy line corresponds to M −1 µν (q), which can be viewed as the dressed bosonic propagator (see Appx. F) summing the particle-hole bubbles, for s ∈ {ch, sp, orb, so, orb*, so*}, or the particle-particle bubbles, for s = P , to the infinite order.

B. Landau Fermi-liquid parameters
We can now calculate the Landau Fermi-liquid parameters for the considered three-orbital degenerate model from Eq. (64).
For each channel s ∈ {ch, sp, orb, so, orb*, so*}, we have: where we applied R = R 0 I and Z = R 2 0 for the degenerate model considered here. The scattering amplitude for each particle-hole channel s can be evaluated from where we introduce the Fermi surface average OsOs (0) is the density of state at the Fermi-level, which coincides with the bare susceptibility χ (0) OsOs . The Fermi-liquid parameters F s can be extracted from the scattering amplitude (see Appx. G) From the definition of the quasiparticle susceptibility Eq. (59) and Eq. (68), we obtain the random phase approximation (RPA) like expression for the susceptibilities for s ∈ {ch, sp, orb, so, orb * , so * }. Note that we have applied the Fermi-surface average over k and k . The divergence of the quasiparticle susceptibilities and the scattering amplitudes can be determined from the condition F s (q) = −1. Although Eq. (69) has an RPA-like form, the Fermi-liquid parameters are renormalized by the correlation effect for different q, which provides a more accurate description for strongly correlated systems.
C. Pairing interaction from the particle-particle channel The reducible pairing vertex in the orbitalantisymmetric spin-triplet pairing channel s = P can be computed by projecting the particle-particle scattering vertexΓ P (Eq. (64)) onto the orbital-antisymmetric spin-triplet pairing basis h P (see Appx. (C 2)): where we applied R = R 0 I and Z = R 2 0 for the degenerate model considered here and restrict the pairing at q = 0. The diagrammatic representation for Eq. (70) is shown in Fig. 3(a). In this scattering process, only the particle-particle fermionic bubbles and the local multiplets fluctuation between different particle number sectors in M −1 µν are involved (the fluctuation basis h P and h P in Eqs. (E19)-(E26) selects the fluctuation that does not conserve the particle number. ) We can now derive the RPA-like form for the quasiparticle susceptibility. From Γ sc pp , we compute the reducible pairing interaction by averaging the k and k over the Fermi surface The irreducible pairing interaction Γ irr pp can be extracted from (see Appx. G) From the definition of the quasiparticle susceptibility Eq. (59) and Eq. (72) , we obtain the RPA-like expression for the pairing susceptibility The divergence of the pairing susceptibilities and vertex can be determined from the condition Γ irr pp χ D. Pairing interaction from the particle-hole channel Besides the s-wave pairing induced from the particleparticle vertex, the particle-hole vertices can also induce the local and the non-local pairing through the charge and spin-fluctuation mechanism [63][64][65][66]. To compute the irreducible pairing vertex for the orbital-antisymmetric spin-triplet pairing, we again project the particle-hole vertices onto the pairing basis h P : where the charge, spin, orbital, and spin-orbital scattering verticesΓ s are defined in Eq. (65). The diagrammatic representation for Eq. (74) is shown in Fig. 3(b), where the M −1 rsrs , M −1 rsls , and M −1 lsls contain the summation of the particle-hole bubbles to the infinite order (see Appx. F), and we include both the direct and the exchange (crossing) diagrams. The irreducible pairing interaction from the particle-hole channel can be computed from: where we assume an s-wave pairing to compare with the local pairing fluctuation mechanism in the previous section.

A. Superconducting phase diagram
In this subsection, we apply our RISB saddle-point approximation and fluctuation approach to the degenerate three-orbital Hubbard-Kanamori model with Hund's coupling J = U/4, which serves as an effective model for Hund's metals. We will focus on the order parameter Ô P computed from Eq. (19) and the pairing susceptibility χ P computed from Eq. (42). Figure 4 (a) shows the intensity plot of the spin-triplet pairing order parameter Ô P at T = 0.0005t. The peak of the order parameters locates at the so-called Hund's metal crossover, where the quasiparticle weights Z decrease significantly, as shown in Fig. 4 (c) and (d) for selected fillings n = 1.6, 2.0, 2.4, and 2.8. The faster the decrease in Z, the stronger the enhancement in the pairing order parameters Ô P . The normal state in the superconducting regime can be viewed as Hund's metals, where the quasiparticle weight is small, and the local multiplet is populated with high spin states, favoring the local spin-triplet pairing [16,53,[67][68][69].
We also show the uniform pairing susceptibility χ P evaluated from the fluctuation technique in Fig. 4 (b). The pairing susceptibility is initially positive at small Coulomb interaction U and diverges at the critical point. Then, the pairing susceptibility turns negative, indicating the instability towards the s-wave spin-triplet ordering state. The phase boundary determined from the divergence of the pairing susceptibility is shown in Fig. 4 (a), which agrees with the onset of the mean-field order parameters indicating the consistency of our approach. We also compare our phase diagram with the DMFT results on a Bethe lattice at T = 0.04t rescaled to the 2D bandwidth W = 8t in Fig. 4 (a). While the RISB superconducting regime is broader than the DMFT results, the overall phase diagram agrees qualitatively with the DMFT [53].
We now turn to the finite-temperature phase diagram for the s-wave spin-triplet pairing state. Figure 5 (a) shows the intensity plot of the s-wave spin-triplet order parameters Ô P at U = 8t as a function of electron filling n and temperature T . The superconducting region has a dome shape structure, where the maximum T c locates around n = 2.5. Figure 5 (b) shows the uniform pairing susceptibility χ P computed from the fluctuation approach for filling n = 2.0, 2.4, and 2.8 as a function of temperature T . With decreasing T , the pairing susceptibility increases and diverges at the critical temperature T c . The critical temperature obtained from the divergence of the pairing susceptibility agrees with the onset of the mean-filed order parameters, as shown in Fig. 5 (a). We also compare our phase diagram with the DMFT results on a Bethe lattice in Fig. 5 (a) corresponding to U = 6t rescaled to the 2D bandwidth W = 8t consid- ered here. [53]. Both methods generate a dome shape structure where the peak in RISB is closer to half-filling. Figure 5 (c) shows the intensity plot of the s-wave spin-triplet pairing order parameters Ô P as a function of Coulomb interaction U and J = U/4 at filling n = 2.7. The critical temperature T c peaks around U = 6t, which is around the Hund's metal crossover. Figure 5 (d) shows the corresponding uniform pairing susceptibility χ P computed from the fluctuation approach for U = 5t, 6t, and 12t. The pairing susceptibility diverges at T c and turns negative, indicating the instability towards the s-wave spin-triplet pairing states. The T c obtained from the divergence of the susceptibility again agrees with the onset of the mean-field order parameters, as shown in Fig. 5 (c). We also compare our phase diagram with the DMFT results on a Bethe lattice in Figure 5 (c) at n = 2.0 to match with our critical U c at T = 0.0005t. The phase diagrams obtained from both methods are again similar with a dome shape structure where the T c peaks around the Hund's crossover.
Note that there are two main reasons for expecting qualitative agreement (but quantitative agreement) between our RISB results and the DMFT results of Ref. [53]. The first reason is that RISB (equivalently GA) is essentially a variational approximation to DMFT, in the sense that it is variational in the limit of infinite dimension [70], where DMFT is exact. Also, RISB can be  viewed as an approximation to DMFT, from a quantum embedding perspective, where the uncorrelated bath has the same number of orbitals as the impurity (while the bath is infinite in DMFT). Hence, RISB is expected to be less accurate (but more efficient) compared to DMFT. Nevertheless, we note that, in this work, we assumed a 2D square lattice, while a Bethe lattice was used in Ref. [53]. In fact, it is known that different lattice structures can lead to quantitative differences in the results, but the qualitative behaviors are generally similar [71].

B. Landau parameter and pairing interaction
For studying the pairing mechanism, it is instructive to investigate the quasiparticle interaction vertex in the spin, charge, orbital, spin-orbital, and pairing channel. To obtain these quantities, we applied the Fermi-liquid approximation in Sec. IV, which reproduces the exact physical susceptibility, as shown in Appx. J.
Let us first discuss the charge, spin, orbital, and spinorbital fluctuation, encoded in the Landau parameters F s . The Landau parameters F s in each channel are shown in Fig. 6. We found that the Landau parameters in the charge F ch and orbital F orb (orb*) channels show a peak around the Hund's crossover and diverges at the Mott transition at n = 3. The kink in F ch corresponds to the possible phase separation instability found in the previous slave-spin study [72]. Moreover, we found the instability towards the ferromagnetic ordering F sp = −1 for a wide range of electron filling. Consequently, F sp is the dominant fluctuation in the particle-hole channel. In addition, the spin-orbital channel F so(so*) also shows a subleading instability at n = 3.
We now turn to the irreducible pairing vertex in the particle-particle channel Γ irr pp originated purely from the local pairing fluctuation describing the superconducting instability. Figure 7 (a) shows the behavior of the pairing interaction Γ irr pp in the particle-particle channel as a function of Coulomb interaction U . The condition Γ irr pp χ O P = −1 indicates the divergence in the pairing susceptibility. In the weak-coupling limit, i.e., U t, Γ irr pp follows the bare pairing interaction U − 3J for all the electron filling n. With increasing U , the effective interactions for different electronic filling are renormalized to smaller values and eventually become negative signalizing the instability towards the pairing states. The pairing instability determined from Γ irr pp locates around the Hund's metal crossover as discussed in the previous subsection. On the other hand, as shown in Fig. 7 (b), the pairing instability determined from the particle-hole scattering channel Γ irr ph takes place at a much lower U below the Hund's metal crossover. Consequently, the particle-hole spin-fluctuation mechanism cannot explain the pairing instability around Hund's metal crossover. The strong attraction in Γ irr ph is, however, related to the ferromagnetic instability, as shown in Fig. 6(b).

VI. CONCLUSIONS
On the basis of the fluctuation approach around the RISB normal state saddle-point, we developed an efficient method to compute general susceptibilities, quasiparticle interaction vertex, Fermi-liquid parameters, and pairing interaction for the multiorbital Hubbard model. The method has an RPA-like efficiency and a similar accuracy compared to DMFT for correlated systems.
We applied our method to the degenerate three-orbital Hubbard-Kanamori model to investigate the origin of the s-wave orbital-antisymmetric spin-triplet pairing in Hund's metal, previously found in the DMFT studies [53]. We showed that, in agreement with DMFT, the pairing susceptibility of the s-wave spin-triplet pairing states diverges around the Hund's metal crossover. The phase diagram is in good qualitative agreement with DMFT. By computing the pairing interaction considering the particle-particle and the particle-hole scattering channel, we identified that the origin of the superconducting pairing around Hund's crossover arises from the particle-particle channel, containing the local electron pair fluctuation between different particle-number sectors of the local Hilbert space. The pairing interaction is strongly renormalized in the incoherent Hund's metal regime and becomes negative. On the other hand, the particle-hole spin-fluctuation mechanism induces an s-wave pairing instability already for a smaller value of Coulomb interaction, before entering the Hund's regime.
The local interorbital pairing mechanisms revealed in this work can be applied to the s-wave orbitalantisymmetric spin-triplet pairing states proposed for Sr 2 RuO 4 [59,60,[73][74][75] and KFe 2 As 2 [76,77], where the interplay between the Hund's rule coupling and the spinorbital coupling leads to intriguing gap structures on the Fermi surface. Our approach provides an efficient route for investigating the pairing mechanism for these materials, with the combination of density functional theory. The general formalism that we presented is also applicable for different purposes. For example, it could be utilized for investigating the response functions in the correlation-induced topological materials, e.g., the topological Kondo and topological Mott insulators [22,[78][79][80], and the recently proposed topological iron-based superconductors [81,82]. In addition, the diagrammatic approach proposed in this work may serve as a basis for the non-local extensions beyond RISB, similarly to the diagrammatic approaches beyond DMFT [83]. Finally, our formalism can be applied to the NIB-DMET and other similar quantum embedding methods [84][85][86]. We use the following convention for the Gell-Mann matrices where λ 1 , λ 2 , λ 3 describe the symmetric interorbital interactions or pairings; λ 4 , λ 5 , λ 6 describe the antisymmetric interorbital interactions or pairings; λ 7 , λ 8 , λ 0 describe the intraorbital interactions or pairings. This set of matrices is the most general basis that parameterizes the 3 × 3 quadratic operators in the orbital space for three-orbital models. In the degenerate three-orbital Hubbard-Kanamori model, the O(3) symmetry implies that the order parameters corresponds to the symmetric interorbital fluctuations λ 1 , λ 2 , and λ 3 are identical to each other. Similarly, the order parameters corresponds to the antisymmetric interorbital fluctuations λ 4 , λ 5 , and λ 6 are identical to each other.

Appendix B: Rotationally-invariant slave-boson Nambu formalism
In this section, we outline the basis of the RISB Nambu formalism. We start from a generic multiorbital Hubbard model in the Nambu notation: is the energy dispersion in the Nambu basis.
We also define the Nambu spinor Ξ † where M is the total number of orbitals. The H loc contains the generic local one-body and two-body interactions.
Within RISB framework, the physical operators Ξ iα is mapped to the product of a renormalization matrix and a quasiparticle Nambu spinor: where the quasiparticle spinor is , and the renormalization matrix has the following form [62,87,88]: corresponds to the local quasiparticle Nambu density matrix, and Φ i An is the slave-boson amplitude matrix. We also define the matrices Ξ iα AB = A|Ξ iα |B and Ψ iα nm = n|Ψ iα |m for the fermionic operator in the arbitrary local many-body basis |A and the local Fock basis |n , respectively [9,21]. The local interactions can be expressed in terms of the bosonic amplitudes as [9] H loc = where H loc AB = A|Ĥ loc |B . In order to select the physical states out of the enlarged boson and quasiparticle Hilbert space, one has to enforce the following RISB constraints [9,62] The first constraint limits the Hilbert space to the singleboson states, while the second constraint ensures the rotational invariance of the quasiparticle density matrix under the gauge transformation (see Appx. H).
With the RISB representations and constraints (Eqs (B3)-(B8)), the RISB Lagrangian for the generic Hubbard model (Eq. (B1)) can be expressed as: where the original kinetic hopping term in Eq. (B1) is described by the quasiparticle Hamiltonian: i are the Lagrange multipliers enforcing the RISB constraints (Eqs. (B7) and (B8)) and the structure of the R i matrix (Eq. (B4)). Note that all these single-particle matrices contains the particle, hole, and anomalous sector defined as follows: The Λ i , Λ c i and ∆ i are Hermitian matrices, and the R i and D i are non-Hermitian matrices. These singleparticle matrices are parameterized by Eqs. (21)-(25) utilizing the matrix basis h s andh s , whose structure (for the three-orbital degenerate Hubbard-Kanamori model) is discussed in Sec. C 2.
The slave-boson amplitude can be constructed from the symmetry adaptive basis φ ip : and the matrix basis commutes with all the symmetry operation in the group G of the given problem, i.e., [φ ip , R(g)] = 0 ∀g ∈ G. The procedure for determining φ ip is discussed in Appx. C 1.

Embedding mapping
We now introduce the embedding wavefunction [21] |Φ i = where U PH is the particle-hole transformation on the bath site and N n is the particle number of Fock state |n . Substituting the following identities to Eq. (B9): and 1 is the identity matrix, we obtained the RISB Lagrangian in terms of |Φ i in Eq. (6) in the main text.

Appendix C: Variational basis
In this section, we describe the construction of our variational many-body basis φ p and the single-particle basis h s andh s of our fluctuation approach to the degenerate three-orbital Hubbard-Kanamori model.

Many-body basis
For the charge, spin, orbital, and spin-orbital fluctuations, we construct the many-body basis in Eq. (B16)   using the symmetry adapted basis. The procedure can be found in Ref. [89]. On the other hand, for the pairing state, we construct the many-body variational basis following the procedure in Ref. [62]. First, since the Hubbard-Kanamori interaction (Eq. (3)) can be written into the local Hamiltonian is diagonalized in the Γ = (N, L, S) basis. The σ is a vector of Pauli matrices, and αβγ is the Levi-Civita symbol, which can be expressed in terms of Gell-Mann matrices λ 4 , λ 5 , and λ 6 . Therefore, the slave-boson amplitude can be significantly reduced to where E Γ and |Γ is the eigenvalue and the eigenstate of Eq. (C1), respectively. Comparing Eq. (C5) to Eq. (B16), we identify that the many-body basis for the normal state part is with the corresponding slave-boson c p = Φ(E Γ ), and the pairing part are with the corresponding slave-boson amplitudes c p = Φ(E Γ ; 2q) and c p = Φ(E Γ ; −2q), respectively. In the end, we have 43 bosonic amplitudes listed in Tab. I.

Single-particle basis
The single-particle basis h s andh s , parameterizing Eqs. (B11)-(B15), are block matrices, where the component h s corresponds to the normal part and h s corresponds to the anomalous part of the matrix. The components for each fluctuation channel, in the degenerate three-orbital model, are as follow: for the normal part, and for the anomalous part, where the basis is chosen to be normalized, i.e., Tr h s h † s = 1. We see that h P describes the pairing fluctuation, while h ch , h sp , h orb , h so , h orb * , and h so* describes the charge, spin, orbital, and spinorbital fluctuations, respectively. The linear response for a generic operator is given by the following equation: Note again that µ runs through all the variational variables in x (Eq. (26)), and we use the variational parameters as the subscripts. To evaluate Eq. (D1), it is necessary to calculate dxµ dξ ξ=0 , which can be determined by taking the total derivative of Eq. (36) with respect to ξ, as follows: where M is the fluctuation matrix defined in Eq. The second part M qp , which involves the partial derivatives of the quasiparticle term of the Lagrangian L qp with respect to r s and l s , is computed from the following equations: and the other unlisted components of M qp are zero. We also defined k = (ω n , k) and k ≡ k ωn . Note that since we consider degenerate three-orbital model, at the normal-state saddle-point, the renormalization matrix, the local potential, the quasiparticle energy dispersion, and the Green's functions are all degenerate and diagonal matrices, i.e., where E qp k = R 2 0 k + l 0 and I is the 6 × 6 identity matrix. The Matsubara summation for the fermionic Green's function convolutions in M qp rr , M qp rl , and M qp ll can be evaluated analytically from the Lindhard function. For example, the particle-hole convolution: (E14) and the particle-particle convolution: The analytical continuation to real frequency can be achieved by replacing iΩ n → ω + i0 + . The third part M emb involves the partial derivatives of the embedding term of the Lagrangian L emb with respect to D s , l c s , and ζ s , which can be evaluated as follows. First, we evaluate the first order derivatives using the Hellmann-Feynman theorem: (E18) Then, we can compute the second order derivatives from the following equations: where the other unlisted components of M emb are zero. The above second-order derivatives and Eqs. (38)-(39) can be evaluated using the linear response theory. We apply a perturbation to the embedding Hamiltonian where E n is the n-th excited state energy ofĤ emb and |n is the n-th excited state wavefunction ofĤ emb . Beside the method proposed in Eqs. (E28), one can also use the finite difference method to evaluate the partial derivatives in Eqs. (E19)-(E26). Note that both methods requires the diagonalization of the embedding Hamilto-nianĤ emb , which is the most time-consuming part of the linear-response calculations. With the current state-ofthe-art, we can easily study the f -electron materials, containing 7 correlated orbitals, using exact-diagonalization and machine learning techniques [90]. For the systems with more correlated orbital, one may also utilize the density matrix renormalization group or auxiliary-field quantum Monte Carlo methods [91].

Appendix F: Fluctuation matrix as a bosonic propagator
Here we discuss how the fluctuation matrix can be interpreted as the propagator for the fluctuations of the bosonic variables x i . Let us expand the Lagrangian, Eq. (46), to the second order in δx t i = (δr ch , δl ch , δd ch , δD ch , δl c s , δζ ch , ..., δr s , δl s , δd s , δD s , δl c s , δζ s , ..., δr P , δl P , δd P , δD P , , δl c P , δζ P ) (F1) around the normal-state saddle-point: whereΛ αβµ are the three-leg vertices defined in Eq. (60) and (61) and G(k) is the Nambu propagator. We also introduce the four-leg vertex: (F3) We immediately see that the q independent part of the fluctuation matrix: can be viewed, in the Gaussian fluctuation sense [38], as the inverse of the bare bosonic propagator D −1 0 . It is important to note that D 0 describes the local multiplet fluctuations because it contains the embedding susceptibilities shown in Eqs. (E19)-(E26). We see that, for the pairing channel s = P in Eqs. (E19)-(E26), the multiplet fluctuation selects the basis h P that increases and removes electron pairs from the saddle-point wavefunction. Therefore, it describes the local fluctuation with pair excitations. On the other hand, for channel s ∈ {ch, sp, orb, so, orb*, so*}, the particle number is conserved. Consequently, they describe the corresponding local charge, orbital, and spin fluctuations. G Particle-hole channel: Particle-particle channel: We now discuss the role of M qp (q). By integrating out the fermionic field Ξ kα in Eq. (F2) to the one-loop order, we found the self-energy correction is related to the fluctuation matrix through π(q) ≡ −M qp (q). Therefore, we can write the total fluctuation matrix in terms of the Dyson equation: The total fluctuation matrix corresponds to the dressed bosonic propagator with the self-energy correction summing the fermionic bubbles to the infinite order. From Eq. (E7)-(E9), we see that M qp contains only the particle-particle bubbles for the pairing channel s = P , and the particle-hole bubbles for the other channels s ∈ {ch, sp, orb, so, orb*, so*}. Figure 8 shows the diagrammatic representation of the Dyson equation for the particle-hole and the particle-particle channels.
Appendix G: Random phase approximation for the interaction vertex In this section, we derive the random phase approximation for the interaction vertex at q = 0. Therefore, we suppress the q dependent ofΓ,λ, D, and π in the following derivation. The interaction vertex has the following form (see Eqs. (64) and (F5)): whereΛ µ is the three-leg vertex and D ≡ M −1 is the bosonic Green's function defined in Eq. (F5). We want to obtain an RPA like form for the vertex: (G2) after averaging k and k over the Fermi surface, where F s is the Landau parameter.
For the model considered in this work, where we restricted the variational variables x to real numbers (Eqs. (21)-(25)), the U (1) gauge degrees of freedom in the charge, spin, orbital, and spin-orbital channels are fixed. However, we are left with one gauge degree of freedom relating to the Nambu pseudo-spin rotation generator: where τ i is the Pauli matrix corresponding to Nambu pseudospin. From the definition of the gauge transformation (Eqs. (H1)-(H3)), we derive the kernel K: where the K vector is only non-zero in the pairing channel. We can then construct the vector space v i,µ using the Gram-Schmidt process and compute the susceptibilities through Eq. (I11)-(I13).

Appendix J: Validity of the Fermi-liquid approximation
In this section, we show the pairing susceptibility χ P computed from the equation without enforcing the "quasiparticle constraint" (Eq. (42)) and the equation with the "quasiparticle constraint" (Eq. (56)) in Fig. 9. The χ P obtained from the two approaches are identical for all the parameter regime, indicating the validity of the Fermi-liquid approximation described in Sec. IV. We also compute the bare pairing interaction in the particle-hole channel defined as follows: In this case, the summation of the fermionic particlehole bubbles in Fig. 3 (b) are ignored and only the bare interaction (Landau parameters F s at q = 0) is considered. Figure 10 shows the bare pairing interaction in the particle-hole channel. We found that the bare pairing   pairing state, where the superscript N and sc corresponds to the energy in the normal state and the superconducting state, respectively. The energetic in both Fig. (a) and (b) shows a typical weak-coupling to strongcoupling crossover behavior [71,92,93], where the energy gain is dominated from the potential energy in the weakcoupling limit, and from the kinetic energy in the strongcoupling limit. Interestingly, we find this crossover locates around the Hund's metal crossover where the quasiparticle weight drops rapidly and the superconducting order parameter shows a pronounce peak.
where the fluctuation matrix M is given in Appx. E and we have to enforce R = I in each element. We have also introduced the following susceptibilities: where O is the single-particle matrix representation of a generic operator. The Green's function has the following form: and G ωn,k is the Green's function evaluated at ξ = 0. We also introduce the quasiparticle Hamiltonian (lowlevel mean-field Hamiltonian): where Λ corresponds to the correlation potential in NIB-DMET. For the degenerate model considered in this work, the susceptibility can be written as: χ OsOs (q) = χ Finally, we comment on the advantages and the disadvantages between RISB and NIB-DMET. One advantage of RISB with respect to NIB-DMET is the presence of the renormalization matrix R. It allows the description of the Mott transition within the single-site approach [94], while in the standard NIB-DMET, one has to use at least a two-site cluster to capture the Mott transition [24]. On the other hand, the additional determination of R in RISB may require more self-consistency iterations with respect to NI-DMET, leading to more diagonalization of the embedding HamiltonianĤ emb . Nevertheless, the performance and the accuracy of the two methods are similar [26]. Note that our approach does not apply to the "interacting bath" construction of DMET (IB-DMET), which produces more accurate results than the NIB-DMET [24,85,95]. The extension of our approach to IB-DMET will be an interesting future topic.