METHODS

Electrical transport measurements

High-quality ReO3 single crystals were synthesized via the chemical vapor transport method [39]. To comprehensively evaluate the intrinsic superconductivity and its sensitivity to pressure homogeneity, the high-pressure transport measurements were systematically conducted using different apparatuses and pressure-transmitting media (PTM). Measurements in the lower pressure range were performed by using palm-type CAC with liquid glycerol as the PTM (s1). The multi-anvil triaxial compression of the CAC provides a highly isotropic and nearly ideal hydrostatic environment, allowing the emergence of superconductivity to coincide precisely with the structural transition pressure. Higher pressure measurements were carried out by using BeCu-type DACs. To compare hydrostatic effects, two different PTM were employed in the DAC experiments: soft solid KBr (s2) and liquid Daphne 7373 (s3). While Daphne 7373 offers relatively good hydrostaticity under uniaxial compression, solid KBr provides relatively poor hydrostatic conditions, generating significant internal pressure gradients and local shear stresses. For sample s2, local regions with higher pressures show the superconducting transition at lower pressures. Pressure values inside the sample chamber were initially determined at room temperature by the ruby fluorescence method [40] or from the first-order diamond Raman peak shift [41]. The pressure variations upon cooling down in both DAC and CAC have been characterized by comparing the pressure calibration curves at room temperature and low temperatures.

X-ray diffraction measurements

The high-pressure SXRD experiments were performed at the BL15U1 station of the Shanghai Synchrotron Radiation Facility (SSRF) (λ = 0.6199 Å). Here, the ethyl-alcohol mixture with a ratio of 4:1 was used as a PTM. Using the geometric parameters calibrated by a CeO2 standard, the DIOPTAS software package [42] was used to integrate the powder diffraction patterns and then convert them to one-dimensional profiles. Rietveld refinements of XRD patterns were conducted using GSAS with EXPGUI packages [43]. The P-V relationship is fitted by using the Birch-Murnaghan equation:

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {{P}}\left( {{V}} \right) &=& \frac{{3{{{B}}}_0}}{2}\left[ {{{\left( {\frac{{{{{V}}}_0}}{{{V}}}} \right)}}^{7/3} - {{\left( {\frac{{{{{V}}}_0}}{{{V}}}} \right)}}^{5/3}} \right]\\ &&\times \left\{ {1 + \frac{3}{4}\left( {{{B}}_0^{\prime} - 4} \right)\left[ {{{\left( {\frac{{{{{V}}}_0}}{{{V}}}} \right)}}^{2/3} - 1} \right]} \right\},\\ \end{eqnarray*}\end{document}

where _V_0 is the cell volume at AP, V is the volume under pressure _P, B_0 is the bulk elastic modulus, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{B}}_0^{\prime} = 4$\end{document} is the derivative of the bulk elastic modulus at zero pressure. If we assume the rotation is rigid, the rotation angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{\theta }}$\end{document} can be estimated by the change of a/c via the following relationship:

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {{\theta }} = \arccos \left( {\sqrt 6 {{a}}/{{c}}} \right).\ \end{eqnarray*}\end{document}

WHH two-band model

The Werthamer-Helfand-Hohenberg (WHH) two-band model in the dirty limit [25]:

(3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} &&{{{a}}}_0[ {\ln {{t}} + {{U}}( {{{{D}}}_1{{h}}} )} ][ {\ln {{t}} + {{U}}( {{{{D}}}_2{{h}}} )} ]\\ &&+ {{{a}}}_2[ {\ln {{t}} + {{U}}( {{{{D}}}_2{{h}}} )} ] + {{{a}}}_1[ {\ln {{t}} + {{U}}( {{{{D}}}_1{{h}}} )} ] = 0,\\ \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{t}} = {{T}}/{{{T}}}{{ c}}$\end{document} is the reduced temperature, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{h}} = \hbar {{{\mu }}}0{{{H}}}{{{c}}2}/2 {\phi}0{{{k}}}{{B}}{{T}}$\end{document} is the dimensionless magnetic field with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\phi}0$\end{document} the magnetic flux quantum, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{D}}}{{n}}$\end{document} is the charge carriers diffusion coefficient. The coefficients a_0-a_2 are defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}0 = 2{{\omega }}/{{{\lambda }}}0$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}1 = 1 + {{{\lambda }}} - /{{{\lambda }}}0$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}2 = 1 - {{{\lambda }}} - /{{{\lambda }}}0$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{\lambda }}} - = {{{\lambda }}}{11} - {{{\lambda }}}{22}$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{\lambda }}}0 = \sqrt {{{\lambda }} - ^2 + 4{{{\lambda }}}{12}{{{\lambda }}}{21}} $\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{\omega }} = {{{\lambda }}}{11}{{{\lambda }}}{22} - {{{\lambda }}}{12}{{{\lambda }}}{21}$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{\lambda }}$\end{document} is the superconducting coupling constants matrix. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{U}}( {{x}} ) = {{\psi }}( {1/2 + {{x}}} ) - {{\psi }}( {1/2} )$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{\psi }}( {{x}} )$\end{document} is the di-gamma function. The dimensionless parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}_0$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}_1$\end{document}, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{a}}}_2$\end{document} describe the effective intraband and interband pairing interactions, while the scaled diffusivities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{D}}}_1$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{D}}}_2$\end{document} reflect the distinct carrier mobilities and scattering rates for the two active bands. The pronounced disparity in the extracted diffusivities (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{D}}}_1 \gg {{{D}}}_2$\end{document}) signifies a strong interband asymmetry, which effectively suppresses the pair-breaking effect in the stronger band and leads to the robust upper critical field observed in the R-I phase. All extracted fitting parameters at selected pressures are summarized in Table S4.

Resistance fitting

To probe the lattice dynamics, we analyzed the normal-state resistance by the following equation [18,26]:

(4) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {{R}}\left( {\boldsymbol{T}} \right) &=& {{{R}}}_0 + {{{R}}}_{{\rm D}}{\left( {\frac{{\boldsymbol{T}}}{{{{{\Theta }}}_{{\rm D}}}}} \right)}^5\\ &&\times \mathop \int \nolimits_0^{{{{\Theta }}}_{{\rm D}}/{{T}}} \frac{{{{{z}}}^5{{dz}}}}{{\left( {{{{e}}}^{{z}} - 1} \right)\left( {1 - {{{e}}}^{ - {{z}}}} \right)}}\\ &&+\, {{{R}}}_{{\rm E}}{\left[ {\frac{{{T}}}{{{{{\Theta }}}_{{\rm E}}}}{{\sinh }}^2\left( {\frac{{{{{\Theta }}}_{{\rm E}}}}{{2{{T}}}}} \right)} \right]}^{ - 1}, \end{eqnarray*}\end{document}

where the first term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{R}}}0$\end{document} is the residual resistance, and the second and third terms are the Bloch-Grüneisen and Einstein expressions, which represent the electron scattering by acoustic and optical phonons, respectively. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{R}}}{{\rm D}}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{R}}}{{\rm E}}$\end{document} are resistance constant and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{\Theta }}}{{\rm D}}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${{{\Theta }}}_{{\rm E}}$\end{document} are Debye and Einstein temperatures, respectively. The fitting results of ambient sample and CAC s1 are shown in Fig. S2a and c.

DFT calculations

The system was theoretically calculated by employing the first principles method as implemented in QUANTUM ESPRESSO (QE) package [44], and the equilibrium structure, phonon dispersions, EPC, and superconductivity of the compound were obtained as a function of the external hydrostatic pressure from 0 to 45 GPa. The plane-wave kinetic-energy cutoff and the energy cutoff for charge density were set as 100 and 400 Ry, respectively. The Brillouin zone (BZ) was sampled by a k-point mesh of 8 _× _ 8 _× _ 8, and a Methfessel-Paxton smearing width of 0.02 Ry was adopted for calculating the self-consistent electron density. The dynamic and EPC matrix elements were computed with an 8 _× _ 8 _× _ 8 q mesh. The phonon and EPC were evaluated quantitatively by employing the density functional perturbation theory (DFPT) and Eliashberg theory [45,46].

The electronic structure, including the DOS, was calculated using the Vienna Ab initio Simulation Package (VASP) [47]. The exchange-correlation interaction was described using the Perdew-Burke-Ernzerhof functional within the generalized gradient approximation [48], and the electron-ion interaction was treated using the projector-augmented-wave method [49]. The plane-wave energy cutoff was set to 520 eV. The Brillouin zone was sampled using automatically generated k-point meshes with a target grid density of 400 k-points per Å−3 of reciprocal-cell volume. COHP and integrated COHP (ICOHP) analyses were performed using the LOBSTER code [50] to evaluate the pressure evolution of the covalent bonding strength of the Re-O bonds. The wave functions used for the LOBSTER calculations were obtained from the corresponding VASP calculations. To ensure the convergence of the wave functions, the electronic self-consistent-field cycles were converged until the total-energy difference between successive steps was less than 10−7 eV, and no symmetry operations were imposed in these calculations.

Magnetic susceptibility measurement

Magnetic susceptibility measurements were performed in an Magnetic Property Measurement System from Quantum Design. For the initial measurements (sample s4), a miniature non-magnetic BeCu-type DAC with a 400 μm culet and a BeCu gasket pre-compressed to 5 GPa were used, with the initial sample chamber thickness being ∼ 60 μm. To quantitatively evaluate the bulk nature of pressure-induced superconductivity in the R-I phase of ReO3, we estimated the superconducting volume fraction (f) by using the ZFC magnetic susceptibility data. For sample s4 at 23 GPa, the ZFC diamagnetic signal was measured to be m = −1.5 × 10−6 emu under an applied magnetic field of H = 10 Oe. The diameter of the sample was D ≈ 175 μm. Under 23 GPa, the sample thickness was corrected to 46 μm based on the equation of state. This yields an estimated geometric volume of V ≈ 1.11 × 10−6 cm3. Because the sample is approximately a flat disk (d/D ≈ 0.26) positioned perpendicular to the applied magnetic field, the demagnetization factor N is ∼ 0.69. The intrinsic superconducting volume fraction f in cgs units is derived as:

(5) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {{f}} = \frac{{\frac{{4{{\pi m}}}}{{{{HV}}}}}}{{ - 1 + {{N}}\times \frac{{4{{\pi m}}}}{{{{HV}}}}}}. \end{eqnarray*}\end{document}

The superconducting volume fraction is deduced to be ∼ 78%.

To further strictly verify the bulk superconductivity at the maximum _T_c dome while minimizing non-hydrostatic effects, we performed a second set of measurements (sample s5) using a DAC equipped with a 300 μm culet and a rhenium (Re) gasket. For this sample, solid KBr was used as the pressure transmit medium. Due to the smaller culet size and the inclusion of the pressure medium, the sample dimensions were strictly constrained, and a larger excitation magnetic field of H = 50 Oe was applied to ensure a sufficient signal-to-noise ratio.

For sample s5 at 25 GPa, the measured ZFC diamagnetic signal was m = −1.2 × 10−6 emu. The sample diameter was D ≈ 75 μm, and the thickness was corrected to d ≈ 27 μm, yielding a cylindrical volume of V ≈ 1.19 × 10−7 cm3. With an aspect ratio of d/D ≈ 0.36, the demagnetization factor N is calculated to be ∼ 0.62. Applying the identical demagnetization correction formula, the intrinsic superconducting volume fraction f for s5 at 25 GPa reaches an astonishing ∼ 98% (107% for 30 GPa). This robust macroscopic diamagnetic response provides direct evidence of bulk superconductivity in the R-I phase.