Uncovering origins of heterogeneous superconductivity in La3Ni2O7

NV ground state Hamiltonian

Each NV centre features an S = 1 electronic spin ground state, described by the Hamiltonian \(H_=D_\mathrmgsS_Z^2\), where \(\widehatZ\) corresponds to the quantization axis along the NV axis^51,52. In the absence of external perturbations, the zero-field splitting, D_gs = (2π) × 2.87 GHz, captures the energy difference between the m_s = 0⟩ spin sub-level and the two degenerate m_s = ±1⟩ sub-levels (Fig. 1d); initialization to the m_s = 0⟩ state is achieved by optical pumping using a 532 nm laser^40,53,54,55.

Both magnetic field and stress alter the relative energies of the spin sub-levels (Fig. 1d). In particular, magnetic fields B couple through the Zeeman Hamiltonian H_B = B · S (refs. ^56,57,58,59), where γ is the NV spin gyromagnetic ratio, and stress \(\overleftrightarrow\sigma \), which is a rank-2 tensor with six independent elements, couples through the stress Hamiltonian, \(H_\sigma =\varPi _ZS_Z^2+\varPi _X(S_Y^2-S_X^2)+\varPi _Y(S_XS_Y+S_YS_X)\) (refs. ^18,60,61). Here, (Π_X, Π_Y, Π_Z) are functions of the stress tensor that transform according to the irreducible representations of the NV point group: \(\varPi _X=-(b+c)(\sigma _YY-\sigma _XX)+(\sqrt2b-\fracc\sqrt2)(2\sigma _XZ),\,\varPi _Y=-(b+c)(2\sigma _XY)\,+\)\((\sqrt2b-\fracc\sqrt2)(2\sigma _YZ)\) and Π_Z = (a₁ − a₂)(σ_XX + σ_YY) + (a₁ + 2a₂)σ_ZZ, where a₁, a₂, b, c correspond to the stress susceptibilities of the NV centre^18,62.

Measurement protocols

As superconductors support metastable currents, they often show hysteretic behaviour. Thus, it is crucial to keep track of the thermodynamic history of the sample. We use the following measurement protocols, which always begin at T > T_c. Each of these steps is also shown graphically in Extended Data Fig. 1b.

Diamagnetism

For the B/H measurements performed in both the main text and in the Extended Data, we use the following protocol:

1. 1.

   ZFC: we begin by cooling the sample below the transition temperature T_c without any applied field (H = 0). Electrical transport is recorded during this step.

2. 2.

   The field H is turned on. We typically choose fields H of the order 100–150 G, applied along \(\widehatz\).

3. 3.

   Field warming: the temperature is then raised in steps, and the transport and magnetometry measurements are taken simultaneously during this warming procedure.

Flux trapping

We use two different protocols to measure flux trapping. The first protocol, detailed below, is used in the flux trapping study of La₃Ni₂O₇ (Fig. 2e and Extended Data Fig. 7b):

1. 1.

   ZFC: we perform ZFC to 20 K, and then perform NV magnetometry to collect a reference ZFC spectrum.

2. 2.

   Field cooling: a 150 G field is turned on at 150 K, where the sample is in the normal state. Then, we cool below T_c to 20 K.

3. 3.

   The field is turned off, and NV magnetometry is performed to collect the field-cooled spectrum.

This method is further detailed below in the ‘Flux trapping and analysis’ section.

In the next section, we perform benchmarking experiments on Bi₂Sr₂CaCu₂O_8+δ. For this material, we use a slightly different flux trapping protocol, in which we force the flux into the system while it is already in the superconducting state:

1. 1.

   ZFC: the sample is cooled below T_c without any applied field.

2. 2.

   The field is turned on. The field should be at a magnitude great enough that it begins to penetrate through the superconductor.

3. 3.

   The field is turned off, and NV magnetometry is performed. Again, without the presence of the external field, any remnant non-zero magnetic field must be present as a result of trapped flux.

Field cooling

The last measurement protocol that we use is the field cooling protocol (Extended Data Fig. 6d):

1. 1.

   Similar to the first flux trapping protocol, the field is turned on while the sample is in the normal state. Here, the sample is cooled in a sequence of temperature increments, during which NV magnetometry is performed. Here, due to the flux trapped as the system transitions into the mixed phase, we expect a different B/H response than in the ZFC-FW protocol.

Notation

We reserve uppercase indices I, J and coordinates X, Y, Z for the lab frame defined by the culet and imaging system. More precisely, \(\widehatZ\) is perpendicular to the culet and \(\widehatX,\widehatY\) are defined by the orientation of the CCD imaging plane, see Extended Data Fig. 1a. This is the most natural frame for mechanical considerations.

The response of the NV centre is most naturally described in the frame adapted to its symmetry axis, as is conventional¹⁸. As there are four NV orientations, there are four such frames, which we label as g = 0, 1, 2, 3, when it is necessary to disambiguate (Extended Data Fig. 8b). We adopt lowercase indices i, j and coordinates x, y, z for NV frames. Thus, the NV^(g) group has symmetry axes \(\widehatz^(g)\). We define the global NV frame to be that of the g = 0 group and typically suppress the g label for quantities in that frame.

Each NV group contributes two resonance peaks to the ODMR spectrum observed in a pixel. We define D_g to be the frequency of the midpoint of the two peaks associated with group g and E_g to be their half-difference. We note E = Δν, as defined in Fig. 1d, and we use the two representations interchangeably.

Benchmarking the NV-DAC platform

To validate the NV-DAC system for characterizing high-pressure superconductivity, we benchmark the method using the well-characterized cuprate superconductor Bi₂Sr₂CaCu₂O_8+δ (BSCCO)⁶³. Our strategy combines conventional transport measurements with spatially resolved magnetic imaging using NV centres embedded directly in one of the culets of the DAC. Similar to the study of La₃Ni₂O₇ in the main text, this integrated approach allows us to precisely determine the superconducting transition and directly image the associated Meissner effect, demonstrating the effectiveness of NV centres for probing superconductivity under extreme conditions.

The sample is mounted inside a DAC equipped with four electrical leads for transport measurements and an off-sample microwave antenna for NV-based magnetic field sensing, with rubies to calibrate pressure (Extended Data Fig. 2a). Transport measurements show a superconducting transition at approximately 92 K at 2.1 GPa (Extended Data Fig. 2d).

Zero-field-cooled diamagnetism

Following ZFC to a base temperature of 7.7 K, we apply a magnetic field of 158 G. We then record a wide-field map of B/H by imaging the local magnetic responses from the shallow NV centres in the diamond culet beneath the BSCCO sample. Extended Data Fig. 2b shows representative raw spectra from a one-dimensional cross-section (indicated in Extended Data Fig. 2a), analogous to Fig. 1. Far away from the sample (at the top of the waterfall plot), the [111] NV resonances are split exactly by the expected value for H = 158 G. Approaching the sample edge, this splitting sharply increases, consistent with the compression of magnetic field lines near the boundary of a superconductor. Moving further into the sample interior, the resonances converge, indicating diamagnetic suppression. By fitting these resonances pixel by pixel, we obtain a 2D B/H map illustrating BSCCO in the Meissner state (Extended Data Fig. 2c). The NV centres at the sample centre exhibit B/H ≈ 0, confirming the expected bulk superconductivity. The enhancement due to field deformation at the sample edge is correspondingly strong.

Field warming

With the field H = 158 G fixed, we then warm the sample stepwise up to 95 K. Selected temperature-dependent B/H maps are shown in Extended Data Fig. 2c. The diamagnetic response gradually fades as temperature increases, with the Meissner effect disappearing around 90 K, consistent with the transport-derived T_c of 92 K.

Flux trapping

Following another ZFC to 7.7 K, we apply a magnetic field H = 200 G (above the lower critical field of BSCCO at low temperature, H_c1 ≈ 100 G; ref. ⁶⁴) and then quench the field. Imaging the remnant magnetic field B (Extended Data Fig. 2e) shows that up to 60 G is trapped within the sample region, especially near the boundaries of the sample, with no field trapped outside the sample. This result confirms strong magnetic flux trapping in our pressurized BSCCO sample.

Higher pressure

Increasing the pressure to 7.7 GPa and performing a resistance measurement shows a slight reduction in T_c to approximately 89 K (Extended Data Fig. 2d). At 7.7 K and an applied field of 158 G, we again observe diamagnetic suppression within the sample and magnetic field enhancement at its edges after ZFC. However, the overall response is weaker than at 2.1 GPa. Furthermore, the diamagnetic response exhibits a striped texture, possibly because of pressure-induced micro-fractures in the thin sample flake, which were not observable under optical microscopy, highlighting the unique diagnostic abilities of NV magnetometry.

These benchmarking experiments with BSCCO demonstrate that NV-DAC magnetometry reliably detects superconducting transitions and directly images superconducting diamagnetism under high-pressure conditions. Thus, this integrated method provides a robust and sensitive platform for exploring superconductivity in new materials under extreme conditions.

Electronic transport for S1, S2 and S3

Electrical resistance was measured using the four-probe method with platinum leads pressed onto the sample (Extended Data Fig. 3c). In the main text, we show select transport measurements from S1, S2 and S3. Here, we present an expanded set of resistance curves. In Extended Data Fig. 3a, we list resistance curves taken over the full range of measured pressures of S1. The first signatures of a resistance drop appear at \(\overline\sigma _ZZ=16\,\mathrmGPa\). Extended Data Fig. 3b shows resistance measurements taken at \(\overline\sigma _ZZ=23\,\mathrmGPa\) with different permutations of transport leads. The transport channel that probes the resistance of the sample between leads 2 and 3 shows essentially no superconducting transition. The region spanned by these leads is precisely the region of S1, which shows the least diamagnetism at this stress point (Figs. 1c and 2a). By contrast, the regions between leads 3 and 4, and likewise between 4 and 5, show more meaningful diamagnetism. Correspondingly, there is an appreciable transport signature in resistance measurements taken across these regions. This emphasizes the ability of our method to locally diagnose superconductivity in a multimodal fashion.

We show the transport for all pressures for S2 (Extended Data Fig. 3c) and S3 (Extended Data Fig. 3d). Although S2 shows superconducting behaviour starting at \(\overline\sigma _ZZ=13\,\mathrmGPa\), the transport of S3 shows no visible superconducting transition over the entire pressure range.

Low-pressure ZFC-FW B/H maps for S1, S2 and S3

ZFC-FW maps of diamagnetism in the main text were shown for pressures above the superconducting threshold for all samples. Here, we show low-pressure maps of B/H for S1, S2 and S3 (Extended Data Fig. 4a–c). The maps are measured at temperatures both below and above the high-pressure T_c of about 80 K. All high-temperature and low-pressure, B/H maps show B/H ≈ 1 over the sample region, until the normal stress exceeds 15 GPa and temperature drops below T_c (for example, S1 B/H map at 16 GPa and 35 K and S3 B/H map at 20 GPa and 7.5 K), which is in agreement with the phase diagram shown in the main text. We also show the EDX mapping of S2 in Extended Data Fig. 4d.

Extended temperature sweep ZFC-FW B/H maps for S1

In the main text, we selected five temperature points to show the evolution of the magnetic response of S1 with temperature. Here, we focus on S1 at \(\overline\sigma _ZZ=21\,\mathrmGPa\) and present a series of 18 maps of B/H from 20 K to 101 K (Extended Data Fig. 5). All regions of non-trivial magnetic response B/H ≠ 1 vanish above the temperature at which a kink appears in the electronic transport (Fig. 2c) at T_c ≈ 80 K.

Field cooling magnetic response for S1

The magnetic response maps shown in the main text were obtained in the ZFC-FW modality. We have also mapped B/H while field cooling from T > T_c at fixed H = 150 G at \(\overline\sigma _ZZ=23\,\mathrmGPa\) in sample S1 (Extended Data Fig. 6a,c). The figure also shows a complete set of ZFC-FW data for comparison in the same stress environment (Extended Data Fig. 6b). Broadly, the most salient feature of this comparison is that the maps at the same temperature in the two protocols look similar, but the magnetic response is overall weaker in the field cooling measurement. That is, if a pixel shows a dia- or paramagnetic response that turns on below T_c during the ZFC-FW protocol, it also does in the field cooling protocol with the same sign but typically smaller magnitude. This is consistent with the behaviour in dirty superconductors in which flux penetration in the mixed state becomes hysteretic because of pinning. We note that the electrical resistance curves measured during the two protocols are indistinguishable, suggesting that the extra flux penetrating the sample in the field cooling configuration is not significantly modifying the network carrying supercurrent.

Flux trapping and analysis

The field-cooled flux-trapping protocol, which was previously described, was used for this experiment. Here, we note that the trapped flux that was frozen in across the superconducting transition during the field cooling procedure led to only a small signal. Challenges in directly measuring the flux included: (1) inhomogeneities in E₀ due to stress and electric-field disorder, which made the magnetic field contributions difficult to identify, and (2) in the small-field regime, the NV groups experience minimal E₀, leading to reduced distinguishability of all eight groups. As a result, the relevant peaks, which correspond to the NV⁽⁰⁾ group, may be hard to fit. To rectify this issue, a differential approach was adopted, in which we compare ODMR spectra taken after the following steps (Extended Data Fig. 7a):

1. 1.

   Perform ZFC to T = 20 K and H = 0 G and collect the ZFC spectrum.

2. 2.

   Turn the field on and set H = 150 G.

3. 3.

   Cycle back above T_c to T = 150 K and perform field cooling back to T = 20 K.

4. 4.

   Turn off the field and collect the field-cooled spectrum.

The first ZFC spectrum provides a reference measure of the local stress and electric field contributions to the peak splitting, and ultimately, we compute the trapped field by subtracting this in quadrature¹⁸ from the final field-cooled spectrum:

This approach succeeds for 85% of the pixels in sample S1. There are two failure modes: first, some of the low-field spectra simply do not have identifiable peaks (13.3% of the pixels; Extended Data Fig. 7b, purple). Second, a small number of pixels (1.5%) have peaks that apparently shift inwards in the field-cooled final spectra by more than our noise floor of 5.6 MHz, because of electric field disorder in the absence of fields and stress; these pixels are marked in red. We display raw ODMR data comparing the ZFC and field-cooled spectra in Extended Data Fig. 7c.

In the main text, we discussed how regions of strong flux trapping in S1 correlated with regions of strong diamagnetism as observed in ZFC-FW. At first, this seems counterintuitive: we expect flux to be preferentially trapped at crystalline defects, suggesting that regions with weak diamagnetism (that is, less pure or homogeneous superconductivity) should correlate with flux trapping. However, because La₃Ni₂O₇ superconductivity is already quite heterogeneous, we expect all parts of the crystal to have sufficient defects to trap flux. In this filamentary limit, we are more concerned with superconductivity being present to trap flux. So, we expect—and observe—that regions of greater diamagnetism trap more flux.

Sum of flux

Magnetostatics dictates that the total magnetic flux through the entire sample chamber should remain the same, regardless of whether a sample is present or not. That is, ∑_iB_i ≈ NH, where the index i runs over all N pixels in a spatial map, B_i is the field measured at a given pixel and H is the applied field. Equivalently, \(\sum _i\fracB_iNH\approx 1\). To this end, we explicitly calculate the value of \(\sum _i\fracB_iNH\) (across the entire sample chamber) for all three samples and show the values in Extended Data Table 1. For all temperatures and pressures, across all three samples, the value of \(\sum _i\fracB_iNH\) is extremely close to unity. To ensure that this is not just a quantitative artefact, we also compute \(\sum _i\fracB_iNH\) over just the superconducting regions and observe values consistent with an approximately 2−4% diamagnetic suppression in the superconducting region.

High-pressure, optical and EDX experimental details

Sample preparation

Single-crystal sample growth is detailed in refs. ^9,65. To prepare the sample for the high-pressure NV-DAC measurement, a piece of a La₃Ni₂O₇ single crystal was shattered by compression between two glass slides. Shards were subsequently chosen for prescreening by EDX, based on their size and shape relative to the DAC sample chamber. We then selected samples with varying degrees of stoichiometric inhomogeneity for the studies of S1, S2 and S3.

Creation of NV centres

For the high-pressure NV measurement, we use 16-sided standard design diamond anvils with [111]-crystal cut culets (Almax-easyLab). These anvils are polished from synthetic type Ib ([N] ≤ 200 ppm) single-crystal diamonds. The culet diameters are in the range of 300–400 μm. Similar to previous work, we perform ¹²C⁺ ion implantation with multiple energy shots up to 450 keV with dosage up to 9.9 × 10¹¹ cm⁻² (CuttingEdge Ions) to generate a layer of vacancies of about 500 nm from the culet surface²¹. Following implantation, we perform vacuum annealing of the diamond anvils (at P < 10⁻⁶ mbar) in a custom-built furnace at a temperature >850 °C for 12 h. During the annealing step, the vacancies become mobile and probabilistically form NV centres. The NV diamond is glued to the cylinder side of a miniature DAC machined out of beryllium copper (Cu–Be). A Pt foil with a width of 20 μm is incorporated near the culet to apply microwave radiation for spin-state control. We further use Pt foil to serve as the transport connection pads by compression against the sample crystal surface under pressure. On the opposing (piston) side, we glue a type Ia diamond anvil with [100]-crystal cut. A double-sided insulating gasket, fashioned out of rhenium (Re) foil and cubic boron nitride (cBN) mixed with epoxy, is used to isolate the transport leads and the microwave wire.

Wide-field imaging in the AttoDRY800

The NV magnetic and stress imaging of sample S1 and S2 are performed in a dry, closed-cycle cryostat (Attocube AttoDRY800) integrated with a custom-built wide-field fluorescence microscope. The continuous excitation beam from a 532-nm laser (Sprout Solo-10W) is diffused and refocused at the back of an objective lens (20× Mitutoyo Plan NA 0.42). The field of view is in the range of 200–300 μm. The fluorescence collected by the objective is separated from the excitation path by a dichroic mirror and collected by an electron multiplying charge-coupled device (CCD) (Princeton Instrument PROEM-HS: 512B eXcelon). A data acquisition device (National Instruments USB-6343) is used for fluorescence counting as well as modulation of a microwave source (SRS SG384) for continuous wave ODMR measurements. The microwaves are amplified (Mini-Circuits ZHL-25W-63+) and subsequently channelled to the sample using coaxial cables and a custom-built circuit board. A vector external magnetic field is applied by running currents from a d.c. power supply (Keithley 2200-30-5 and Keysight e36154a) into custom-built electromagnet coils that sit outside the cryostat shroud. The \(z\prime \)-coil is circular and closer to the shroud. The maximum applied field is 240 G. The \(x\prime \) and \(y\prime \) coils are rectangular Helmholtz coils, which are further away from the shroud. Their maximum field is 50 G. A lock-in amplifier (Stanford Research Systems SR860) in combination with a voltage-controlled current source (Stanford Research Systems CS580) is used to measure the four-point resistance of the sample.

Wide-field imaging in the AttoDRY2100

NV magnetic and stress imaging of sample S3 was performed in a helium-gas-cooled closed-cycle cryostat (Attocube AttoDRY2100) integrated with a custom-built wide-field fluorescence microscope. A continuous 532 nm excitation beam (Gem SMD12-INV-v7) was spatially homogenized and focused through the back aperture of a long-working-distance objective (LT-APO/ULWD/VISIR, NA 0.35), yielding a field of view of 100–300 μm. Fluorescence emission was separated from the excitation path by a dichroic beam splitter and collected by an electron-multiplying CCD camera (Princeton Instruments ProEM-HS:512BX eXcelon). Signal acquisition and microwave modulation for continuous wave ODMR were performed using a National Instruments USB-6343 data acquisition card. Microwaves from an SRS SG384 source were amplified (Mini-Circuits ZHL-16W-43+) and delivered to the sample via coaxial cables and a custom-designed printed circuit board. A superconducting vector magnet integrated within the cryostat provided independently controlled magnetic fields up to 1 T along the \(x\prime \)-, \(y\prime \)– and \(z\prime \)-directions. Four-point resistance measurements were carried out with a Zurich Instruments MFLI lock-in amplifier, in combination with a Stanford Research Systems CS580 voltage-controlled current source.

Stoichiometry measurements by EDX spectroscopy

The EDX spectroscopy measurements were performed using a Gemini 360 field emission scanning electron microscope and an Oxford Instruments Ultim Max detector. The electron energy is fixed at 15 keV, and the beam current at 60 mA. The measured elements are calibrated by well-characterized standards La₂O₃ and Ni.

Calibrating σ
_ZZ and ruby pressure

In the main text, we stress the importance of using \(\overline\sigma _ZZ\) as a marker of sample stress, as opposed to ruby pressure. Here, we show a reasonably linear tracking of ruby pressure with \(\overline\sigma _ZZ\) across all three samples (Extended Data Fig. 1c). This emphasizes our point that although ruby pressure is a reasonable proxy for the average stress conditions of the sample, it is imprecise in two key ways. First, ruby pressure is lower by a factor of about 0.73 compared with \(\overline\sigma _ZZ\). This reflects the non-hydrostaticity of the salt pressure medium; the in-plane stress is lower than the stress along the loading axis as we might expect for a solid. Second, ruby pressure fails to capture the several-GPa local variations of normal stress present across the samples.

Calibrating magnetic fields

Field calibration in the AttoDRY800

For the S1 and S2 measurements performed in the Attocube AttoDRY800 cryostat, we use a set of electromagnets that are roughly pointed in orthogonal directions, denoted as \(x\prime \), \(y\prime \) and \(z\prime \).

All frames are shown in Extended Data Fig. 1b and are discussed above in the section ‘Notation’. The vector magnetic field B produced at the sample is linearly related to the three currents \(J_i^\prime \) applied to each of the electromagnets,

where we work in the NV⁽⁰⁾ frame (repeated indices summed). To calibrate the matrix C, we perform a series of ODMR measurements on NV centres in the facet of the diamond, in which we expect an ambient pressure environment, to mitigate any stress background in the calibration. Furthermore, this is done in situ, such that, after calibration, the cell does not move within the cryostat and the calibration remains accurate within each experimental run.

We perform the calibration measurements by finding the current configruations that produce magnetic fields close to alignment along the four NV orientations. We start by roughly orienting the magnetic field along the [111] crystal direction of diamond (roughly along the DAC axis). By varying the applied currents \(J_i^\prime \), we adjust the applied field until the ODMR spectrum exhibits two widely spaced peaks because of the Zeeman splitting of the NV⁽⁰⁾ group and two peaks at 1/3 the spacing of the outer pair. The degeneracy of the peaks of the other three groups corresponds to them all sharing an angle with cosine −1/3 with the magnetic field.

This current configuration defines \(\bfB\propto \widehatz\).

We then repeat the above alignment procedure for each of the remaining three NV groups g = 1, 2, 3 and thereby obtain reference current configurations that produce \(\bfB\propto \widehatz^(g)\).

At this point, we have an over-constrained set of linear equations to solve for the nine components of C. We use three of the measurements to solve for C, and use the fourth one to estimate the error in the alignment:

where Δν is the Zeeman splitting of the widely spaced peaks from the fourth NV group. The typical error is below 1%, which corresponds to a misaligned angle below about 0.6°. We use the calibrated C to determine the electromagnet current values required for any desired vector magnetic field B.

Field calibration in the AttoDRY2100

Calibration of sample S3 in the Attocube AttoDRY2100 cryostat follows essentially the same procedure as for samples S1 and S2, but with additional refinement. As above, we seek the transformation matrix that enables calculation of fields in the NV frame. The AttoDRY2100 has three built-in orthogonal superconducting magnets, which simplifies our task. We perform the method described above to get an initial solution for C, and then refine these values with the following automated procedure. Thirty randomly oriented 100 G fields are applied to the sample using roughly aligned currents. For each, we obtain a corresponding ODMR spectrum in which the peaks are fitted using a peak finder that is seeded with the expected positions of the peaks from the initial calibration C computed from the NV ground state Hamiltonian. We then use the measured peak locations to obtain a best-fit B from the model Hamiltonian.

We assemble the 30 applied current vectors \(J_i^\prime \) and best-fit measured field vectors B_i into 30 × 3 matrices A and F, respectively. We then compute the Moore–Penrose pseudo-inverse A⁺ and form the transfer matrix

which formally gives the optimal solution to this over-constrained linear system (here, m runs over the 30 measurements). This method allows us to refine our calibration with arbitrarily many test fields and can achieve misalignments <0.1°.

Diamond anvil miscut

Polishing (111) cut anvils is technically challenging, and deviations between the culet plane and the ideal crystal plane up to θ < 5° are to be expected. That is to say, the [111] NV axis \(\widehatz\) has a miscut angle θ with respect to the culet normal \(\widehatZ\). This is shown in Extended Data Fig. 1a, in which we define a lab frame with XYZ axes aligned so that \(\widehatZ\) is orthogonal to the culet and a rotated NV frame with xyz axes aligned conventionally with respect to the [111] NV symmetry axis.

The magnetic field B and stress tensors \(\overleftrightarrow\sigma \) extracted from NV splittings are most naturally constructed in the NV frame. Neglecting the miscut generally leads to negligibly small errors in B_Z measurements and normal stress σ_ZZ components, as these are only corrected at O(θ²). However, it generally leads to more important corrections to the transverse components of these quantities. For example, the miscut rotates the large normal ZZ stresses into the smaller shear components at O(θ¹). Thus, to correctly estimate the shear on the sample, given by the σ_*Z components in the lab frame, we must take the miscut into account.

To rotate measured stress tensors into the lab frame from the NV frame we require the rotation matrix,

where the Euler angles ϕ and θ specify the orientation and magnitude of the miscut. The final rotation ψ leaves the culet normal invariant and can be used to orient the lab XY-axes with the axes of the image plane in the CCD.

We measure the miscut angle with two independent methods. First, we perform single-crystal X-ray diffraction on the diamond and resolve the tilting of the diamond lattice with respect to the culet. Second, we load the DAC with N₂ pressure medium and pressurize to 9 GPa. The negligible ability of N₂ to sustain shear⁴³ means that σ_XZ, σ_YZ ≈ 0 across the culet. We can, therefore, evaluate σ^Lab = Rσ^NVR^T and solve for the rotation matrix R such that the \(\sigma _XZ^\mathrmLab\approx \sigma _YZ^\mathrmLab\approx 0\).

The N₂ method gives a miscut magnitude of θ = 3.5(6)° and an azimuthal orientation of ϕ = −5(6)° for S1. XRD gives a miscut magnitude of θ = 3.02°, but disambiguating the azimuthal angle proved difficult. We chose to use the miscut parameters derived from the N₂ method, noting that the XRD method is also more susceptible to small misalignments between the culet plane and the detector owing to imperfections in the XRD sample mounting. Regardless, we note that the two miscut angles agree within error, and these results are within the specifications of the diamond polishing method. We apply this rotation to measured stress tensors to extract the lab frame stress.

The effect of this rotation is rather small. As mentioned above, the adjustment is first order for shears, and results in a correction of about 1 GPa, and second order for normal stress, with a correction of about 0.1 GPa. We have checked that arbitrary (incorrect) miscut corrections with θ < 5° do not qualitatively change the shape of the phase diagrams and shear–superconductivity correlations presented in the main text, establishing that miscut corrections are minimal.

ODMR peak fitting

A single wide-field image has 10⁵ pixels, each containing an ODMR spectrum. We use a highly parallelized fitting strategy to extract magnetic field and stress parameters. Each ODMR spectrum has up to eight NV resonances: a pair from each of the four NV groups. Rather than directly working with the two peak frequencies for each group g, it is more convenient to work with their average frequency D_g and their half-difference E_g. For example, extracting B_z requires finding the E₀ splitting, whereas a stress tensor solution generally necessitates finding the D_g for all groups g. These tasks, naturally, require that we correctly identify which resonances come from which NV group, and fit their frequencies accurately. This is complicated by several factors: stress and field gradients in the sample broaden and weaken the resonances; stress-induced symmetry breaking of the point group of the NV centre reduces peak height; and the laser and microwave polarizations are less effective at polarizing certain NV orientations. We use a custom-built MATLAB program to address these challenges. The operational procedures are as follows:

1. 1.

   In each pixel, we deploy the standard MATLAB peak-finding routine to roughly identify the position of all NV resonances after first denoising the data with a Savitzky–Golay filter.

2. 2.

   We then pair the ±1 resonances of a given NV group by their proximity to a regionally seeded value of either D_g or E_g using a maximum bipartite matching algorithm. As an example, let us focus on the [111] group. If a field of 100 G is applied along the [111] NV, we expect the resonances to have an E₀ ≈ 280 MHz. Because B varies at most by about 15%, and because the other NV groups will experience a splitting of about 1/3 of what is experienced by the [111] NV, the pair of peaks that has an E splitting nearest to this seed is reliably the [111] NV.

   The limited varisation of B makes E_g seeding well suited for La₃Ni₂O₇ mapping; BSCCO (Extended Data Fig. 2) is better suited to D_g seeding as in that sample B varies greatly, whereas stress (which controls D_g) is relatively uniform.

3. 3.

   We impose constraints on peak height, width, D_g and/or E_g to further enhance the accuracy in partnering resonances. Resonances from the same NV group are expected to have similar height and width by symmetry constraints.

4. 4.

   In <1% of spectra within the sample chamber, this partnering fails because of the poor signal-to-noise ratio of one of the peaks in the pair. Thankfully, knowledge of either D_g or E_g allows the fit to proceed even if only one resonance is available. If we seek to fit E_g, we may extract D_g from a dataset at a different temperature or magnetic field (which influences D_g trivially), in which the formerly obscured peak is recovered. Likewise, fits of D_g can use E_g values obtained at a different temperature.

5. 5.

   Once peaks are identified, resonance values are further refined by quadratic interpolation around the frequency points identified in the peak-finding step to mitigate sampling error (we sample our ODMR every 2 MHz).

6. 6.

   The processes in steps 1–5 are repeated in parallel for all pixels, allowing for a mapping of D_g or E_g for the desired NV resonance.

Our automated method generally completes in <10 s.

From here, several further refinements are possible:

1. 1.

   Our software contains a customized version of the QDMLab GPU fitting code developed in ref. ⁶⁶. The method above identifies only the peak of the resonances, so for benchmarking, we use the QDMLab-based module for curve-fitting the whole resonances with Gaussian line-shapes, seeded by the results of the peak-finding procedure. In general, re-fitting with QDMLab does not significantly change the results of the peak-finding method for any of the data presented here.

2. 2.

   In the case of high stress and low temperature, we often observe the inversion of contrast in the non-[111] resonances (see section ‘Positive contrast’). In some high-stress datasets, we find several pixels in which these inverted resonances are the only ones present, and that both [111] resonances are weak. As we typically rely on the [111] resonances to determine B_z, this is a challenge. Knowing the projection of the magnetic field along any of the non-[111] NV subgroups is sufficient to extract B_z, given knowledge of the local orientation of the field. We can determine this orientation locally; adjacent pixels in which both [111] and non-[111] resonances are visible allow us to determine the relative projection of B onto these groups. Generally, B is essentially parallel to H (which is along \(\widehatz\)) owing to the rather small diamagnetism (1 > B/H > 0.97) in these regions.

3. 3.

   Even after all these steps, we find that of order tens of pixels fail to automatically fit in typical magnetic B/H maps. We manually fill in these rare pixels by assuming the B_z is harmonic in the 2D plane and visually checking that the ODMR spectrum has plausible, if weak, peaks at those positions.

4. 4.

   For maps of stress, when the automated peak fitting procedure fails, we can often resolve peaks manually by reseeding the fitting procedure using rough peak positions determined by continuity with nearby pixels. Again, when we need to do manual reseeding, we confirm peak assignments visually within each ODMR spectrum.

Positive contrast

Generally, NV resonances generate a fluorescence dip in ODMR spectra, referred to as negative contrast. However, we frequently observe inverted positive contrast peaks (for example, in Fig. 1e) at low temperature and high pressure (\(\bar\sigma _ZZ\gtrsim 15\,\mathrmGPa\)). We can follow the inversion of peaks as a function of decreasing temperature and increasing stress and field, and thereby identify them with their parent negative contrast resonances. At higher stress points, we make use of these resonances in extracting local magnetism and stress.

The phenomenon of positive contrast in pressurized NVs has been reported by multiple groups^21,60, without a comprehensive explanation for their origin. In a forthcoming paper⁶⁷, we provide a microscopic model for this intriguing and useful phenomenon.

Constructing the local stress tensor from NV measurements

In the absence of any applied magnetic field, the average shift D_g and splitting E_g of the ODMR peaks of the NV centres in group g are determined by the NV symmetry-preserving and symmetry-breaking stresses, respectively. With these eight quantities in hand, it is possible to reconstruct the local stress tensor \(\overleftrightarrow\sigma \), whose six components are over-constrained. In practice, at zero magnetic field, the ODMR peaks because of the various groups g overlap, and it is impossible to accurately resolve the peak positions and/or assign them to the correct groups.

The standard method for stress tensor extraction thus makes use of two kinds of ODMR measurements at a given point¹⁸, both of which exploit carefully aligned applied magnetic fields. First, the ‘D’ measurements are performed by applying a longitudinal magnetic field along the NV^(g) axis to split the resonances of group g from those of the other three groups.

Second, the ‘E’ measurements require a magnetic field that is perpendicular to the relevant NV^(g) axis. This now splits away the three other groups and allows resolution of the two central peaks of interest.

However, there are confounding effects. To understand these, we recapitulate the effective stress Hamiltonian governing the NV spin sub-levels, written in the local NV frame^18,60,61:

where

and a₁, a₂, b, c are the stress susceptibilities of the NV centre^18,62. The Π_z term shifts the two ODMR peaks corresponding to the 0 ↔ ±1 transitions together, directly contributing to D, whereas the Π_x, Π_y parts of the stress couple −1 ↔ +1 and contribute to the relative splitting E.

Extracting the stress contribution Π_z from the ‘D’ measurement is relatively straightforward, as \(D=\varPi _z+D_\mathrmgs(T)+\frac(\gamma B_\perp )^2D_\mathrmgs(T)\), including both zero-field splitting D_gs and leading order magnetic effects. The temperature-dependent zero-field splitting D_gs(T) is well characterized⁶⁸ and easy to subtract. Moreover, as it is about 3 GHz, the contribution from transverse magnetic fields γB_⊥ < 3 MHz (which are generated from field misalignment) is negligible. As such, Π_z is robust against magnetic fields and may be reliably extracted.

By contrast, the ‘E’ measurements are much more challenging. First, the sizable transverse magnetic fields (about 100 G) required to split away the non-group g peaks also tend to mix the NV spin sub-levels, which reduces the ODMR contrast. This is compounded by a similar reduction of contrast because of NV symmetry-breaking stress, which is non-trivial in the pressure range of this work. Even when the group g peaks can be resolved, E is a complicated function of Π_x, Π_y, B and contributions due to inhomogeneous internal electric fields that produce an ambient (no stress) effective splitting \(\mathcalE_\) (ref. ¹⁸). The \(\mathcalE_\) contribution combines with stress (Π_x, Π_y) by an elliptic integral; this composite quantity combines in quadrature with B_z and trigonometrically with (B_x, B_y). Extracting the stress contributions from the E measurement is highly nonlinear and is very sensitive to slight misalignments in the applied magnetic field. Finally, we note that there is further experimental uncertainty in the stress susceptibilities (b, c) by which (Π_x, Π_y) couple because of the aforementioned confounding effects of electric and magnetic fields (which demonstrably affect their calibration¹⁸).

In short, although wide-field D_g maps for each of the four NV groups can be extracted reliably, the E_g maps are generally much more challenging because of reduced ODMR contrast, tight magnetic field orientation tolerances and uncertainties in the relevant susceptibilities. We can take advantage of two simplifications even in the absence of E_g data. First, only three components of the stress tensor are physically relevant to the sample—the traction σ*_Z—and we need not reconstruct the full stress tensor. Second, the stress field must satisfy the constraints of mechanical equilibrium. In particular, the stress at the surface of the diamond anvil is dominantly σ_ZZ and even in the presence of local surface shears, we expect the internal shear components of the stress to be very small. Building on this assumption, we define a method to reconstruct the normal stress σ_ZZ and shear τ from the four D_g values alone.

To test the approximate 4D method, we benchmark it against full reconstruction using a standard 4D + 2E method at a collection of sampled points at which the relevant E peaks were readily resolved. We find that the two methods produce linearly correlated results for the shear traction τ, with a slope that characterizes the error of the accepted assumptions, which we use for rescaling all shear stresses obtained by the 4D method.

The 4D method

The 4D method is motivated by the observation that the stress field throughout a linear elastic medium, such as the diamond anvil, is determined by the three components of the surface traction \(\boldsymbol\tau ,\sigma _ZZ\). Finite-element simulations of diamond suggest that the components σ_XX−σ_YY and σ_XY generated by these tractions are small. Thus, we make the simplifying assumption that these two components of the stress tensor vanish—and thus the four measured D_g values—are sufficient to fix the four remaining components of \(\overleftrightarrow\sigma \).

To recombine the four measured D_g values into an estimate of the stress tensor in the lab frame, it is convenient to re-express the NV stress susceptibility in terms of tensorial quantities:

Here, the D_g susceptibility Y^(g) is a symmetric rank-2 tensor, which is invariant under the point group of NV group g and \(\widehatz^(g)\) is the z-axis of group g. We take the susceptibilities a₁, a₂ = 2π × 4.86(2), −3.7(2) MHz GPa (ref. ⁶⁹). Given the four measured values of D_g at a point, and the vanishing internal shear assumption, it is straightforward to solve for the stress tensor σ_IJ in the lab frame.

We note that several physically important components of the stress tensor are completely determined even without the assumption of vanishing internal shear. These include the pressure \(p=\frac13\sigma _ii\) and the normal stress in any of the NV-axis directions, notably σ_zz, which, up to a small miscut correction, is the normal component of the stress on the sample chamber. In particular, noting that \(\sum _g\widehatz^(g)\widehatz^(g)=4\mathbbI/3\), we obtain p = ∑D_g/(12a₁).

Benchmarking with the 4D + 2E method

We validate the 4D method by comparing it with the standard stress tensor extraction strategy, using, apart from the 4D values, 2E values to fully constrain the last two components of \(\overleftrightarrow\sigma \). We note that the normal component of the boundary traction σ_ZZ is the same in both methods (up to small miscut corrections). In Extended Data Fig. 8a, we randomly sample points on S1 at which the ODMR spectrum for the E measurements showed clearly resolvable peaks and calculate the shear traction τ using both methods. We observe a strong linear correlation, τ_4D ≈ 1.8τ_4D+2E (Extended Data Fig. 8b). The fact that the slope is not one and that there is some offset may be attributed to both uncertainties in the standard method (several GPa) and uncertainties in the susceptibilities, especially b, c (ref. ¹⁸). We consider τ_4D+2E the true shear stress and rescale all τ obtained by the 4D method by this factor of 1.8.

The shear traction τ in the 4D method is independent of a₁ and thus depends inversely on the susceptibility a₂. Changing this susceptibility uniformly rescales the traction components. Various groups have reported different values of a₂, such as a₂ = −6.5 MHz GPa⁻¹ (ref. ⁷⁰) and a₂ = −2.51 MHz GPa⁻¹ (ref. ⁶¹).

We use the susceptibilities calculated in ref. ⁶⁹. Future work is needed to clarify the true value of this susceptibility, but adjusting our results will simply require the multiplication of some constant to all shear values in our work (and a similarly simple correction to σ_ZZ) and not qualitatively change the shape of the phase diagrams we extract.

Extraction of normal component

For low-stress datasets in S1 (\(\bar\sigma _ZZ < 17\,\mathrmGPa\)), we did not take separate measurements of D₁₋₃. However, when H is applied along \(\widehatz\), we see that the non-[111] resonances are strongly overlapping, as they are degenerate under this magnetic field. This indicates that D₁, D₂ and D₃ are approximately equal and given by the average of the two visible peaks. This permits the construction of σ_ZZ.

Pixel registration

As pressure varies in a DAC, the sample may slightly change shape and/or position because of the pressure gradient. Furthermore, temperature cycling can lead to subtle translations in the sample due to piezo drift. Thus, it is crucial to register the multimodal information to a common set of pixel coordinates.

Here, we show how to register maps from different pressure points and from different probes.

To register between different probes at the same pressure (NV fluorescence, EDX and white light), we manually select key features such as the ruby, edge of the sample and the leads as control points (no fewer than 10 features) and use these to calculate best-fit projective transformations between the images.

To register between different pressure points, we rely on white-light images, as these have the most visible structure to correlate. In this case, we select invariant key features using the Speeded-Up Robust Features algorithm in MATLAB. We then use these features as control points to calculate the best-fit affine transformation between different pressure points. We restrict to affine transformations, rather than the more general projective transformations permitted between probes, as the white-light imaging optics are fixed across all pressures.

Constructing the pixel-based phase diagram

In the main text, we present phase diagrams in σ_ZZ–τ–T space created from aggregated NV data across all samples. At a high level, we begin by selecting all the pixels that fall within a defined sample-spanning ‘region of interest’ across stress and temperature for each sample. Then, we define a vector of stresses, temperature and B/H for each pixel by the collation of the aligned magnetic and stress maps. Finally, we plot these points in the respective coordinate spaces to generate the phase diagrams.

Now, we will detail the specific procedures for making the two 2D phase diagrams and the full 3D phase diagram, as well as the stoichiometric phase diagram.

σ
_ZZ compared with T

For S1, S2 and S3, maps of σ_ZZ and B/H are constructed at all stress and temperature points.

For each sample, these maps are aligned with respect to each other, as previously discussed. The goal of this phase diagram is to depict the effect of σ_ZZ variation on B/H, so we mask out those parts of each sample that never conclusively superconduct at any point in our dataset. The samples contain substantial regions that never superconduct, for reasons that are evidently beyond the scope of stress (for example, chemical inhomogeneity in S3 and possible oxygen inhomogeneity or layer stacking disorder in S2). We could, in principle, generate this phase diagram by collecting pixels from the whole sample region, but in practice this generates a qualitatively identical phase diagram with a worse signal-to-noise ratio because of the integration over tens of thousands of ‘dead’ pixels. Similarly, we mask out sample regions with high shear (above about 2 GPa).

Doing this across all datasets yields pixels of the order of one million. We plot these points in (σ_ZZ, T) space with their colour controlled by B/H and find excellent coverage of phase space—emphasizing how inhomogeneity can be leveraged as a metrological resource. For denoising purposes, we perform k-means clustering of these points (k = 250) after first temporarily rescaling the temperature axis (by a factor of 1/20, so as to roughly equate the typical spacing of points along both axes before clustering). We set a superconductivity threshold, B/H = 0.98, which used to define the diamagnetic suppression, below which we consider a local region to exhibit superconductivity; specifically, we use a sigmoidal low-passer (k = 20) centred around this value. We note that we set B/H to unity for all regions with B/H > 1; these regions are non-superconducting and correspond, for example, to regions of magnetic-field enhancement. Using these clustered points, we then interpolate along both axes to fill the phase diagram. We emphasize that our results are relatively insensitive to choices of k or axis rescaling.

σ
_ZZ compared with τ

Similarly, we generate vectors from our maps of σ_ZZ, τ and B/H. We restrict ourselves to data with T ≤ 35 K: as in the previous section, the diagram that integrates over all temperatures is qualitatively identical, but with a worse signal-to-noise ratio. We generally lack shear data for low-stress datasets (\(\bar\sigma _ZZ\) less than 16 GPa, 10 GPa and 3 GPa for S1, S2 and S3, respectively), but these datasets have very weak to no magnetic signatures, so in practice their absence does not influence the constructed superconducting dome. We perform k-means clustering (k = 85) and use the same sigmoidal low-passer, with no axis rescaling necessary. We linearly interpolate to fill the phase diagram.

Three-dimensional σ
_ZZ –τ–T phase diagram

In this plot, we eschew B/H colouration for visual simplicity and instead construct a phase boundary by plotting T_c for each (σ_ZZ, τ) coordinate. To do this, we take the list of pixels from the previous step and bin them into (σ_ZZ, τ) cells of side length 2 GPa by 1 GPa. For each cell, we fit a sigmoidal function \(B/H=a+(1-a)\left(1-\frac11+\exp (k(T-T_\rmmid))\right)\) to extract the superconducting-normal crossover T_c = T_mid + 1/k. In the main text, we define T_c as the onset of superconductivity in transport and diamagnetism, which we model here as one characteristic sigmoid width (1/k) above the crossover temperature. These T_c values are splined together along both axes, and then the spline curves are linearly interpolated to generate the full surface. We note that this interpolation is done down to 0 on all axes, which requires extrapolation. However, we emphasize that this 3D phase diagram is in excellent accordance with the two 2D phase diagrams, although it was generated in a different way.

La:Ni compared with T

In Fig. 5c, we use EDX to determine the correlation of the La:Ni ratio with diamagnetism. This was done by binning pixels by stoichiometry intervals 0.2 wide, and calculating the average B/H of the bin. Here, we extend this analysis across the temperature axis; we plot a 2D La:Ni compared with T phase diagram in Extended Data Fig. 5e. Similar to the other 2D phase diagrams, we linearly interpolate the binned points in this phase space. We use a sigmoidal low-passer threshold of B/H ≈ 0.996, reflecting the broadly weaker diamagnetism of S3. As expected, this phase diagram reproduces the two qualitative features observed in the main text: (1) T_c ≈ 60 K and (2) an optimal La/Ni ratio of about 1.5, but also provides more quantitative information in the temperature–stoichiometry plane.

Estimating superconducting volume fraction from stray field measurements

Our local measurements of B/H provide an alternate route to extract the superconducting volume fraction of the sample.

Our NV centres (embedded in the anvil culet) are about 500 nm from the sample, which is significantly smaller than the sample dimensions (about 10–100 μm). Thus, we treat the NVs as measuring the stray field at the surface of the sample. To extract the superconducting volume fraction f from our experimentally measured average diamagnetic suppression, ⟨B/H⟩, we work with the following simple model. Assume that a fraction f of the sample is a perfect superconductor with χ_SC = −1, and the remaining fraction of the material has a paramagnetic response, that is, χ_PM ≈ 0. The effective magnetic susceptibility of the sample is then given by

As our NVs effectively measure the stray field at the surface of the sample, we can approximate χ_eff ≈ ⟨B/H⟩ − 1 and therefore \(|\chi _\mathrmeff|\approx 1-\langle B/H\rangle \). Thus, for our three samples, we obtain a superconducting volume fraction of 3%, 0.3% and 0.5% for S1, S2 and S3, respectively. We may compare these volume fractions to those extracted in the literature using conventional methods involving a.c. or d.c. susceptibility measurements. Broadly, the highest volume fractions (>40%) are seen for the most hydrostatic pressure media, such as helium⁹. This aligns with our finding that shear quenches superconductivity; more hydrostatic pressure media are less able to support shear.

Wohlleben effect

S1 exhibits intriguing magnetic behaviour: strong regions of paramagnetism that seem spatially separated from superconducting diamagnetism, but nevertheless vanish simultaneously with it above T_c. Although we expect, as in BSCCO, to see a halo of B/H > 1 surrounding diamagnetic regions (Extended Data Fig. 2c) as a result of field expulsion, the paramagnetism in S1 is unusual in its asymmetry (appearing predominantly on the left edge of the sample) and in its relative disconnection from the superconducting regions. In the main text, we connect this behaviour to the Wohlleben effect, which we will explore further here. This effect has been extensively observed and documented in granular superconductors^46,48. For example, some models rely on Josephson coupling between discrete superconducting grains, leading to spontaneous currents, which generate the paramagnetism⁷¹. The fact that granularity seems to be crucial for observing this type of paramagnetic enhancement seems to point in the right direction for our La₃Ni₂O₇ samples, which exhibit strong heterogeneity at the μm scale.

Finally, although the Wohlleben effect has most often been examined under field-cooled protocols in previous works, some studies report a small paramagnetic signal in zero-field-cooled measurements a few kelvin below T_c (ref. ⁷²). By contrast, we observe a clear ZFC paramagnetic response across the entire sub-T_c range, with a magnitude that is approximately proportional to the diamagnetic signal.

At a high level, some version of the Wohlleben effect seems like the most reasonable explanation for the data we see, given its robustness across pressure points and samples and its appearance in our granular La₃Ni₂O₇ system and not in other clean systems that we have studied such as BSCCO.

Gradients of stress

Another way hydrostaticity may be violated is through stress gradients. Here, we ascertain if these stress gradients correlate with diamagnetic textures. In Extended Data Fig. 9a, we reproduce the diamagnetic map of S1 at \(\bar\sigma _ZZ=21\,\mathrmGPa\) and 20 K. In Extended Data Fig. 9b,c, we plot the gradients in σ_ZZ and τ, respectively. We then draw a white oval around the regions of maximal σ_ZZ and a yellow oval around the regions of τ, and also plot these ovals on Extended Data Fig. 9a. We observe that large gradients in either σ_ZZ or τ do not suppress superconductivity. So, it seems that stress gradients do not correlate with the absence of superconductivity.