Introduction

The discovery of superconductivity in CeCu2Si21 marked the beginning of decades-long intense research into unconventional superconductivity2,3, encompassing cuprate4,5,6, iron-based7,8,9,10, and heavy fermion superconductors11,12,13. The proximity of superconductivity to antiferromagnetic (AF) quantum critical points (QCP) in these systems implicate AF fluctuations that proliferate at the AF QCP as the pairing glue, leading to unconventional superconductivity with a sign-changing superconducting order parameter2,3,14.

Experimental evidence for magnetically driven superconductivity in these systems include: (i) reduction of the magnetic exchange energy is much larger than the superconducting condensation energy15,16,17,18,19,20,21,22; (ii) the observation of a spin resonance mode, which in the spin-exciton scenario indicates a sign-changing superconducting order parameter2,23,24,25; and (iii) the persistence of two-dimensional (2D) high-energy AF fluctuations that resemble spin waves in magnetically ordered parent compounds19,26,27,28,29,30,31,32,33.

From an empirical perspective, it is important to identify whether these features are common in different unconventional superconductors, so that ingredients for a unified pairing mechanism may be established. Of the above observations, (i) is model-independent and has been verified for cuprate, iron-based and heavy fermion superconductors15,16,17,18,19,20,21,22; (ii) the spin resonance mode has been found in cuprate, iron-based, and heavy fermion superconductors, but their spin-excitonic nature needs to be separately tested and requires a quasiparticle origin of the magnetic excitations2,23,24; while (iii) has been established for cuprate and iron-based superconductors19,26,27,28,29,30,31,32,33, magnetic excitations in heavy fermion superconductors such as CeCu2Si2 are strongly three-dimensional (3D) at low energies18, and it is unclear whether they become 2D at higher energies.

CeCu2Si2 (S-type) is an archetypal heavy fermion unconventional superconductor, and is naturally located near a 3D AF QCP34,35. Upon the introduction of slight Cu deficiencies the system can be tuned to AF order (A-type), with an ordering vector τ ≈ (0.22, 0.22, 0.53) [Fig. 1a]36. Similar dispersive paramagnons up to E ≈ 1 meV were found to stem from τ in both A-type and S-type CeCu2Si218,37. While these dispersive AF fluctuations were discussed in terms of an effective Heisenberg model with short-range magnetic couplings18, magnetic order in A-type CeCu2Si2 was suggested to result from Fermi surface nesting36, and AF fluctuations in CeCu2Si2 were found to exhibit critical slowing down consistent with being near a spin-density-wave QCP35. In the superconducting state of CeCu2Si2, a spin gap forms with spectral weight built up just above it18, consistent with the formation of a spin resonance mode with an energy Er≈0.2 meV, which in the spin-excitonic scenario suggests magnetic pairing38,39. The recent discovery of fully-gapped superconductivity in CeCu2Si240,41,42,43,44,45,46,47,48 challenges the role of magnetic excitations in its superconducting state and calls for a more detailed understanding of its magnetism, including the origin and high-energy properties of its magnetic excitations.

Fig. 1: Crystal structure, phase diagram, and reciprocal space of CeCu2Si2.
figure 1

a Schematic phase diagram of CeCu2Si2. b Crystal structure of CeCu2Si2, with Ce3+ ions forming a body-centered tetragonal lattice. c [H, H, L] scattering plane of CeCu2Si2, with the shaded area representing a Brillouin zone for the Ce3+ ions. The black dots indicate magnetic Bragg peaks in A-type CeCu2Si236, and the dashed ovals are schematic intensity contours of low-energy AF excitations in our S-type CeCu2Si2. d Slice of the CeCu2Si2 Fermi surface with kx = ky, and the arrow represents the intraband nesting vector. The electronic structure used in this work is identical to that in ref. 46, and 3D plots of the Fermi surfaces are shown in Supplementary Fig. 7.

Full size image

In this work, by carrying out detailed inelastic neutron scattering measurements over large energy and momentum ranges, we uncover magnetic fluctuations up to E ≳ 5 meV in CeCu2Si2. While magnetic fluctuations below E ≈ 1.5 meV are strongly 3D and dispersive18,35, they become increasingly 2D with increasing energy and form an almost dispersionless column in energy. By comparing with theoretical calculations, we find magnetic excitations in CeCu2Si2 can be accounted for by intraband scattering of quasiparticles associated with the heavy electron band [Fig. 1d], and therefore allowing us to estimate the band renormalization by matching our calculations with experimental data. We expect this method to be broadly applicable in heavy fermion metals near magnetic criticality. The agreement between our experimental and theoretical results suggests that despite signatures of non-Fermi-liquid behavior34,35, the magnetic excitations in the normal state of CeCu2Si2 are reasonably captured by a LDA+U band structure with additional mass renormalization. Our discovery of almost 2D high-energy magnetic excitations in CeCu2Si2 is reminiscent of similar findings in cuprate and iron-based superconductors, and favors magnetic pairing with a sign-changing superconducting order parameter.

Results

Inelastic neutron scattering

Large single crystals of S-type CeCu2Si2 with Tc ≈ 0.5 K were grown using a vertical floating zone method49. Multiple crystals with a total mass of ≈12 g were co-aligned in the [H, H, L] scattering plane using the E3 four-circle neutron diffractometer at the Chalk River Laboratory (see Methods for details). Inelastic neutron scattering measurements were carried out using the Multi-Axis Crystal Spectrometer (MACS)50 at the NIST Center for Neutron Research, with fixed outgoing neutron energies Ef = 3 or 5 meV. We reference momentum transfer Q = (H, K, L) in reciprocal lattice units, with H = aQx/(2π), K = bQy/(2π), and L = cQz/(2π) (a = b ≈ 4.1 Å and c ≈ 9.9 Å). Ce3+ ions in CeCu2Si2 form a body-centered tetragonal lattice [Fig. 1b], and the corresponding Brillouin zone in the [H, H, L] plane is shown in Fig. 1c. Isotropic background intensities were estimated from regions with H = 0 and 1, and have been subtracted from our data.

Maps of the [H, H, L] plane at T = 1.6 K are compared in the left column of Fig. 2 for different energies, with the corresponding cuts along (H, H, 1.48) and (0.22, 0.22, L) shown in the middle and right columns. For E = 0.5 meV, magnetic excitations are relatively sharp, with spectral weight asymmetrically located around τ [Figs. 1c, 2a–c]. With increasing energy, the magnetic excitations gradually broaden, while maintaining the asymmetric distribution of spectral weight around τ [Fig. 2d–o]. Given the same asymmetric distribution is observed for multiple momentum and energy transfers, it is an intrinsic effect rather than a result of instrumental resolution. Two dispersive branches were previously observed at low energies in CeCu2Si218, with the branch closer to the zone boundary being increasingly dominant in intensity with increasing energy18. In our results, a single branch is resolved and can be identified as the dominant branch in previous work, while the weaker branch is unresolved and likely shows up as a shoulder in intensity. Clear modulation of magnetic intensity along (H, H, 1.48) persists up to the highest measured energy (E = 5.5 meV), with little or no magnetic intensity at the (0, 0, L) and (1, 1, L) positions; on the other hand, intensity along (0.22, 0.22, L) becomes weakly L-dependent for E ≥ 2 meV. These observations suggest that although CeCu2Si2 is close to a 3D AF QCP34,35, its high-energy magnetic excitations are almost 2D (see Supplementary Note 1 and Supplementary Fig. 1 for theoretical evidence of quasi-2D magnetic excitations), similar to cuprate and iron-based superconductors.

Fig. 2: Wave vector dependence of spin excitations in CeCu2Si2 as a function of energy.
figure 2

Constant-energy maps of the [H, H, L] plane for different energies are shown in the left panels, corresponding cuts along the [H, H, 1.48] direction obtained by binning data with 1.38 < L < 1.58 are shown in the middle panels, and corresponding cuts along the [0.22, 0.22, L] direction obtained by binning data with (0.17, 0.17) < (H, H) < (0.27, 0.27) are shown in the right panels. Data in this figure were measured using Ef = 5 meV. All vertical error bars represent statistical errors of 1 s.d.

Full size image

Dispersion of the magnetic excitations can be directly visualized in the energy-(H, H, 1.48) map in Fig. 3a. By fitting scans along (H, H, 1.48) using Gaussian peaks symmetrically positioned around (0.5,0.5,1.48) as shown in the middle column of Fig. 2, the magnetic dispersion along (H, H, 1.48) can be quantitatively extracted from E = 0.5 meV to 5.5 meV [Fig. 3a] (see Methods section for details). Consistent with previous observations18,35,37, the magnetic excitations are dispersive for E ≲ 1.5 meV, but at higher energies they form a column in energy away from the zone boundary. Such an evolution from dispersive to columnar magnetic excitations is unexpected for a local-moment magnetic system, but has been observed in itinerant magnetic systems, including heavily hole-doped iron pnictides51, Fe-doped MnSi352, and MnSi53.

Fig. 3: Dispersion of magnetic excitations in CeCu2Si2.
figure 3

a Energy-[H, H, 1.48] map at T = 1.6 K obtained by binning data with 1.38 < L < 1.58. b Theoretical energy-[H, H, 0.5] maps obtained from RPA calculations using our LDA+U band structure. c Comparison of experimental and theoretical dispersions, with the theoretical dispersion scaled down by the band renormalization factor r. d Dispersion of electron-like bands along the [kx, ky = kx, kz] direction. Data in this figure were measured using E = 5 meV. All horizontal error bars are least-square fit errors of 1 s.d.

Full size image

Theoretical calculations

To understand the origin of the magnetic excitations in CeCu2Si2 which extend up to at least E = 5.5 meV, we calculated magnetic excitations within the random-phase approximation (RPA) [Fig. 3b] using a LDA+U band structure (see Methods for details). As can be seen in Fig. 3b, despite extending over a larger energy scale, the calculated magnetic excitations are in good agreement with our experimental results, with the extracted dispersion consisting of a dispersive part at low energies and a columnar part at high energies (see Methods for details). By introducing an overall band renormalization factor r that scales our LDA+U band structure, excellent agreement between experimental and theoretical dispersions can be achieved [Fig. 3c]. We fit the dispersions with a linear-dispersing part that intersects a columnar part at Ecross, with the experimental dispersion further constrained to stem from τ = (0.22, 0.22, 0.53). Scaling the theoretical dispersion so that Ecross is identical for theoretical and experimental dispersions leads to r ≈ 40. For comparison, the normal state-specific heat coefficient from our LDA+U band structure is C/T ≈ 50 mJ/mol K2, which requires r ≈ 20 to match the experimental normal state value of ≈1.0 J/mol K2 for T → 0 K35.

The Ce f electrons in CeCu2Si2 participate in the formation of an electron and a hole Fermi sheet, both indispensable for its superconducting state46. Compared to specific heat measurements which contain contributions from all Fermi sheets, magnetic excitations measured by neutron scattering are sensitive to bands that exhibit good nesting properties. Therefore, the larger value of r inferred from our neutron scattering results suggests a larger renormalization factor for the well-nested heavy electron band, relative to the hole band. Such band-dependent renormalization effects have been discussed in the context of iron-based superconductors54,55,56, with bands exhibiting markedly different renormalization factors for different Fe 3d orbitals. In addition, while specific heat measurements are only sensitive to states within ~kBT (~0.1 meV for T = 1.6 K) of the Fermi level, magnetic fluctuations probed in our experiments involve states on the order of several meVs within the Fermi level. Therefore, stronger renormalization effects compared to our LDA+U calculations for states away from the Fermi level (relative to those within kBT of the Fermi level) can also contribute to the larger r values extracted from our inelastic neutron scattering measurements. Our findings illustrate the utility of neutron scattering measurements in extracting renormalization factors with band-specificity, complementing specific heat measurements. It should be noted that such a band-specificity is limited to the well-nested band, which dominates the magnetic excitation spectra, while other bands are effectively not probed. This method hinges on the fact that magnetic excitations arising from quasiparticles encode information on the band structure57,58,59, and can be especially useful in heavy fermion metals, for which angle-resolved photoemission spectroscopy measurements are challenging due to the small energy scales involved.

The structure of the heavy electron band along the (kx, ky = kx, 2π/c) direction [Fig. 3d] offers an intuitive understanding of the unusual magnetic dispersion in CeCu2Si2 [Fig. 3a, b]. The excellent agreement between our experiment and calculations demonstrates that the AF fluctuations in CeCu2Si2 are well-described as particle-hole excitations, in which quasiparticles below the Fermi level are excited to unoccupied states above the Fermi level. In our LDA+U calculations (without the renormalization by r), the band bottom is around −15 meV. Therefore, the crossover from dispersive to columnar behavior that occurs at around 60 meV is dominated by states above the Fermi level. Comparing states just above the Fermi level (E ≲ 30 meV) with those well above the Fermi level (E > 50 meV) reveals that the latter are much lighter, characteristic of conduction bands from none-f orbitals [Fig. 3d]. Therefore, the crossover from dispersive excitations at low-energies to columnar excitations at high-energies in CeCu2Si2 reflects the electron band’s reduction of f-orbital content above the Fermi level (see Supplementary Fig. 2 for the partial density states of CeCu2Si2 from our LDA+U calculations). For other values of kz, the band experiences a similar loss of f-orbital content for E > 50 meV, although more complex behaviors are seen closer to the Fermi level.

Temperature evolution of low-energy magnetic excitations

Maps of the [H, H, L] plane for E = 0.3 meV at T = 0.3 K (T < Tc) and 1.6 K (T > Tc) are compared in Fig. 4a, b, with their difference shown in Fig. 4c. Since E = 0.3 meV is above the energy window of the spin resonance mode in CeCu2Si2 (Er ≈ 0.2 meV)18, magnetic excitations in the superconducting and normal states are similar. Examining the difference of excitations measured at T = 0.3 K and 1.6 K nonetheless reveals a subtle shift of magnetic spectral weight toward the Brillouin zone center along the (H, H) direction upon cooling. The systematic presence of such behavior across multiple Brillouin zones [Fig. 4c] demonstrates this behavior to be an intrinsic property of CeCu2Si2. Cuts along the (H, H, 1.5) direction are compared in Fig. 4d for T = 0.3 and 1.6 K, and their difference is shown in Fig. 4e. Such a temperature-dependent shift is similar to the shift of ordering vector36 and magnetic excitations35 in A-type CeCu2Si2, which also move towards the Brillouin zone center upon cooling. These observations can be naturally understood now we have shown that magnetic excitations in CeCu2Si2 arise from heavy quasiparticles, and result from a combination of the intrinsic asymmetry of the magnetic dispersion and depletion of the electronic density of states near the Fermi level upon cooling (see Supplementary Note 2 and Supplementary Fig. 3 for details). Such a depletion occurs in A-type CeCu2Si2 due to a spin-density-wave gap, and in S-type CeCu2Si2 due to a superconducting gap.

Fig. 4: Temperature dependence of the E = 0.3 meV magnetic excitations across Tc.
figure 4

E = 0.3 meV constant-energy maps of the [H, H, L] plane for a T = 0.3 K, b T = 1.6 K, and c their difference. d Cuts along (H, H, 1.5) from maps in (a) and (b), obtained by binning data with 1.4 < L < 1.6. e The difference of cuts in (d). Data in this figure were measured using Ef = 3 meV. All vertical error bars represent statistical errors of 1 s.d.

Full size image

Discussion

Our experimental observation of magnetic excitations extending up to at least E = 5.5 meV in CeCu2Si2 demonstrates quasi-2D magnetic fluctuations with an energy scale much larger than the superconducting pairing energy to be a common feature in unconventional superconductors. As the bandwidth of magnetic excitations is captured by effective magnetic interactions, which in turn determine the saving of magnetic exchange energy in the superconducting state ΔEmag, the high-energy magnetic excitations observed in our work suggest a ΔEmag that is at least as large as previously reported18. The commonality of a much larger ΔEmag compared to the superconducting condensation energy ΔESC [Table 1] and the presence of a spin resonance mode2,23,24 across different families of unconventional superconductors, suggest a common pairing mechanism and favors a sign-changing superconducting order parameter such as d + d44 or s±45,46 for CeCu2Si2, rather than s++47. To conclusively distinguish between these scenarios, it is important to study the dispersion of the spin resonance mode in comparison with theoretical results under different pairing symmetries to test its spin-excitonic nature60,61,62. Our work presents a model that captures the normal state magnetic excitations of CeCu2Si2, and is consistent with magnetically driven superconductivity in CeCu2Si2.

Given that the magnetic excitations in A-type CeCu2Si2 and the superconducting and normal states of S-type CeCu2Si2 are similar for E ≳ 0.4 meV35,37, we expect the observed high energy excitations in our S-type CeCu2Si2 to also be present in A-type CeCu2Si2, as well as compositions in between. In addition, columnar spin excitations near τ are also seen in CeNi2Ge2, although compared to CeCu2Si2 no low-energy dispersive features were reported63. Since CeNi2Ge2 is paramagnetic and relatively far away from an AF QCP, this suggests that compared to the low-energy dispersive excitations, the columnar excitations at high energies are more robust upon tuning away from the QCP. In both cuprate and iron-based superconductors, the quasi-2D high energy magnetic excitations remain robust when tuning towards the superconducting state, while the low-energy excitations may change dramatically. Therefore, the prevalence of robust high-energy magnetic excitations suggests short-range 2D magnetic correlations provide a backdrop from which unconventional superconductivity emerges, while low-energy AF fluctuations and electronic structure can range from strongly 2D to having significant 3D features and are more tunable. This is analogous to conventional superconductors, in which a large Debye cutoff energy provides a backdrop for potential high-temperature superconductivity.

Table 1 The relationship between magnetic exchange and superconducting condensation energies in unconventional superconductors.16,17,18,19,20,21.
Full size table

Methods

Sample preparation and inelastic neutron scattering measurements

Several rod-shaped S-type CeCu2Si2 single crystal samples were grown using a vertical optical heating floating zone method [Supplementary Fig. 4a]49. To avoid Si excess which results in A-type CeCu2Si2, we used a high-pressure Ar atmosphere and a relatively small overheating of the floating zone beyond the melting temperature, which effectively reduces Cu evaporation. Using inductively coupled plasma atomic-emission spectroscopy, we determined the atomic percentage of our samples to be Ce:20.1(1)%, Cu:40.3(1)%, and Si:39.6(1)%. This stoichiometry is found to be consistent across several pieces of our samples, indicating they are dominantly S-type64, in agreement with previous transport measurements49. Specific heat measurements were carried out for several pieces of our CeCu2Si2 samples [Supplementary Fig. 4c], all exhibiting a specific heat jump below Tc ≈ 0.5 K, different from A-type and A/S-type CeCu2Si2 samples36,65. The magnitude of the specific heat jump exhibits some sample dependence, possibly due to parts of the samples being nonsuperconducting. The specific heat measurements evidence our CeCu2Si2 samples are dominantly S-type, without prominent signatures of antiferromagnetism. While the presence of a minority phase of A-type or A/S-type CeCu2Si2 is difficult to rule out, the high-energy magnetic excitations uncovered in our work should also be present in A-type and A/S-type CeCu2Si2, as discussed above. Therefore, the possible presence of such minority phases will not affect the conclusions of our work.

We cut the rod-like samples into segments of a few centimeters and used the E3 four-circle neutron diffractometer to identify single-grain pieces by mapping the ϕ and χ rotation angles, with the scattering angle 2θ adjusted to the scattering angle of an intense structural Bragg peak. We then co-aligned four such segments in the [H, H, L] plane, as shown in Supplementary Fig. 4b.

Inelastic neutron scattering measurements were carried out using the MACS50 at the NIST center for neutron research in Gaithersburg, MD. Our measurements were carried out using fixed Ef = 3 or 5 meV, with Be filters placed after the sample for both Ef and before the sample for Ef = 3 meV. MACS consists of 20 spectroscopic detectors, and by rotating the sample and shifting the detectors, a map of the scattering plane at a fixed energy transfer can be efficiently constructed. The double-bounce analyzers are vertically focused, while the monochromator is doubly focused. Instrumental energy resolutions at the elastic line are ΔE ≈ 0.14 meV for Ef = 3 meV, and ΔE ≈ 0.35 meV for Ef = 5 meV. Sample alignment is confirmed on MACS for the (110) and (002) structural Bragg peaks. As shown in Supplementary Fig. 5, our samples are reasonably well-aligned for sample arrays used in inelastic neutron scattering measurements.

Extraction of experimental and calculated dispersions

To extract the experimental dispersion of magnetic excitations in CeCu2Si2, cuts along (H, H, 1.48) were obtained by binning data with 1.38 < L < 1.58 and fit using two Gaussian peaks

$$I(x)={a}_{1}\exp \left(-\frac{{(x-\delta )}^{2}}{2{c}^{2}}\right)+{a}_{2}\exp \left(-\frac{{(x-1+\delta )}^{2}}{2{c}^{2}}\right).$$
(1)

The same expression with an additional constant term is used to extract the calculated dispersion. In these fittings, x = 0 corresponds to (0, 0), x = 1 corresponds to (1, 1) and δ is the fit peak center position. Representative fits to our experimental data are shown in the middle column of Fig. 3, and representative fits to our theoretical results are shown in Supplementary Fig. 6.

LDA+U band structure

The LDA+U band structure calculations were performed using the full-potential augmented plane-wave plus local orbital method as implemented in WIEN2k66. The Perdew–Burke–Ernzerhof exchange-correlation energy67 was used with spin-orbit coupling and an effective on-site Coulomb interaction U = 5 eV68. The orbital characters were obtained using WANNIER90 code69 via WIEN2WANNIER interface70. Our LDA+U band structure was used previously to study the pairing symmetry of CeCu2Si246, and is similar to band structures in previous LDA+U calculations40,45 and from the renormalized band approach71. As the f-electrons are itinerant in our LDA+U calculations, the obtained Fermi surfaces are “large”.

We note that there is a subtle difference in the band structures of refs. 45 and46, with the latter used in the calculations of this work. Comparing the heavy electron Fermi surfaces in these two works, there is an extra ring-like Fermi surface in ref. 46 around the Γ point (Supplementary Fig. 7). This difference results from details in implementing the calculations, and affects neither key features of the band structure nor the expected physics, which is dominated by the cylindrical heavy electron Fermi surface common to refs. 45 and46.

The partial density of states of CeCu2Si2 from our LDA+U calculations is shown in Supplementary Fig. 2. As can be seen, the Ce-f5/2 density of states is mainly located just above the Fermi level, and decreases rapidly above 50 meV, becoming increasingly small around 100 meV. Such an evolution of the partial density of states is consistent with the notion that a change of character of the band states causes the crossover from dispersive to columnar magnetic excitations. However, as the partial density of states contains contributions from all the electronic states, and not just the well-nested regions that give rise to the magnetic excitations, the signatures for such a change is not as clear in the partial density of states compared to Fig. 3d.

Calculation of magnetic excitations in CeCu2Si2

The bare magnetic susceptibility with four indices is:

$${[{\chi }_{0}]}_{\nu \nu ^{\prime} }^{\mu \mu ^{\prime} }({\bf{q}},i{\omega }_{n})=-\mathop{\sum}\limits_{{\bf{k}},{n}_{1},{n}_{2}}{a}_{{n}_{1}}^{\nu }({\bf{k}}){a}_{{n}_{1}}^{\nu ^{\prime} * }({\bf{k}}){a}_{{n}_{2}}^{\mu ^{\prime} }({\bf{k}}+{\bf{q}}){a}_{{n}_{2}}^{\mu * }({\bf{k}}+{\bf{q}})$$
(2)
$$\times \frac{1}{\beta }\mathop{\sum}\limits_{i{\omega }_{m}}\frac{1}{i{\omega }_{m}-{\epsilon }_{{n}_{1}}({\bf{k}})}\frac{1}{i{\omega }_{m}-i{\omega }_{n}-{\epsilon }_{{n}_{2}}({\bf{k}}+{\bf{q}})}$$
(3)

where \({a}_{n}^{\mu }({\bf{k}})\), ϵn(k) are the unitary matrices diagonalizing H0 and the energy dispersion, respectively. The sum over n is taken over the entire band index. Using the Matsubara frequency sum rule and the Fermi-Dirac function \({n}_{F}(\epsilon )=\frac{1}{{e}^{\beta \epsilon }+1}\), we get:

$${[{\chi }_{0}]}_{\nu \nu ^{\prime} }^{\mu \mu ^{\prime} }({\bf{q}},i{\omega }_{n})=-\mathop{\sum}\limits_{{\bf{k}},{n}_{1},{n}_{2}}{a}_{{n}_{1}}^{\nu }({\bf{k}}){a}_{{n}_{1}}^{\nu ^{\prime} * }({\bf{k}}){a}_{{n}_{2}}^{\mu ^{\prime} }({\bf{k}}+{\bf{q}}){a}_{{n}_{2}}^{\mu * }({\bf{k}}+{\bf{q}})\times \frac{{n}_{F}\left[{\epsilon }_{{n}_{2}}({\bf{k}}+{\bf{q}})\right]-{n}_{F}\left[{\epsilon }_{{n}_{1}}({\bf{k}})\right]}{i{\omega }_{n}+{\epsilon }_{{n}_{2}}({\bf{k}}+{\bf{q}})-{\epsilon }_{{n}_{1}}({\bf{k}})}$$
(4)

using iω → ω + iη, we have

$${[{\chi }_{0}]}_{\nu \nu ^{\prime} }^{\mu \mu ^{\prime} }({\bf{q}},\omega )=-\mathop{\sum}\limits_{{\bf{k}},{n}_{1},{n}_{2}}{a}_{{n}_{1}}^{\nu }({\bf{k}}){a}_{{n}_{1}}^{\nu ^{\prime} * }({\bf{k}}){a}_{{n}_{2}}^{\mu ^{\prime} }({\bf{k}}+{\bf{q}}){a}_{{n}_{2}}^{\mu * }({\bf{k}}+{\bf{q}})\times \frac{{n}_{F}\left[{\epsilon }_{{n}_{2}}({\bf{k}}+{\bf{q}})\right]-{n}_{F}\left[{\epsilon }_{{n}_{1}}({\bf{k}})\right]}{\omega +{\epsilon }_{{n}_{2}}({\bf{k}}+{\bf{q}})-{\epsilon }_{{n}_{1}}({\bf{k}})+i\eta }$$
(5)

The transverse RPA susceptibility for a multiband system is:

$${\hat{\chi }}_{{\rm{RPA}}}=\frac{{\hat{\chi }}_{0}}{1-\hat{U}{\hat{\chi }}_{0}}$$
(6)

This is in fact a Bethe-Salpeter equation where \(\hat{U}\) is a n2 × n2 matrix with orbital number n.

$${\hat{U}}_{rr}^{rr}=U,\,{\hat{U}}_{ss}^{rr}=U^{\prime} ,\,{\hat{U}}_{rs}^{rs}=J,\,{\hat{U}}_{sr}^{rs}=J^{\prime} \,(r\,\ne \,s).$$
(7)

On-site interactions \(U=U^{\prime}\)=0.25 eV, \(J=J^{\prime} =0\) eV, and a 20 × 20 × 20 k-mesh were used in our RPA calculations, similar to previous work45. The U values used in our RPA calculations are smaller than those in our LDA + U calculations to avoid the divergence of RPA magnetic susceptibility. The key features of the RPA susceptibility are mostly determined by the bare susceptibility χ0 (Supplementary Fig. 8), which already contains the essential features of the magnetic susceptibility, and are not strongly affected by the interaction term U.

Disclaimer: The identification of any commercial product or trade name does not imply endorsement or recommendation by the National Institute of Standards and Technology.