Observation of Caroli-de Gennes-Matricon Vortex States in YBa$_2$Cu$_3$O$_{7-\delta}$

TL;DR

This study combines experimental tunneling measurements and theoretical modeling to identify vortex core states in YBa2Cu3O7−δ, suggesting that copper oxide superconductivity is conventional but masked by minority carriers.

Contribution

It provides the first positive identification of vortex states in copper oxide superconductors through combined experimental and theoretical approaches.

Findings

Vortex core states are present but weak in YBa2Cu3O7−δ.

Variations in vortex spectra are explained by vortex interactions and disorder.

Superconductivity in copper oxides may be conventional with minority carriers.

Abstract

The copper oxides present the highest superconducting temperature and properties at odds with other compounds, suggestive of a fundamentally different superconductivity. In particular, the Abrikosov vortices fail to exhibit localized states expected and observed in all clean superconductors. We have explored the possibility that the elusive vortex-core signatures are actually present but weak. Combining local tunneling measurements with large-scale theoretical modeling, we positively identify the vortex states in YBa$_2$Cu$_3$O$_{7-\delta}$. We explain their spectrum and the observed variations thereof from one vortex to the next by considering the effects of nearby vortices and disorder in the vortex lattice. We argue that the superconductivity of copper oxides is conventional, but the spectroscopic signature does not look so because the superconducting carriers are a minority.

Figures

Full text

Type-II superconductors immersed in a magnetic field let quantized flux tubes perforate them: the Abrikosov vortices. This remarkable property underlies and often limits many applications of superconductors. In 1964, Caroli, de Gennes, and Matricon used the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity to predict that vortices in type-II superconductors host a collection of localized electrons bound to their core Caroli et al. (1964). The direct observation of these localized states 25 years later by scanning tunneling spectroscopy (STS) is a spectacular verification of the BCS theory Hess et al. (1989). The formation of vortex-core bound states is an immediate consequence of the superconducting condensate being composed of electron pairs, while excitations in the vortex, being unpaired, have a different topology. Core states are, therefore, a robust property of superconductors, like edge states in topological insulators, irrespective of the origin and symmetry of the force that glues the electrons into pairs. In spectroscopy, they appear in the clean limit $\ell\gg\xi$ as a zero-bias peak in the local density of states (LDOS) at the vortex center, where $\ell$ and $\xi$ are the electron mean free path and superconducting coherence length, respectively, or as a structureless LDOS in the dirty limit $\ell\lesssim\xi$ Renner et al. (1991). Next to NbSe2 Hess et al. (1989); Suderow et al. (2014), the core states were seen by STS in several superconducting materials DeWilde et al. (1997); Eskildsen et al. (2002); Nishimori et al. (2004); Guillamón et al. (2008); Zhou et al. (2013); Du et al. (2015), including the pnictides, which are believed to host unconventional pairing Yin et al. (2009); Shan et al. (2011); Song et al. (2011); Hanaguri et al. (2012); Fan et al. (2015). The high-$T_{c}$ cuprates stand out as the only materials in which the vortex-core states have been looked for but not found. In YBa2Cu3O7-δ (Y123), discrete finite-energy structures initially believed to be vortex states Maggio-Aprile et al. (1995); Shibata et al. (2003); Shibata-2010 were recently shown to be unrelated to vortices Bruér et al. (2016). In Bi2Sr2CaCu2O8+δ* (Bi2212), the vortex cores present no trace of a robust zero-bias peak, but instead very weak finite-energy features apparently related to a charge-density wave order Renner et al. (1998); Hoogenboom et al. (2000a); Pan et al. (2000); Matsuba et al. (2003); *Matsuba-2007; Levy et al. (2005); Yoshizawa et al. (2013).

The absence of vortex states in cuprates is challenging the existing theories. Because these states are topological they are robust Wang and MacDonald (1995); Franz and Tešanović (1998), and they survive modifications of the BCS theory like strong-coupling extensions Berthod (2015) that do not change the nature of the condensate. To explain the cuprate vortex phenomenology, one needs either to leave BCS theory Arovas et al. (1997); Himeda et al. (1997); Andersen et al. (2000); Kishine et al. (2001); Berthod and Giovannini (2001), or to extend it by including additional order parameters that condense inside the vortex cores and gap out the zero-bias peak Zhu and Ting (2001); Maśka and Mierzejewski (2003); Fogelström (2011). To date, none of these approaches has given a satisfactory account of the phenomenology observed by STS. The discovery that the low-energy structures in Y123 do not belong to vortices Bruér et al. (2016) suggests that these theoretical efforts have been misguided.

The electronic structure of the cuprate high-$T_{c}$ superconductors (HTS) is notoriously complex, as manifested in a rich phase diagram. This complexity reveals a competition of different effective interactions, from which a variety of individual and collective modes emerges progressively as the temperature is lowered towards the ground state in which the system freezes at absolute zero, and which continues to keep the secret of the most stable superconductivity ever observed. It is generally believed that the phenomena taking place close to the Fermi surface—charge, spin, pairing orders, and their fluctuations—all derive from a single band Zhang and Rice (1988) or a small subset of bands Emery (1987) in the CuO2 layer(s). Consistently, the interpretations of STS spectra Fischer et al. (2007) have postulated that all electrons contributing to the measured LDOS are excited out of the superconducting condensate, in agreement with Leggett’s theorem Leggett (1998). On the other hand, it is well known that for all cuprates, at any doping, the superfluid density remains much smaller than the electron density Uemura et al. (1989); Bernhard et al. (1995); Božović et al. (2016). Our recent high-resolution STS experiments on Y123 have also revealed that only a fraction of the signal recorded on the sample surface is of superconducting origin Bruér et al. (2016). Early specific-heat measurements have given a similar hint Junod et al. (1999). We are therefore lead to a new paradigm, in which the low-energy electronic state of the HTS involves a minority superconducting channel in parallel with nonsuperconducting majority charges guilty for the pseudogap and the associated orders. Here we show that the minority carriers are fairly conventional in the superconducting state, showing Caroli–de Gennes–Matricon states in the vortex cores as predicted by the BCS theory for $d$-wave superconductors.

STS measurements with a normal metal tip do not discriminate the superconducting (SC) and nonsuperconducting (NSC) channels and collect electrons from each. Our working hypothesis is that the tunneling conductances originating from the SC and NSC channels are additive. An earlier report in which vortex-induced changes were tracked in Bi2212 rested on a different scenario, namely, the vortices would suppress locally the condensate and reveal a competing magnetic order Hoffman et al. (2002). As the NSC is not known, a simple subtraction to reveal the SC is not feasible. Yet, inhomogeneities of the SC like those induced by vortices can be singled out by subtracting the tunneling conductance outside vortices from that inside. If the formation of a vortex leaves the NSC unchanged, this procedure eliminates the NSC and permits a comparison of the SC inhomogeneities with the BCS theory. Figure 1(a) shows the LDOS predicted by the BCS theory at and near the center of an isolated vortex in a two-dimensional superconductor with electronic structure similar to that of Y123 111See Supplemental Material [url] for a presentation of the model and calculation methods, which includes Refs. Jorgensen et al., 1990; Atkinson, 2009; Schabel et al., 1998a; *Schabel-1998b; Ichioka et al., 1996; Torquato and Stillinger, 2003, 2003; Weismann et al., 2009.. The prominent feature is the zero-energy peak localized at the core center. With increasing distance from the core, the peak is suppressed with or without a splitting depending on the direction. Note that this LDOS anisotropy is unrelated to the $d$-wave gap anisotropy. In the quantum regime $k_{\mathrm{F}}\xi\sim 1$ relevant for Y123, the vortex size is comparable with the Fermi wavelength and the Fermi-surface anisotropy determines the vortex structure Uranga et al. (2016); Berthod (2016). In the Supplementary Material, Fig. S8 indeed shows that the zero-bias LDOS is locked to the crystal rather than gap-node directions. For a meaningful comparison with experiment, the LDOS far from the vortex must be subtracted, as done in Fig. 1(b).

Figure 2(a) shows a $90\times 90~{}\mathrm{nm}^{2}$ area on the surface of Y123, where 19 inhomogeneities can be identified as vortex cores. Details about the sample preparation and measurement technique were reported in Ref. Bruér et al., 2016. Because of the large NSC, the contrast due to vortices is weak; it is maximized by mapping the ratio of the STS tunneling conductance at 5 and 17 mV bias [Fig. 2(a)]. Local increases of the zero-bias conductance are also seen in the raw data [Fig. 2(b)], and correlate well with the vortex positions determined by the best contrast. In order to remove the NSC, we delineate a small region in-between vortices, calculate the average spectrum in this region [inset of Fig. 2(b)], and subtract this average from all spectra in the map. After subtraction, the spectral traces show the expected vortex signature with a maximum at zero bias in the cores. This is demonstrated in Figs. 2(c) and 2(d) with two traces running from one vortex core along the two directions shown in Fig. 2(a). The peak developing locally at zero bias after subtracting the same background from the spectra of each trace clearly shows there is a larger local density at low energy near the vortex cores. It is not an artifact of the subtraction. Similar results are found in all vortices.

A comparison of Figs. 2(c) and 2(d) with Fig. 1(b) reveals evident similarities, but also differences. Similarities include the zero-bias peak, which has its maximum at the core center and is suppressed with increasing distance from the center, the absence of superconducting coherence peaks in the core leading to symmetric dips at $\pm 17$ meV in both experiment and theory, a striking spatial anisotropy, with the peak extending farther along path no2 than along path no1, reminiscent of the different decay lengths observed in the theory between directions 10 and 11. Among the differences, one notices the central peak being taller in Fig. 1(b) than in the measurements, and the spectrum appearing locally reinforced or split at intermediate distances along path no1 and no2, respectively, while Fig. 1(b) shows a monotonic evolution. We now demonstrate that these differences can be explained by considering that (i) the vortices are not isolated in the experiment, (ii) the relative orientation of the vortex and crystal lattices influences the LDOS, and (iii) the vortex lattice is disordered. Incorporating (i) in the theory reduces drastically the calculated peak height; (ii) means that the LDOS anisotropy depends upon the positions of nearby vortices; finally, (iii) implies that each vortex sits in a specific local environment and presents spectra different from its neighbors.

Figure 3 displays a series of spectral traces along various paths connecting either nearest-neighbor [Fig. 3(a)] or next-nearest neighbor [Fig. 3(b)] vortices. The notion of nearest- and next-nearest neighbor refers to a local fourfold coordination generally observed among the vortices, despite the long-range disorder in their arrangement. A trend is systematically observed: along a line connecting nearest-neighbor vortices, the zero-bias peak remains visible along the whole path, while it disappears along paths joining next-nearest neighbors. Considering the peak anisotropy as it is predicted by theory (Fig. 1), this trend suggests that the locally fourfold-coordinated vortex lattice tends to align along the crystalline axes, as we will confirm by a detailed modeling. In the absence of atomic resolution imaging, we infer the lattice orientation based on optical images of twin boundaries. Another lesson of Fig. 3 is that all vortices, although similar, are different: some show a single peak, others show a split peak; the height of these peaks is also varying. This variability reflects disorder in the vortex positions, resulting in irregular distributions of supercurrents around each core.

We have undertaken large-scale simulations of the LDOS in disordered vortex configurations, in order to study how this modifies the core spectra with respect to the isolated vortex, and thus better understand the observations made by STS. Thanks to a new approach described in Ref. Berthod, 2016, we are able to compute the LDOS in a disordered vortex lattice with the same accuracy as in the isolated core. All relevant details are given as Supplemental Material Note (1), and we only briefly review here the key ingredients. The structure of the vortices in a finite field is deduced from the self-consistent solution in the ideal vortex lattice. We find that the LDOS changes dramatically with varying the orientation of the vortex lattice relative to the crystal axes. The system size required in order to compute the LDOS with our target resolution (3 meV) contains half a million unit cells and extends well beyond the area of our STS experiments, which contains roughly 54 500 Cu sites. Inside the STS field of view, we locate vortices at the positions indicated in Fig. 4(a). Outside, we generate vortex positions randomly, however, constraining the intervortex distance to be at least 16 nm. The resulting configurations show no orientational order, but short-range coordination similar to what is seen in the STS image. We do not imply that the actual vortex distribution outside the field of view has no orientational order—in fact, it probably has some Maggio-Aprile et al. (1995)—but the available data prevent us from inferring such an order. We also vary the orientation of the crystal lattice with respect to the vortices. For each of 600 generated configurations, we calculate the LDOS along the paths no6 and no8 of Fig. 3 and compare with the experimental traces. With the configuration giving the smallest difference with experiment on these two traces, we recalculate the LDOS on the whole domain covered by STS with a resolution of 1 nm and deduce the theoretical map and traces shown Fig. 4(b).

The simulation confirms that quasiparticle scattering off nearby vortices reduces the central peak in each vortex, that the LDOS is different in all cores, and that it depends on the configuration of the vortices outside the STS field of view. We also find that the agreement with experiment is systematically better if the orientation of the microscopic lattice is such that the traces no6 and no8 are close to a nodal direction, in agreement with the twin-boundary directions. Although our search for a good vortex configuration has focused on two traces connecting next-nearest neighbor vortices, the resulting model reproduces the difference between these traces and those connecting nearest-neighbor vortices. It is also striking that the model correctly predicts the contrast of the STS image and the apparent size of the vortex cores without any further adjustments. Figure 4(b) presents less spectral variations from core to core than Fig. 4(a), but this appears to be a compromise by which this configuration achieves the best overall agreement. Other vortex configurations that compete closely do show variations comparable with experiment, including split peaks in some of the vortices (Fig. S7). We emphasize that our procedure does not deliver the best fit, which would require us to optimize systematically the vortex positions rather than trying random configurations. The vortex positions within the field of view should be optimized as well. Figure 4(b) indeed reveals that the LDOS maxima are in general displaced with respect to the points where the order parameter vanishes: the LDOS gets polarized by asymmetries in the supercurrents Berthod (2013, 2016). Without any doubt, vortex configurations could be found that improve the agreement in Fig. 4, but we feel that such a costly optimization is unlikely to reveal new physics. Interestingly, the simulation presents vortex cores that appear to be split, e.g., at the beginning of path no6 (see also Fig. S8). Reference Hoogenboom et al., 2000b reported a similar observation in Bi2212, which was ascribed to interaction with pinning centers. In Fig. 4(b) the splitting is merely a LDOS distortion due to an irregular distribution of the supercurrent.

In a clean BCS superconductor, the zero-temperature superfluid density $n_{s0}$ equals the total electron density $n$ Leggett (1998). Estimates based on the penetration depth show instead that $n_{s0}/n=23\%$–$34\%$ in Y123 222At an optimal hole doping of 0.16, the nominal carrier density of the bilayer Y123 is $n=2\times 0.84/(abc)=9.7~{}\mathrm{nm}^{-3}$ with the unit-cell parameters $(a,b,c)=(3.82,3.89,11.68)$ Å. With an effective mass of typically 2–3 times the electron mass Padilla et al. (2005) and a zero-temperature penetration depth $\lambda=1600$ Å Basov et al. (1995), the estimated superfluid density $n_{s}=m^{*}/(\mu_{0}e^{2}\lambda^{2})$ is 2.2–3.3 nm*-3*.. The analysis of STS data yields similar relations between the SC and NSC, with some model dependence, 14% in Ref. Bruér et al., 2016, 25% in the present work Note (1). In a one-channel picture, the property $n_{s0}\ll n$ requires large pairing fluctuations that would invalidate a mean-field description. Our results show that the mean-field theory works well in vortices, though. Another possibility would be that superconductivity emerges in a population of low-energy quasiparticles carrying only a small fraction of the spectral weight: in the RVB theory, the quasiparticle weight behaves as $2x/(1+x)$ as a function of doping $x$ Anderson et al. (2004), as observed in optical data of underdoped cuprates Mirzaei et al. (2013), which yields a value close to 25% at optimal doping. In this scenario, the NSC should disappear with overdoping and the vortices should present a clear zero-bias peak, a challenge for future STS experiments. Alternatively, one can imagine two-channel scenarios. It has been recently proposed that the cuprate superconductivity is lead by oxygen $2p$ rather than copper $3d$ holes Rybicki et al. (2016). We speculate that the copper holes remain localized and form the NSC (in the specific case of Y123, where the zero-bias conductance is large, we cannot rule out a parasitic surface channel as another contribution to the NSC), while the oxygen holes are responsible for the recently uncovered universal Fermi-liquid signatures Barišić et al. (2013) and enter the SC condensate. While the origin of pairing in this condensate remains mysterious, its spectroscopic properties are well described by the BCS theory, as we have demonstrated by unveiling the Caroli–de Gennes–Matricon-like states in the vortex cores.

We thank T. Giamarchi and D. van der Marel for their comments on the manuscript, N. Hussey and G. Deutscher for discussions, and G. Manfrini and A. Guipet for their technical assistance. This research was supported by the Swiss National Science Foundation under Division II. The calculations were performed in the University of Geneva with the clusters Mafalda and Baobab.

I Microscopic model for the superconducting channel of YBa2Cu3O7-δ in the mixed state

I.1 Electronic structure

YBa2Cu3O7-δ (Y123) has two CuO2 planes in the unit cell, CuO chains running along the $b$ axis, and a small orthorhombic distortion with inequivalent $a$ and $b$ axes ($b>a$) Jorgensen et al. (1990). We ignore the CuO chains and the distortion, absent in other cuprates, and therefore irrelevant for high-$T_{c}$ superconductivity. Some effects of the chains on the vortex-core spectra have been studied in Ref. Atkinson, 2009. We furthermore ignore the bilayer splitting for simplicity and represent the CuO2 layers as a one-band tight-binding model on a perfect square lattice with parameter $a=3.85$ Å. We have also performed calculations for a two-band system including bilayer splitting, and found only inessential quantitative differences in the vortex cores. The one-band model is more convenient for large-scale simulations. We use the tight-binding parameters $t_{1}=-281$ meV, $t_{2}=139$ meV, and $t_{3}=-44$ meV determined by photoemission in Ref. Schabel et al., 1998a; *Schabel-1998b-s for the first, second, and third neighbor hopping, respectively, ignoring $t_{4}$ and $t_{5}$ for simplicity. The chemical potential is set to $\mu=-356$ meV for an electron density $n=0.84$, corresponding to optimally hole-doped Y123 with 0.16 hole per unit cell. The dispersion relation measured from the chemical potential is $\xi_{\bm{k}}=2t_{1}[\cos(k_{x}a)+\cos(k_{y}a)]+4t_{2}\cos(k_{x}a)\cos(k_{y}a)+2t_{3}[\cos(2k_{x}a)+\cos(2k_{y}a)]-\mu$. The average group velocity on the Fermi surface is $\langle v_{\mathrm{F}}\rangle=4.11\times 10^{7}$ cm/s. The Fermi surface is shown in Fig. S1(a). Due to a van Hove singularity at $-376$ meV, the DOS has a negative slope in the low-energy region, with more weight for the occupied states [Fig. S1(b)]. We set the amplitude of the $d$-wave order parameter to $\Delta_{0}=19$ meV. In the uniform superconductor, the gap $\Delta_{\bm{k}}=(\Delta_{0}/2)[\cos(k_{x}a)-\cos(k_{y}a)]$ has its maximum at the point $(\pi/a,0.74/a)$ of the Fermi surface, giving coherence peaks at $\pm 17$ meV [Fig. S1(c)]. This amplitude and symmetry of the order parameter follow self-consistently from the Bogoliubov-de Gennes equations for an (instantaneous) attractive interaction $V=-247$ meV on nearest-neighbor bonds.

Figure S1(c) shows the base-temperature zero-field tunneling spectrum of Y123 Bruér et al. (2016) and a possible decomposition in two channels. The superconducting channel (SC) is calculated with the tight-binding model, and the non-superconducting channel (NSC) is the difference. The relative weight of the two contributions is adjusted in such a way that the NSC has no structure—peak or dip—at the edges of the superconducting gap. The resulting NSC has an overall positive slope. The latter is sensitive to the choice of hopping parameters in the SC, and is therefore not a robust feature of the analysis. This slope is irrelevant for our study of vortex cores focusing on energies $\lesssim 50$ meV. In contrast, the dips near $\pm 50$ meV, the peaks near $\pm 30$ meV, and the subgap peaks near $\pm 5$ meV are robust properties of the spectrum measured by STS in regions where superconductivity is suppressed Bruér et al. (2016). The dashed blue line shows the spectrum expected in such a region where the superconducting gap is closed, assuming that the relative weight of the two channels remains unchanged. This is fully consistent with the spectra measured in non-superconducting regions Bruér et al. (2016). A more elaborate modeling was introduced in Ref. Bruér et al., 2016, that involved bilayer splitting as well as an interaction with the spin fluctuations. The main effect of these additional ingredients is to assign (part of) the dips at $\pm 50$ meV to the SC rather than the NSC. As these energies are not our main concern and these sophistications are impractical in view of large-scale vortex calculations, here we disregard them.

I.2 Isolated vortex, self-consistent solution

We solve the Bogoliubov-de Gennes equations self-consistently with a single vortex at the origin using the method described in Ref. Berthod, 2016. The reader is referred to Ref. Berthod, 2016 for all practical details, while here we only give the elements specific to the present calculation. As our Hamiltonian extends up to third neighbors, it spreads the wave function on the lattice with a diamond shape. We therefore consider a finite system with diamond shape and linear size $M$, having $1+2M(1+M)$ lattice sites. We use $M=200$ (80 401 sites) for calculating the self-consistent order parameter and $M=500$ (501 001 sites) for calculating the local density of states (LDOS). The order parameter requires a smaller system because the coherence length imposes a spatial cutoff. With $M=200$, a Chebyshev expansion order $N=6000$, and termination using the Jackson kernel Berthod (2016), the calculation retrieves the exact order parameter within 0.1%. For the LDOS, the spatial cutoff would be set by the mean free path, which is infinite in our model. With $M=500$ and an expansion order $N=2000$, we reach an energy resolution of roughly $3$ meV without perturbations associated with the system’s boundaries.

The self-consistent order parameter is plotted in Fig. S2(a). The quantity $|\Delta(\bm{r})|$ is defined as the sum of the order-parameter modulus on the four bonds surrounding the site $\bm{r}$. It is well fitted by the isotropic ansatz $\Delta(r)=\Delta_{0}/[1+\xi_{0}/r\exp(-r/\xi_{1})]$ with $\xi_{0}=17a$ and $\xi_{1}=29.5a$. The difference between the exact and approximate data is negative along the $x$ and $y$ directions and positive along the diagonals. The self-consistent order parameter indeed displays a small in-plane anisotropy unlike the ansatz, and relaxes faster to its asymptotic value along the diagonals than along the lattice axes, as already found in similar calculations Ichioka et al. (1996); Berthod (2016). A “core size” $\xi_{c}$ may be defined by the condition $\Delta(\xi_{c})=\Delta_{0}/2$, yielding $\xi_{c}=11.5a=4.4$ nm. This agrees very well with the BCS expression of the coherence length $\xi=\hbar v_{\mathrm{F}}/(\pi\Delta_{0})=4.5$ nm if the Fermi-surface average of the velocity is substituted for $v_{\mathrm{F}}$.

The LDOS plotted in Figs. S2(b), S2(c), and S2(d) displays a zero-energy peak that is maximum at the vortex center and decays differently in all directions. The in-plane anisotropy of the LDOS is not a consequence of the in-plane anisotropy of the order parameter, as illustrated in Fig. S2(d), where it is seen that the typical angular dependence of the LDOS remains unchanged if an isotropic order parameter is used. It is also not a consequence of the $d$-wave gap symmetry: repeating the calculation for an $s$-wave gap leads to the same anisotropic LDOS with an un-split peak along the antinodal directions and a split peak along the nodal ones. Note that the LDOS peak in Fig. S2 is a genuine continuum, not the superposition of densely packed discrete core levels as in $s$-wave superconductors Franz and Tešanović (1998); Berthod (2016). In fact, the LDOS anisotropy relates to the dispersion and Fermi-surface anisotropies, which tend to favor low-energy LDOS structures in the directions normal to the Fermi surface. The two traces in Figs. S2(b) and S2(c) span different distances from the core; Figure 1(a) of the main text allows one to compare these traces along the same distance. The specific signature of the core for an isolated vortex [Fig. 1(b) of the main text] is obtained by subtracting the zero-field spectrum (green curve in Fig. S1) from the vortex LDOS.

I.3 Ideal vortex lattices, self-consistent solutions

The self-consistent solution for an ideal square vortex lattice with inter-vortex distance $d=54a$ is displayed in Fig. S3(a). This corresponds to a density of 19 vortices on $90\times 90$ nm2 as observed in the experiment (Fig. 2 of the main text). The vortex lattice is aligned with the microscopic lattice, the nearest-neighbor vortices being found along the $x$ and $y$ directions. The order parameter is maximum in the center of the squares formed by four nearest-neighbor vortices and has saddle points with $|\Delta(\bm{r})|=12$ meV on the lines joining them, leading to a significant spatial modulation. Note that the solution has been constrained to have a maximum gap of 19 meV like in zero field for simplicity—and also because no measurable reduction of the gap size is observed experimentally at this field—requiring a slight increase of the interaction to $V=-260$ meV. $|\Delta(\bm{r})|$ has a rounded shape in the cores, which is captured by the ansatz (generalized for vortex lattices Berthod (2016)) with an increased value $\xi_{0}=115a$ and a reduced value $\xi_{1}=8a$ with respect to the isolated vortex. The core size defined as $\Delta(\xi_{c})=\Delta_{0}/2$ is slightly increased to $\xi_{c}=15.9a=6$ nm compared with the isolated vortex, consistently with previous studies in the quantum regime Berthod (2016).

The corresponding LDOS is shown in Figs. S3(b), S3(c), and S3(d). Although the relative difference between the exact and ansatz solutions is 12% at maximum, the LDOS curves calculated with both order parameters are almost undistinguishable. This is the justification for using the non-self-consistent ansatz when studying the LDOS in disordered vortex configurations, for which a full self-consistent calculation is impractical. The zero-energy LDOS peak is considerably suppressed and broadened in the core with respect to the isolated vortex. We have checked that the vortex-lattice calculation correctly reproduces the isolated-vortex limit as the distance $d$ is increased: both spectra differ by less than 5% for $d\gtrsim 170a$ ($B\lesssim 0.5$ T). At the field considered, however, both the spectral and spatial signatures of the vortex cores differ markedly from those of the isolated vortex seen in Fig. S2. We have also verified that the suppression of the zero-energy LDOS peak is not due to the increased value of $\xi_{c}$ and more rounded order parameter in the core: the spectra shown in Fig. S3 remain qualitatively unchanged if we use the ansatz order parameter with the values of $\xi_{0}$ and $\xi_{1}$ corresponding to the isolated vortex. The reason for a suppressed vortex-core peak can be understood by comparing the low-energy LDOS in Figs. S2(b) and S3(d). Because for the isolated vortex the LDOS extends farther along the (10) direction than along the (11) direction, in the vortex lattice the wavefunctions in different cores strongly overlap and the core states get delocalized. This overlap is suppressed when the vortex lattice is not precisely aligned with the (10) direction and/or the vortex positions are disordered, such that the zero-energy LDOS peak is restored in these situations (see below).

For comparison, we show in Fig. S4 the self-consistent order parameter and LDOS for a square vortex lattice rotated 45*∘* with respect to the tight-binding lattice with an inter-vortex distance $d=38\sqrt{2}a$, which corresponds to the same field as in Fig. S3. There are significant differences between the two vortex-lattice orientations (hereafter I and II), both in the self-consistent order parameter and in the LDOS. While in I the gap has saddle points between nearest-neighbor vortices and maxima between next-nearest-neighbor ones, the situation is reversed in II with the gap maxima between nearest-neighbor vortices. It appears that the order parameter doesn’t move rigidly with the vortex lattice: when rotating the vortex-lattice orientation from I to II, the cores move but the gap maxima and saddle points stay in place. The shape of the core is also quite different in I and II, where a best fit to the ansatz gives $\xi_{0}=9.9a$, $\xi_{1}=15a$, and a core size $\xi_{c}=6.4a=2.5$ nm smaller than in zero field. There is more structure in the case II, because each saddle point is in fact replaced by two saddle points separated by a local minimum. This explains the larger discrepancy between the ansatz and the exact solution, which reaches 32% for II at the local minima. As a result, the differences between the LDOS calculated with the ansatz and self-consistent order parameters are slightly larger in II than in I. These differences remain nevertheless small compared with the qualitative differences between the LDOS in I and II: the zero-energy peak is neither strongly suppressed nor split in II as it is in I; at the position of the local minimum between two saddle points in II, the LDOS has a double peak at zero energy, while at the saddle point in I the LDOS is gapped. We see in Fig. S4(d) that the low-energy LDOS is much more localized in the cores compared with I, which highlights the much weaker wavefunction hybridization along the (10) directions in case II.

The data in Figs. S2, S3, and S4 demonstrate that the vortex-core LDOS is not only a function of field, but also and more importantly a function of the positions of neighboring vortices. Depending on where the neighbors are, the zero-energy LDOS peak may be sharp or not, split or not, etc. Experimentally, one therefore expects variability in the measured vortex-core spectra when the vortex positions are disordered. These figures also show that, in order to investigate theoretically this variability, it is sufficient to work with the ansatz order parameter, whose only inputs are the vortex positions and the values of $\xi_{0}$ and $\xi_{1}$.

II LDOS calculations for disordered vortex lattices

II.1 Disordered vortex configurations

The high sensitivity of the theoretical LDOS to vortex positions prompts us, for a meaningful comparison with the STS experiment, to use in the calculation the vortex positions as they appear under the STM. Three difficulties arise: (i) the LDOS maxima which are accessible experimentally may not sit exactly on the phase singularity points where the order parameter vanishes, due to a polarization of the LDOS by asymmetric supercurrents Berthod (2013, 2016); (ii) we cannot disregard vortices that are outside the STS field of view although we don’t know their positions; and (iii) due to lack of atomic resolution on the Y123 surface, the orientation of the microscopic lattice is only known approximately via the macroscopic twin boundaries. We ignore (i), expected to be a small effect, and locate the vortices inside the STS field of view at the positions of largest $dI/dV$ contrast [see Fig. 4(a) of the main text]. Outside the field of view, we generate vortex positions with the same density as inside, which corresponds to a field of 4.85 T. The positions are chosen at random, however, with a hard-core repulsion constraining the inter-vortex distance to be at least $d_{0}=41a$. This value was selected to be as large as possible: for larger values the random generation process would be stuck, not able to fit in the required number of vortices. The resulting vortex configurations show some degree of order, similar to what is seen in the STS field of view. An example is shown in Fig. S5. The structure factor $S(\bm{k})=\mathscr{N}^{-1}\left|\sum_{n}e^{-i\bm{k}\cdot\bm{R}_{n}}\right|^{2}$, where $\bm{R}_{n}$ are the positions of the $\mathscr{N}$ vortices, is isotropic indicating no orientational order. The angular average $\bar{S}(k)=\mathscr{N}^{-1}\sum_{nm}J_{0}(k|\bm{R}_{n}-\bm{R}_{m}|)$, where $J_{0}$ is the Bessel function, shows oscillations of wavevector $2\pi/d_{0}$ due to the hard-core repulsion. Furthermore, a power-law suppression of $\bar{S}(k)$ is observed for $k<2\pi/d$, where $d=54a$ is the inter-vortex distance in the equivalent ordered square lattice. Such power law is reminiscent of hyperuniformity, i.e., a type of order characterized by the suppression of density fluctuations at long wavelengths Torquato and Stillinger (2003), where $\bar{S}(k)\sim k^{2-\eta}$ with $0<\eta\leqslant 2$ in two dimensions. It is likely that in reality the vortices outside the field of view present more order than the configurations generated by our procedure Maggio-Aprile et al. (1995), but the experimentally available vortex positions are not sufficient for inferring such an order. While certain characteristics of the vortex ordering outside the field of view may influence the LDOS inside, we do not expect this to change any of the conclusions we draw from our analysis.

II.2 Search for a good configuration of vortices

We have generated 600 disordered vortex configurations for various orientations of the microscopic lattice with respect to the STS field of view. For each configuration, we generate the order parameter using the ansatz and the values of $\xi_{0}$ and $\xi_{1}$ corresponding to Fig. S3, and we calculate the LDOS along the four paths displayed in Fig. S6(a), that correspond to paths no6 and no8 in Fig. 4 of the main text. These four traces share a common point [the cross in Fig. S5(a)]: we use the LDOS at this point as the reference spectrum and subtract it from the calculated LDOS along the four paths. The same procedure applied to the experimental data [Fig. S6(b)], allows us to compute a sum of squared differences as the figure of merit for each vortex configuration. We find that the agreement between measurement and simulation is systematically better when the four paths are close to a nodal direction. The best compromise is reached if the STS field of view is rotated by 5.7*∘* relative to the microscopic lattice as shown in Fig. S5(a). This points to a tendency for the nearest-neighbor vortices to align along the crystal axes, and justifies a posteriori our use of the values of $\xi_{0}$ and $\xi_{1}$ corresponding to that orientation rather than that of Fig. S4. We emphasize again that the precise values of $\xi_{0}$ and $\xi_{1}$ play very little role in the LDOS. Figure S5(a) is the best among the 600 configurations; the four traces are compared with the experimental ones in Fig. S6. Using this configuration, we calculate the LDOS in the whole STS field of view with 1 nm resolution. For comparing with Fig. 4(a) of the main text, we compute $dI/dV$ by adding to our theoretical LDOS the non-superconducting channel using the formula quoted in Fig. S1, we evaluate the ratio between the calculated $dI/dV$ at 5 and 17 meV, and thus obtain the map and traces shown in Fig. S6(a) and Fig. 4(b) of the main text.

Figure S7 displays three vortex configurations different from the best one shown in Fig. S5(a) and S6, and the corresponding theoretical spectral traces. These configurations agree reasonably with experiment as well, with a figure of merit within the best 10% out of the 600 considered. One notices, in particular, four spectroscopically very different cores with configuration (a), a reinforcement of the LDOS at intermediate distance in (b), trace $\alpha^{\prime}$, similar to what is seen in Fig. 2(c) of the main text, and split spectra at the center of vortices $\alpha^{\prime}$ and $\delta^{\prime}$ for configuration (c), as seen in the experiment [Fig. S6(b)].

II.3 Mixed-state LDOS and gap symmetry

It is tempting to search signatures of the $d_{x^{2}-y^{2}}$ symmetry of the order parameter in the LDOS around vortices. In the semiclassical regime $k_{\mathrm{F}}\xi\gg 1$, the vortices are slowly varying perturbations of the order parameter compared with atomic distances, and their Fourier components are mostly at low momenta. Therefore, the vortices provide only small momentum transfers and the interaction of the Bogoliubov quasiparticles with the vortex lattice is dominated by forward scattering. In that limit, the details of the Fermi surface are irrelevant and the only source of spatial LDOS anisotropy—apart from the vortex lattice itself—is indeed the order-parameter symmetry. In the quantum regime $k_{\mathrm{F}}\xi\sim 1$ relevant for Y123, however, the vortex lattice scatters Bogoliubov quasiparticles with large momentum transfers of the order of $k_{\mathrm{F}}$ and the LDOS is therefore sensitive to the anisotropy of the Fermi surface. This situation generically leads to the development of LDOS structures along the directions normal to the Fermi surface, a fact well known from studies of impurity scattering Weismann et al. (2009). In Y123, the Fermi surface segments are mainly oriented along the crystallographic directions [see Fig. S1(a)], such that one expects structure in the vortex-lattice LDOS along the (10) and (01) lattice directions, irrespective of the order-parameter symmetry. This is confirmed by our numerical results shown in Fig. S8: for both $d$- and $s$-wave symmetries, the LDOS structures are aligned with the microscopic lattice at all energies. No structure is observed along the (11) and equivalent directions, which are the directions of the gap nodes in reciprocal space. Thus the expectation that the LDOS would “leak” out of the vortices along the directions of the gap nodes is not confirmed. There are nevertheless differences between the LDOS calculated for $d_{x^{2}-y^{2}}$ and $s$ symmetries. The vortex states are more localized in the $s$-wave case, leading to a better contrast, especially at low energy. However, according to these simulations, an unambiguous determination of the order-parameter symmetry based on experimental LDOS maps around vortices appears to be hopeless.