Spin-triplet superconductivity from quantum-geometry-induced ferromagnetic fluctuation

Introduction.— Unconventional superconductivity beyond the canonical Bardeen-Cooper-Schrieffer theory shows rich physical phenomena including high-temperature superconductivity and topological superconductivity. Various fluctuations arising from many-body interactions play the main role in the Cooper pairing for unconventional superconductivity, and low-dimensional fluctuations are particularly favorable. For example, it is argued that high-temperature superconductivity in cuprates is mediated by two-dimensional antiferromagnetic fluctuation Moriya and Ueda (2000); Tsuei and Kirtley (2000); Yanase et al. (2003). Also, in iron-based high-temperature superconductors the extended $s$-wave pairing is mediated by orbital Kontani and Onari (2010); Yanagi et al. (2010); Onari and Kontani (2012) or antiferromagnetic Mazin et al. (2008); Kuroki et al. (2008) fluctuation Hosono and Kuroki (2015); Chubukov (2012); Chubukov et al. (2016).

However, searching for topological superconductivity Qi and Zhang (2011); Alicea (2012); Sato and Fujimoto (2016); Sato and Ando (2017) with Majorana fermion Sato (2009, 2010); Fu and Berg (2010) is an unresolved problem of modern condensed matter physics, which is attributed to the fact that the platform for topological superconductivity is rare in nature. Spin-triplet superconductors are canonical candidates, and it is expected that ferromagnetic fluctuation mediates the spin-triplet Cooper pairing. However, candidate materials are restricted to a few heavy-fermion systems with three-dimensional multiple bands Sauls (1994); Tou et al. (1998); Joynt and Taillefer (2002); Saxena et al. (2000); Aoki et al. (2001); Huy et al. (2007); Ran et al. (2019); Aoki et al. (2022).

In the two-dimensional isotropic continuum models, ferromagnetic fluctuation is not favored because of the constant density of states (DOS), which may imply the absence of two-dimensional spin-triplet superconductivity. Even for the anisotropic lattice systems, most quasi-two-dimensional superconductors do not show ferromagnetic fluctuation and antiferromagnetic fluctuations are rather ubiquitous, as we mentioned above for cuprates and iron-based compounds. Thus, spin-triplet superconductivity from ferromagnetic fluctuation is expected to require peculiar band structures, and the search for such systems is challenging for both materials and theoretical models. In this Letter, nevertheless, we propose a guiding principle for realizing ferromagnetic fluctuation in two-dimensional systems by referring to the quantum geometry of Bloch electrons, which is recently attracting much attention in various fields Marzari and Vanderbilt (1997); Resta (2011); Gao et al. (2014); Gao and Xiao (2019); Lapa and Hughes (2019); Daido et al. (2020); Kitamura et al. (2021); Peotta and Törmä (2015); Liang et al. (2017); Törmä et al. (2022); Rossi (2021); Julku et al. (2021a, b); Liao et al. (2021); Chen and Law (2023); Iskin (2023); Rhim et al. (2020); Watanabe and Yanase (2021); Ahn et al. (2020, 2021).

The importance of quantum geometry in superconductors has recently been recognized as it gives correction to the superfluid weight Peotta and Törmä (2015); Liang et al. (2017); Törmä et al. (2022); Rossi (2021) . In the flat-band systems  Julku et al. (2016, 2020); Xie et al. (2020); Hu et al. (2019); Tian et al. (2023); Törmä et al. (2022); Rossi (2021) the superfluid weight from Fermi-liquid theory vanishes, and the quantum geometric contribution determines the superfluid weight. The quantum geometry also plays essential roles in the monolayer FeSe Kitamura et al. (2022a) and some finite-momentum Cooper pairing states Kitamura et al. (2022b); Chen and Huang (2023); Jiang and Barlas (2023); Kitamura et al. (2023); Chazono et al. (2023). However, how quantum geometry affects the pairing mechanism of superconductivity has not been revealed. This work elucidates a way to create a pairing glue of unconventional superconductivity via quantum geometry.

To show that the quantum geometry enables strong ferromagnetic fluctuation in two-dimensional systems, resulting in spin-triplet superconductivity, we elucidate the criterion for ferromagnetic fluctuation in the multi-band system with SU(2) symmetry. We find that the criterion is given by the generalized electric susceptibility (GES) which is defined as a natural extension of the electric susceptibility to metals. The GES contains the terms obtained by the effective mass and the quantum geometry.

The key physics of quantum-geometry-induced ferromagnetic fluctuation, which is shown below, is nontrivial quantum geometry, especially Fubini-Study quantum metric Provost and Vallee (1980); Resta (2011), from non-Kramers band degeneracy. As shown in this Letter, the dispersive Lieb lattice model with non-Kramers band degeneracy shows strong ferromagnetic fluctuation by this mechanism. Solving the linearized gap equation with the effective interaction calculated by the random phase approximation (RPA), spin-triplet superconductivity is demonstrated.

Criterion for ferromagnetic fluctuation in multi-band Hubbard models.—

We consider the multi-band Hubbard model with SU(2) symmetry, which contains multiple degrees of freedom such as orbitals and sublattices Sup . The SU(2) symmetry means that the spin-orbit coupling and the magnetic field are absent. For the interacting Hamiltonian, we consider the onsite Coulomb interaction $U$ strong enough for the superconducting transition, by assuming strongly correlated materials. We then focus on the momentum dependence of the fluctuation, which mainly determines the superconducting symmetry Yanase et al. (2003). While we consider two-dimensional systems, the following discussions apply to three-dimensional systems.

Throughout this paper, $U$ is treated in the RPA scheme. When the system has only one band, the spin (charge) susceptibility $\chi_{\rm s(c)}(\bm{q},i\Omega_{n})$ can be obtained as $\chi_{\rm s(c)}(\bm{q},i\Omega_{n})=\chi_{\rm s(c)}^{0}(\bm{q},i\Omega_{n})/(1\mp\frac{U}{2}\chi_{\rm s(c)}^{0}(\bm{q},i\Omega_{n}))$ by using the bare spin (charge) susceptibility of noninteracting systems, $\chi^{0}_{\rm s(c)}(\bm{q},i\Omega_{n})$. The interaction does not change the position of peaks in the momentum $\bm{q}$ space. Therefore, also for most multi-band systems, it is expected that the momentum dependence of fluctuations arises from the bare susceptibility. Because the low-frequency spin (charge) fluctuation plays the dominant role in mediating superconductivity, hereafter we focus on the static fluctuations at $\Omega_{n}=0$.

In multi-band systems with SU(2) symmetry, the bare spin/charge susceptibilities hold the relationship $\chi^{0}_{\rm s}(\bm{q})=\chi^{0}_{\rm c}(\bm{q})=2\chi^{0}(\bm{q})$ with the bare susceptibility $\chi^{0}(\bm{q})$. Thus, our main concern is the presence/absence of the peak of $\chi^{0}(\bm{q})$ at $\bm{q}=0$, corresponding to the presence/absence of ferromagnetic fluctuation. The structure of susceptibility $\chi^{0}(\bm{q})$ around $\bm{q}=0$ is determined by the curvature $\lim_{\bm{q}\rightarrow 0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$ with $\mu,\nu=x,y$ qch . As a result, the criterion for the ferromagnetic fluctuation is given by the sign of the curvature (see Fig. 1). Ferromagnetic fluctuation may be present when $\lim_{\bm{q}\rightarrow 0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$ is negative. Otherwise, ferromagnetic fluctuation is prohibited.

The curvature $\lim_{\bm{q}\rightarrow 0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$ itself has a physical meaning. For the discussion, it is useful to consider the charge susceptibility in insulators at zero temperature, instead of the spin susceptibility. Based on the Kubo formula, the curvature expresses the correction to the charge density, $\delta\braket{\hat{n}(\bm{r})}$, by the external electric field $E_{\nu}(\bm{r})$ as, $\delta\braket{\hat{n}(\bm{r})}=-\sum_{\mu\nu}\partial_{r_{\nu}}(\lim_{\bm{q}\rightarrow 0}\frac{1}{2}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi_{\rm c}^{0}(\bm{q})E_{\nu}(\bm{r}))$  Sup . This means that the curvature $\lim_{\bm{q}=0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$ is the electric susceptibility. Thus, by generalizing the concept of the electric susceptibility to metals, we define the generalized electric susceptibility (GES) as $\chi_{\rm e}^{0:\mu\nu}\equiv\lim_{\bm{q}=0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$ Sup .

Formula of GES.— Here, we derive the formula of GES Sup , $\chi_{\rm e}^{0:\mu\nu}=\chi_{\rm e:geom}^{0:\mu\nu}+\chi_{\rm e:mass}^{0:\mu\nu}$,

where $\epsilon_{n}(\bm{k})$ is the energy of the noninteracting Hamiltonian $\sigma_{0}\otimes H_{0}(\bm{k})$, which follows $H_{0}(\bm{k})\ket{u_{n}(\bm{k})}=\epsilon_{n}(\bm{k})\ket{u_{n}(\bm{k})}$ with the Bloch wave function $\ket{u_{n}(\bm{k})}$. Note that $\sigma_{0}$ is the unit matrix of spin space and $n$ is the band index. Thus, GES is given by the two terms, $\chi_{\rm e:geom}^{0:\mu\nu}$ and $\chi_{\rm e:mass}^{0:\mu\nu}$.

The first term $\chi_{\rm e:geom}^{0:\mu\nu}$ named quantum geometric term is determined by the geometric quantities, namely, the Fubini-Study quantum metric $g_{n}^{\mu\nu}(\bm{k})=\sum_{m(\neq n)}A_{nm}^{\mu}(\bm{k})A_{mn}^{\nu}(\bm{k})+{\rm c.c.}$ and the positional shift $X_{n}^{\mu\nu}(\bm{k})=\sum_{m(\neq n)}(A_{nm}^{\mu}(\bm{k})A_{mn}^{\nu}(\bm{k})+{\rm c.c.})/(\epsilon_{m}(\bm{k})-\epsilon_{n}(\bm{k}))$ with the Berry connection $A_{nm}^{\mu}(\bm{k})=i\braket{\partial_{k_{\mu}}u_{n}(\bm{k})}{u_{m}(\bm{k})}$. This term arises from purely interband effects and is absent in single-band systems. In this term, the contributions from the quantum metric and the positional shift are competitive. First, the quantum metric Provost and Vallee (1980); Resta (2011), which is the counterpart of the Berry curvature Berry (1984), represents the distance between two adjacent states and is a positive definite tensor. Therefore, combined with negative $f^{\prime}(\epsilon_{n}(\bm{k}))$, the contribution from the quantum metric is always negative, favoring ferromagnetic fluctuation. Second, the positional shift Gao et al. (2014) means the shift of electrons by the external electric field. In insulators at zero temperature, the contribution from the positional shift corresponds to the well-known formula of electric susceptibility Aversa and Sipe (1995). This term can be rewritten as it is proportional to $F_{nm}(\bm{k})(A_{nm}^{\mu}(\bm{k})A_{mn}^{\nu}(\bm{k})+{\rm c.c.})$ with the integrand of the Lindhard function, $F_{nm}(\bm{k},\bm{q})=(f(\epsilon_{m}(\bm{k}))-f(\epsilon_{n}(\bm{k}+\bm{q})))/(\epsilon_{n}(\bm{k}+\bm{q})-\epsilon_{m}(\bm{k}))\xrightarrow{\bm{q}\to 0}F_{nm}(\bm{k})$. Therefore, this contribution is always positive, which favors antiferromagnetic fluctuation.

Importantly, both quantum metric and positional shift diverge at the non-Kramers band-degenerate point. Therefore, quantum geometry plays an essential role when non-Kramers band degeneracy exists. However, the total geometric term does not diverge because of the cancellation of two contributions Sup .

The effective-mass term $\chi_{\rm e:mass}^{0:\mu\nu}$ of GES is the purely intraband effect and is determined by the band dispersion through the effective mass $[m_{n}^{\mu\nu}(\bm{k})]^{-1}=\partial_{k_{\mu}}\partial_{k_{\nu}}\epsilon_{n}(\bm{k})$. In single-band systems, only this term is finite. This term can be positive and negative. For the hyperbolic dispersion $\epsilon_{n}(\bm{k})=k^{2}/2m$, the effective-mass term is zero because the DOS and effective mass are constants, which means the absence of ferromagnetic fluctuation Sup .

GES with non-Kramers band degeneracy.— Because the non-Kramers band degeneracy enhances the quantum geometry, we focus on the Lieb lattice, which has been realized in ultracold atoms allowing us to tune the strength of $U$ Gross and Bloch (2017); Taie et al. (2015), with the experimental test in mind. The Lieb lattice hosts the flat band with three-fold band degeneracy, and the ground state shows the flat-band ferromagnetism Lieb (1989). To distinguish the quantum-geometry-induced ferromagnetic fluctuation from the flat-band ferromagnetism, we study the dispersive Lieb lattice model in which the second and third-nearest-neighbor hoppings are finite. Unlike the usual Lieb lattice with only the nearest-neighbor hopping, the flat band becomes dispersive and the three-fold band degeneracy at the $M$ point [${\bm{k}}=(\pi,\pi)$] is partially lifted, while the two-fold degeneracy remains protected by the $C_{4}$ rotation symmetry Sup .

The dispersive Lieb lattice model is illustrated in Fig. 2(a). The Fermi surfaces for the chemical potential $\mu_{\rm c}=0.5,0.7$, and $0.9$ are shown in Fig. 2(b), and the band dispersion is in Fig. 2(c). The band-degenerate point lies on the Fermi surface, when $\mu_{\rm c}=0.7$. As shown in Fig. 2(d), the maximum of DOS corresponds to $\mu_{\rm c}=0.7$.

In Fig. 3(a), we show the chemical-potential dependence of GES $\chi_{\rm e}^{0:xx}$. In some regions near the Lifshitz transitions ($\mu_{\rm c}\simeq-0.1$ and $0.9$), the GES shows the dip structure. This structure is induced by the effective-mass term. The effective-mass contribution from each band is proportional to an odd function $f^{(2)}(\epsilon_{n}(\bm{k}))$, and therefore, the effective-mass term tends to cancel out between the states below and above the Fermi energy. However, the cancellation is incomplete for $\mu_{\rm c}$ near the Lifshitz transition point, and thus, the effective-mass term gives a negative GES. This is an understanding of why ferromagnetic fluctuation appears at finite temperatures when the Fermi surface is small, from the viewpoint of the GES.

In contrast, accompanying the band degeneracy on the Fermi surface, we obtain the maximally negative value of GES $\chi_{\rm e}^{0:xx}$ at $\mu_{\rm c}=0.7$, which is dominated by the quantum geometric contribution. As expected from the band degeneracy at the M point, the quantum geometric term of the GES mainly comes from the region near the M point. This is verified by the $\bm{k}$-resolved quantum geometric contribution shown in Fig. 3(b). We find a large negative contribution to the GES from the vicinity of the $M$ point, which in turn induces ferromagnetic fluctuation.

As we have mentioned, the quantum metric gives a negative contribution to the GES, while the positional shift positively contributes. Our results imply that the quantum metric overcomes the positional shift when the band-degenerate point lies on the Fermi surface. This can be intuitively understood from the formula of the quantum geometric term. The quantum metric contributes to the GES with $f^{\prime}(\epsilon_{n}(\bm{k}))$, which is divergent on the Fermi surface at low temperatures, [$f^{\prime}(0)\propto 1/T$]. On the other hand, $F_{nm}({\bm{k}})$ in the positional shift contribution is a regular function. Therefore, the quantum metric becomes significant in the presence of band degeneracy at low energies. Consistent with the intuitive explanation, the geometric term is negatively enhanced at low temperatures owing to the contribution of quantum metric, as shown in Fig. 3(c). The geometric term is well fitted by the scaling $\chi_{\rm e:geom}^{0:\mu\nu}=a/T+b$ with constants $a,b$. Thus, we conclude that the quantum metric on the Fermi surface induces ferromagnetic fluctuation when the non-Kramers band degeneracy lies on the Fermi surface.

However, when the band-degenerate point is slightly off the Fermi surface and temperature decreases so that $T\ll|\mu_{\rm c}-0.7|$, the negative geometric term is suppressed as shown in Fig. 3(d). This is consistent with the fact that the quantum metric contribution is a Fermi-surface term. As the band-degenerate point moves away from the Fermi surface by much more than $T$, $f^{\prime}(\epsilon_{n}(\bm{k}))$ near the $M$ point decays, and the quantum metric contribution is suppressed. At low temperatures, the positional shift contribution overcomes the quantum metric contribution, and the quantum geometric term is positive. Thus, in this case, the ferromagnetic-antiferromagnetic crossover of fluctuation occurs as the temperature decreases.

Quantum-geometry-induced ferromagnetic fluctuation.— Then, to justify the above discussion, we show the bare spin susceptibility defined by $\chi_{\rm s}^{0}(\bm{q})=2\sum_{nm}\int\frac{d\bm{k}}{(2\pi)^{2}}F_{nm}(\bm{k},\bm{q})(1-D_{nm}(\bm{k},\bm{q}))$, where the quantum distance $D_{nm}(\bm{k},\bm{q})\equiv 1-|\braket{u_{n}(\bm{k}+\bm{q})}{u_{m}(\bm{k})}|^{2}$ is closely related to the quantum geometry. Quantum geometry suppresses $\chi_{\rm s}^{0}(\bm{q})$ at $\bm{q}\neq 0$ via nonzero quantum distance $D_{nn}(\bm{k},\bm{q})$, which is expanded as $\sim\sum_{\mu\nu}g_{n}^{\mu\nu}(\bm{k})q_{\mu}q_{\nu}+\cdots$ with the quantum metric. However, $\chi_{\rm s}^{0}(\bm{0})$ is not suppressed, and ferromagnetic fluctuation is relatively enhanced. The van Vleck susceptibility arising from $D_{nm}(\bm{k},\bm{q})$ for $n\neq m$ corresponds to the positional-shift contribution to the GES. For comparison, we also define the bare spin susceptibility without quantum geometry, $\chi_{\rm s:band}^{0}(\bm{q})=2\sum_{n}\int\frac{d\bm{k}}{(2\pi)^{2}}F_{nn}(\bm{k},\bm{q})$, in which magnetic fluctuation is determined by only the effective-mass term. By comparing these two quantities, we can elucidate the effects of quantum geometry.

In Figs. 4(a) and 4(b), we show $\chi_{\rm s}^{0}(\bm{q})$ and $\chi_{\rm s:band}^{0}(\bm{q})$ in the dispersive Lieb lattice model for $\mu_{\rm c}=0.7$. As expected by Fig. 3(a) showing the negative GES, $\chi_{\rm e}^{0:\mu\nu}=\lim_{\bm{q}=0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi^{0}(\bm{q})$, the bare spin susceptibility shows ferromagnetic fluctuation (Fig. 4(a)). However, antiferromagnetic fluctuation is obtained when we neglect the quantum geometry (Fig. 4(b)). Thus, we conclude that the quantum geometry induces the ferromagnetic fluctuation. It is emphasized that the maximum of DOS at $\mu_{\rm c}=0.7$ is not sufficient for the ferromagnetic fluctuation; the relative enhancement of $\chi_{\rm s}^{0}(0)$ compared to $\chi_{\rm s}^{0}(\bm{q}\neq 0)$ by the quantum distance/quantum geometry is essential. Note that the momentum dependence of spin susceptibility plays an essential role in unconventional superconductivity Moriya and Ueda (2000); Yanase et al. (2003). We also show $\chi_{\rm s}^{0}(\bm{q})$ for $(\mu_{\rm c},T)=(0.65,0.05)$ and $(\mu_{\rm c},T)=(0.65,0.01)$ in Figs. 4(c) and 4(d), respectively. Consistent with Fig. 3(d), we confirm the crossover from ferromagnetic to antiferromagnetic fluctuation as the temperature decreases.

Spin-triplet superconductivity.— Finally, we show that quantum-geometry-induced ferromagnetic fluctuation mediates spin-triplet superconductivity. To see this, we set the onsite interaction as $U=0.86$ and solve the linearized gap equation, $\lambda^{\rm t(s)}\Delta_{ll^{\prime}}(\bm{k})=-\frac{1}{N\beta}\sum_{\bm{k}^{\prime}\omega_{n}}\sum_{\{l_{i}\}}V_{ll_{1},l_{2}l^{\prime}}^{\rm t(s)}(\bm{k}-\bm{k}^{\prime})\mathcal{G}_{l_{1}l_{3}}(\bm{k}^{\prime},i\omega_{n})\Delta_{l_{3}l_{4}}(\bm{k}^{\prime})\mathcal{G}_{l_{2}l_{4}}(-\bm{k}^{\prime},-i\omega_{n})$, using the effective interaction obtained by RPA $V^{\rm t(s)}(\bm{q})$, which is $\simeq\frac{-1(3)}{4}U\chi_{\rm s}(\bm{q})U$ in single-band systems Moriya and Ueda (2000); Yanase et al. (2003); Anderson and Brinkman (1973); Miyake et al. (1986); Scalapino et al. (1986) but here extended to multi-band systems Sup . Here, $\mathcal{G}(\bm{k},i\omega_{n})$ is the Green function with the Matsubara frequency, $i\omega_{n}$. The instability of spin-triplet (singlet) superconductivity with the form factor $\Delta(\bm{k})$ is determined by the maximum eigenvalue $\lambda^{\rm t(s)}$. While the mean-field formalism overestimates the transition temperature, the dynamical effect of effective interaction is expected not to alter the superconducting symmetry, as in the cases of ³He Anderson and Brinkman (1973) and cuprates Moriya et al. (1990).

Figure 5(a) shows the spin susceptibility at $\mu_{\rm c}=0.7$ obtained by RPA. Ferromagnetic fluctuation is enhanced by the Coulomb interaction, as we see from the comparison to Fig. 4(a). Eigenvalues of the linearized gap equation are shown in Fig. 5(b) for spin-singlet extended-$s$-wave (orange line) and spin-triplet $p$-wave (blue line) superconductivity Sup . It is revealed that the spin-triplet superconductivity is stabilized around $\mu_{\rm c}\simeq 0.7$ and $0.9$ corresponding to the negative peak of GES in Fig. 3(a).

Especially, we obtain the largest eigenvalue at $\mu_{\rm c}=0.7$ where quantum geometry induces ferromagnetic fluctuation. Combined with the large DOS, the strong ferromagnetic fluctuation enhanced by interaction gives a large eigenvalue for spin-triplet superconductivity. Thus, we conclude spin-triplet superconductivity from quantum-geometry-induced ferromagnetic fluctuation.

Discussion.— In this Letter, we show that quantum geometry induces ferromagnetic fluctuation and results in spin-triplet superconductivity. The Fubini-Study quantum metric on the Fermi surface is an essential quantity for this mechanism of magnetism and superconductivity. Using the dispersive Lieb lattice model, we demonstrated that the non-Kramers band degeneracy on the Fermi surface plays the central role in enhancing the quantum-geometry-induced phenomena. In the diverse studies on unconventional superconductivity, the quantum geometry of electrons coupled to many-body effects has not been focused on. Stimulated by recent developments in the topology and geometry of quantum materials, we shed light on a route to spin-triplet superconductivity and, thereby, topological superconductivity.

A question of interest is whether our theory can be applied to other systems as well. To answer this, we have calculated the GES of Raghu’s model Raghu et al. (2008) for iron-based superconductors Sup . This model has the non-Kramers band degeneracy at the $\Gamma$ point. Also in this model, the quantum geometry induces ferromagnetic fluctuation due to the non-Kramers band degeneracy. In addition, we confirmed the quantum-geometry-induced ferromagnetic fluctuation in other models with the flat band and various band touching including the usual Lieb lattice model gfe . Thus, a wide range of materials with non-Kramers band degeneracy Zhang et al. (2019); Vergniory et al. (2019) are candidates for quantum-geometry-induced ferromagnetism and superconductivity. We expect that future material-specific studies will be stimulated by our work. The exploration of two-dimensional materials with high tunability, e.g., by band engineering through heterostructures, gate voltage, strain, and twist angle is also expected.

Supplemental Materials:
Spin-triplet superconductivity from quantum-geometry-induced ferromagnetic fluctuation

In this section, we show the detailed calculation of the spin, charge, and generalized electric susceptibility. While we focus on two-dimensional systems in the main text, the following discussion is written for systems in any dimension, including two and three-dimensional systems.

We consider the multi-band Hubbard model with SU(2) symmetry,

where $\hat{\bm{c}}_{\sigma}^{\dagger}(\bm{k})=(\begin{array}[]{ccccc}\hat{c}_{1\sigma}^{\dagger}(\bm{k})&\ldots&\hat{c}_{l\sigma}^{\dagger}(\bm{k})&\ldots&\hat{c}_{f\sigma}^{\dagger}(\bm{k})\end{array})$ is the creation operator of electrons with the wave vector $\bm{k}$, spin $\sigma=\uparrow\downarrow$, and the internal degrees of freedom $l$ such as orbitals and sublattices. The dimension of the internal degrees of freedom is represented by $f$. $\hat{n}_{l\sigma}(\bm{R})=\hat{c}^{\dagger}_{l\sigma}(\bm{R})\hat{c}_{l\sigma}(\bm{R})$ is the particle density operator for $l$ and spin $\sigma$ at position $\bm{R}$. The Fourier transform is defined by $\hat{c}^{\dagger}_{l\sigma}(\bm{R})=\frac{1}{\sqrt{N}}\sum_{\bm{k}}e^{-i\bm{k}\cdot(\bm{R}+\bm{r}_{l})}\hat{c}_{l\sigma}^{\dagger}(\bm{k})$. $H_{0}(\bm{k})$ is the matrix representation of the Fourier transform of hopping integrals with the internal coordinate $\bm{r}_{l}$. $U$ is the onsite Coulomb interaction, and $N$ is the volume of the system. The SU(2) symmetry is preserved in this model, since the spin-orbit coupling (SOC) and the magnetic field are absent. We ignore two-body interactions other than the onsite Coulomb interaction, such as an inter-orbital interaction, for simplicity.

The particle density operator for each spin is defined by,

where $\hat{\bm{c}}_{\sigma}(\bm{k},\tau)=e^{\tau\hat{\mathcal{H}}}\hat{\bm{c}}_{\sigma}(\bm{k})e^{-\tau\hat{\mathcal{H}}}$ is the imaginary-time representation with the imaginary time $\tau$. For later calculation, we also define the matrix elements of particle density operators for each spin as,

Using this, the spin susceptibility and the charge susceptibility are defined by,

with the bosonic Matsubara frequency $i\Omega_{n}$ and the inverse temperature $\beta$. $T_{\tau}$ represents the time-ordering product for $\tau$. Here, we used the relations ensured by the SU(2) symmetry $\chi_{\uparrow\uparrow}(\bm{q},i\Omega_{n})=\chi_{\downarrow\downarrow}(\bm{q},i\Omega_{n})$ and $\chi_{\uparrow\downarrow}(\bm{q},i\Omega_{n})=\chi_{\downarrow\uparrow}(\bm{q},i\Omega_{n})$. We also define the matrix representation $\bar{\chi}_{\sigma\sigma^{\prime}}(\bm{q},i\Omega_{n})$ of $\chi_{\sigma\sigma^{\prime}}(\bm{q},i\Omega_{n})$ using the matrix element written by,

The spin (charge) susceptibility of a noninteracting system, namely the bare spin (charge) susceptibility, can be written as,

We used the property of noninteracting systems, $\chi^{0}_{\uparrow\downarrow}(\bm{q},i\Omega_{n})=\chi^{0}_{\downarrow\uparrow}(\bm{q},i\Omega_{n})=0$, which is satisfied by the SU(2) symmetry of Hamiltonian. Here, we define the Green function of noninteracting systems $\mathcal{G}(\bm{k},i\omega_{n})=[i\omega_{n}-H_{0}(\bm{k})]^{-1}$ with fermionic Matsubara frequency $\omega_{n}$. Tr represents the trace for all degrees of freedom except for the spin. The matrix elements are obtained as,

After taking the sum of Matsubara frequency, we get,

with the dimension of the system $d$. We can calculate the band dispersion $\epsilon_{n}(\bm{k})$ and Bloch wave function $\ket{u_{n}(\bm{k})}$ by the eigenvalue equation, $H_{0}(\bm{k})\ket{u_{n}(\bm{k})}=\epsilon_{n}(\bm{k})\ket{u_{n}(\bm{k})}$.

We define the irreducible vertex,

In the random phase approximation (RPA), the spin (charge) susceptibility and its matrix representation are given by,

In this subsection, we derive the generalized electric susceptibility using the Kubo formula and local thermodynamics.

First, we show an alternative way to introduce charge susceptibility by using linear response theory. Based on the Kubo formula, the charge susceptibility of real-space and real-time representation $\chi_{\rm c}(\bm{r},t)$ with a position $\bm{r}$ and real time $t$ is defined by,

Here, $\braket{\hat{n}(\bm{r},t)}$ is the expectation value of the particle density operator, $\hat{n}(\bm{r},t)=\sum_{l,\sigma}\hat{c}_{l,\sigma}^{\dagger}(\bm{r},t)\hat{c}_{l,\sigma}(\bm{r},t)$, namely the charge density, $\phi(\bm{r},t)$ is an external scalar potential, and $\rho_{0}$ is the charge density in the absence of the external field. Also, $\hat{c}_{l,\sigma}^{\dagger}(\bm{r},t)=e^{i\hat{\mathcal{H}}t/\hbar}\hat{c}_{l,\sigma}^{\dagger}(\bm{r})e^{-i\hat{\mathcal{H}}t/\hbar}$ is the Heisenberg representation of creation operator with the Dirac constant $\hbar=1$ for the natural unit. After the Fourier transform with respect to real time, we get the frequency representation of the charge susceptibility,

Here, the infinitesimal $\delta$ ensures that the external field vanishes at $t\rightarrow-\infty$. We can also define the frequency representation of the charge density as,

Since the correlation function should decay away from the external field at $\bm{r}$, the integrand of Eq. (S20) contributes only when $\bm{r}-\bm{r}^{\prime}$ is sufficiently small. In contrast, we assume that the scalar potential spatially modulates on a length scale larger than that of the correlation function. Therefore, we can expand the scalar potential by $\bm{r}-\bm{r}^{\prime}$ as,

where and hereafter, we take the sum of repeated indices, such as $\mu=x,y,z$. For example, Eq. (LABEL:eq:n-rho) up to $n=2$ is explicitly written as,

Through the analytic continuation $\omega+i\delta\rightarrow i\Omega_{n}$, $\chi_{\rm c}(\bm{q},\omega)$ corresponds to Eq. (S7).

Here, we focus on $\chi_{\rm c}(\bm{q},\omega=0)$ for which the system does not depend on time. In other words, we focus on the particle density in equilibrium. When the system is metal and/or at a finite temperature, an equilibrium charge with an external electric field cannot be defined, since the electric current follows in metals or at finite temperatures. Therefore, we consider an insulator at zero temperature. Considering the Lehmann representation ¹¹1The Lehmann representation of charge susceptibility is written by $\displaystyle\chi_{\rm c}(\bm{q},\omega)=-\dfrac{1}{N\mathcal{Z}}\sum_{ab}\left(e^{-\beta E_{a}}-e^{-\beta E_{b}}\right)\dfrac{\bra{a}\hat{n}(\bm{q})\ket{b}\bra{b}\hat{n}(-\bm{q})\ket{a}}{\hbar\omega+E_{a}-E_{b}+i\delta}.$ (S24) Here, $\mathcal{Z}$ is the partition function. $\ket{a}$ and $E_{a}$ are the eigenstate and the eigenvalue of the many-body Hamiltonian $\hat{\mathcal{H}}$. , we see that the charge susceptibility satisfies the relationship $\chi_{\rm c}(\bm{q},\omega)=\chi_{\rm c}^{*}(-\bm{q},-\omega)$ which means $\chi_{\rm c}(\bm{q})=\chi_{\rm c}(-\bm{q})(=\chi_{\rm c}(\bm{q},0))$. Therefore, odd-order derivatives of $\chi_{\rm c}(\bm{q},0)$ with respect to $\bm{q}$ vanish. As a result, up to the second order of $\partial_{r_{\mu}}$, the equilibrium charge density is obtained as,

where $E_{\nu}(\bm{r})=\partial_{r_{\nu}}\phi(\bm{r})$ is the external electric field. In the insulator at zero temperature, the first term vanishes and this directly means that $\dfrac{1}{2}\lim_{\bm{q}\rightarrow 0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi_{\rm c}(\bm{q})$ is the electric susceptibility which gives the correction to the charge density, $\delta\braket{\hat{n}(\bm{r})}=\braket{\hat{n}(\bm{r})}-\rho_{0}$.

Next, to clarify the physical meaning of the quantity $\dfrac{1}{2}\lim_{\bm{q}\rightarrow 0}\partial_{q_{\mu}}\partial_{q_{\nu}}\chi_{\rm c}(\bm{q})$ in metals and/or at a finite temperature, we use the local thermodynamics Shitade et al. (2018, 2019); Gao et al. (2018); Daido et al. (2020); Kitamura et al. (2021). The following discussion is based on Ref. Daido et al., 2020. We consider the scalar potential arising from an inhomogeneous distribution of disorders, structural asymmetry, contact with a substrate, and so on, rather than an applied electric field. The electric field is assumed to be in a small region near $\bm{r}$ compared to the volume $N$. In this setup, while the local particle number depends on $\bm{r}$ due to the spatial variation of $\phi(\bm{r})$, the total particle number is assumed to be constant. Thus, the system is static and the charge current does not flow. The setup is discussed in more detail in Ref. Daido et al., 2020.