Fragile superconductivity in a Dirac metal

In topological semimetals, the electronic band structure features relativistic linearly-dispersive band crossings. This degeneracy gives rise to quasiparticles in condensed matter systems known as Dirac or Weyl fermions, analogous to those found in quantum field theory. Realizing these in bulk materials is an exciting prospect because they bring about topological protection, physical properties beyond those seen in the semiclassical regime, and novel phases of matter. For example, Dirac semimetals are often found to have large linear transverse magnetoresistance, and negative longitudinal magnetoresistance due to the chiral anomaly ^{1, 2}. Furthermore, new phases of matter can be realized when these topological fermions influence other electronic or magnetic properties, and vice-versa. Exotic topological phases like a monopole superconductor are believed to arise in superconducting magnetic Weyl semimetals with inversion symmetry, whose gap functions are given by monopole harmonics ³. On the way to realizing such monopole superconductors, it is natural to study centrosymmetric, non-magnetic superconductors that harbor Dirac fermions to understand the physical phenomena that emerge. One family of interest is the set of square-net materials that contain pnictogens like Sb or Bi, where Dirac fermions arise due to nonsymmorphic symmetry ^{4, 5}. Chamorro et al.⁶ found a large linear magnetoresistance and small effective masses in the square-net material \chLaCuSb2 that were attributed to Dirac fermions. Earlier reports of superconductivity ^{7, 8} in this material makes it an ideal system to study the interplay between Dirac fermions and superconductivity, however superconductivity was not detected by Chamorro et al.

Here we resolve this contradiction, demonstrating that there is an extreme sensitivity of the physical properties to copper stoichiometry in \chLaCuSb2. We have grown single crystals of \chLaCuSb2 with varying copper content to study the effect of disorder on the superconducting state. In optimized samples with the highest superconducting $T_{c}$ and closest to ideal stoichiometry, we used specific heat, susceptibility, magnetization, and muon spin rotation and relaxation ($\mu$SR) to characterize the superconductivity and its critical behavior. We utilize a combination of simple free-electron models, Bardeen-Cooper-Schrieffer (BCS) models, and tight-binding models to understand the fragile and anisotropic superconductivity in \chLaCuSb2.

Despite the anisotropy we observe experimentally, we gain quantitative insight into the superconducting state by calculating thermodynamic quantities in a one-band, isotropic free-electron model as a first approximation. The relevant thermodynamic and superconducting quantities are tabulated in Table 1. In particular, we find a coherence length of $\xi\approx 1\mu$m, whereas the mean-free path is $\ell\approx 0.06\mu$m, suggesting \chLaCuSb2 is in the dirty limit²⁰ since $\xi\gg\ell$. The relevant expression ²¹ for calculating the Ginzburg-Landau (GL) parameter in the dirty limit is $\kappa=(7.49\times 10^{3})\gamma_{nV}^{1/2}\rho_{0}$ (with all quantities in cgs units), where $\rho_{0}$ is the value of the low-temperature residual resistivity plateau. The out-of-plane GL parameter $\kappa_{c}$ depends on currents generated in the $ab$ plane, so using $\rho=\rho_{0a}$ for in-plane resistivity, we find $\kappa_{c}=0.398(7)$. That $\kappa_{c}<1/\sqrt{2}$ is consistent with Type-I superconductivity for $H\parallel c$. Extrapolating our magnetization data to $T=0$, we have $H_{c}=59.8(1.0)$ Oe and $H_{c2}=172(6)$ Oe and thereby estimate²⁰ $\kappa_{a}=H_{c2}/\sqrt{2}H_{c}=2.03(8)$. Next, we use conventional BCS theory to estimate the thermodynamic critical field $H_{c}$. This field is related to the gap function and density of states ²², and in cgs units is given by $H_{c}(0)=1.764\sqrt{6/\pi}\cdot\gamma_{nV}^{1/2}T_{c}$. From this we estimate that $H_{c}(0)\approx 68(1)$ Oe, in remarkable agreement with the measured phase diagram (see Extended Fig. 14). We attribute the small critical field to the low $T_{c}$ and the small density of states in this Dirac material.

Although Type-I superconductivity is mainly found in elemental superconductors and several binary compounds (e.g. noncentrosymmetric \chBeAu ¹⁷ and Dirac semimetal candidate \chPdTe2 ²³), it has also been observed in ternary compounds including \chLaRhSi3 ²⁴ and \chLiPd2Ge ²⁵, among others. However, anisotropic Type-I and Type-II superconductivity is not as widely reported. In the theory of conventional anisotropic superconductors ^{12, 26}, the anisotropy factor $\gamma=\kappa_{a}/\kappa_{c}$ is related to the effective masses as $\gamma^{2}=m_{c}/m_{a}$. The anisotropy of the GL parameter is directly related to the anisotropy of the Fermi surface, as schematized in Fig. 4e. Furthermore, it is at the same time possible to satisfy that $\kappa_{c}<1/\sqrt{2}$ and $\kappa_{a}>1/\sqrt{2}$ when $\gamma\gg 1$, implying a superconductor whose type depends on the applied magnetic field direction ¹². In particular, the direction of the field relative to the crystallographic axes will either favor or disfavor vortices in the superconducting order parameter, based on the free energy in distinct crystallographic directions. While rare, this behavior has been observed in \chC8K ²⁷ and TaN ^{21, 28}, where the angular dependence of the critical field was explicitly studied. However, such an anisotropy seems not to be well-studied in ternary systems. Considering the previous estimates of $\kappa_{a}$ and $\kappa_{c}$, the large anisotropy parameter $\gamma^{2}\approx 26(2)$ is a consequence of the small in-plane effective masses of this Dirac material.

From specific heat, we have pointed to the possibility of multigap behavior or an anisotropic gap, supported by the findings of our tight-binding model. However, determining the temperature-dependence of the gaps using our Hamiltonian requires fitting many free parameters. In Extended Fig. 12, we highlight several toy models that capture the overall features in specific heat and superfluid density, but do not capture all features. One possibility to measure this anisotropic gap is using quasi-particle interference with Scanning Tunneling Microscopy (STM). The quasiparticle interference signal should be significantly enhanced near the gap edge, which may be directly compared to the anomaly in the specific heat. We find it quite compelling to obtain a more complete microscopic understanding of the superconducting state of \chLaCuSb2 using a combination of theory and probes like STM. Overall, despite the fragility of the superconducting state in \chLaCuSb2 to stoichiometry, pressure, and magnetic fields, it is worth exploring related compounds, perhaps ones that contain magnetic moments, to continue the search for exotic topological phases like the monopole superconductor.

Sample synthesis and data analysis on the superconducting state was performed by C.J.L. under the supervision of C.L.B. and T.M.M. The $\mu$SR experiments were carried out by C.J.L., J.G., M.P., S.S., and P.O. under the supervision of C.L.B and G.M.L. The tight-binding analysis was performed by J.Z. under the supervision of Y.L. The DFT calculations for the band structure were performed by S.G. The heat capacity measurements were performed by C.J.L. and T.H. The He3 susceptibility and magnetization measurements were performed by C.J.L., T.B. and J.R.C. Quantum oscillations were measured and analyzed by G.G. under the supervision of B.J.R. Single-crystal X-ray diffraction measurements and refinements were performed by M.A.S. TEM measurements were performed by K.L. The hydrostatic pressure measurements were performed by K.M. and D.K.B. under the supervision of Y.U. and S.N. All authors contributed with comments and edits.

Synthesis: Single crystals of \chLaCuSb2 were grown by the self-flux method. Cut pieces from lanthanum ingot (Ames Laboratory, 99.99%), pieces of copper (Alfa Aesar, 99.999%), and antimony shot (Strem Chemicals, 99.9999%) were weighed with a total mass of about 4 g in various molar ratios and placed in an alumina crucible. An inverted catch crucible with an alumina strainer was placed atop the first crucible, and both were sealed in a quartz tube under partial pressure of argon gas. The ampules were heated to 1070^∘ C at a rate of 100^∘ C/h and held for 12 h, then cooled to 670^∘ C at a rate of 4^∘ C/h before centrifuging. Inspired by previous literature on isostructural \chLaAgSb2, which was reported to produce stoichiometric samples with starting ratio $0.045:0.091:0.864$ (or roughly, $1:2:19.2$) ^{29, 30}, we used varying compositions $1:N_{\text{Cu}}:19.2$ with $1\leq N_{\text{Cu}}\leq 6$ as a parameter to tune the chemical potential of the system and grow crystals with different Cu content.

Polycrystalline samples were also prepared using a simple reaction of the elements in a quartz tube. Powders with ratios $1:\delta:2$ with $\delta=0.8,1.0,1.2$ were synthesized by first melting about 1.0 g of the elements in a quartz tube, using a step furnace at 600^∘ C for 24 hours. This polycrystal was then ground and reheated to $800^{\circ}$C for roughly 12 hours. Further heating of the powder at high temperatures, or low temperatures for extended periods of time, was disadvantageous due to the decomposition of the structure, as evidenced by larger relative percentages of secondary phases \chSb and \chCu2Sb. We compared the lattice constants (in particular, the $c/a$ ratio) of the polycrystalline samples to those of our single crystals, and to those from the literature.

X-ray diffraction: Powder XRD data were collected at room temperature using a laboratory Bruker D8 Venture Focus diffractometer with LynxEye detector in the range from 10-80 degrees. Refinements on the powder XRD data were performed using Topas 5.0 (Bruker). The occupancy of La and Sb ions were constrained to 100%, while the Cu occupancy was refined freely. We also refined the strain parameter and preferred orientation, due to the layered nature of the crystals.

Single-crystal XRD data were acquired at 110(2) K using a SuperNova diffractometer (equipped with Atlas detector) with Mo K$\alpha$ radiation ($\lambda=0.71073$ Å) under the program CrysAlisPro (Version CrysAlisPro 1.171.39.29c, Rigaku OD, 2017). The same program was used to refine the cell dimensions and for data reduction. The structure was solved with the program SHELXS-2018/3 and was refined on F2 with SHELXL-2018/3 ³¹. Analytical numeric absorption correction using a multifaceted crystal model was applied using CrysAlisPro. The temperature of the data collection was controlled using the system Cryojet (manufactured by Oxford Instruments). For all samples in the ZrCuSiAs structure type, the occupancy factor for Cu1 was refined freely. For samples in the defect-\chCaBe2Ge2 structure type, the additional Cu site (denoted Cu1’) was necessary to obtain good fits to the data. The occupancy factor for both Cu1 and Cu1’ were refined freely, and the reported occupancy is $x_{\mathrm{Cu}}=x\mathrm{[Cu1]}+x\mathrm{[Cu1^{\prime}]}$. Crystallographic data tables are included in the Supplementary Information.

Specific heat: Heat capacity measurements were performed in a Quantum Design Physical Properties Measurement System (PPMS) with the dilution refrigerator option. The magnetic field was degaussed above 3.8 K to minimize effects of persistent fields. Measurements were performed on a single crystalline sample of LaCuSb₂ of 10.7(1) mg mass oriented such that the applied magnetic field was along the nominal $c$-axis, with reported fields depicted in Fig. 4 between 0 Oe and 75 Oe. All measurements were performed in a fixed magnetic field and measured upon cooling, with a minimum temperature between $T=0.05-0.1$ K and a maximum temperature of $T=3.8$ K.

Magnetization: Magnetization measurements were performed in a Quantum Design (QD) Magnetic Properties Measurement System (MPMS) with QD iHelium3 He3-insert. To extract the superconducting dome, we measured samples with nominal applied fields of 2 Oe along the $a$-axis, after degaussing at temperatures above $T_{c}$. The samples were cut with a large aspect ratio along the $a$-axis such that the demagnetization corrections were small. Measurements on the optimal sample were performed taking into account the non-zero demagnetization factors. We performed isothermal measurements in applied field by first degaussing at high temperature above $T_{c}$. We then cooled in zero field to the appropriate temperature, and measured from zero to fields well past the point that magnetization vanished to study $4\pi M$ versus $H_{\mathrm{int}}$.

Electrical transport: Shubnikov-de Haas oscillations were measured in a single crystal cut from the same crystal used for other thermodynamic measurements, at temperatures down to 0.3 K and fields up to 16 T. Resistivity and Hall effect measurements above 2K were performed on a PPMS using the AC Transport (ACT) option. For all resistivity measurements we prepared a polished bar-shaped crystal with the current applied along the $a$-axis, using a four-probe configuration consisting of platinum wires and Epo-Tek silver epoxy. For Hall effect measurements we used the five-probe method which allowed the Hall signal to be tuned to zero in zero applied field. Low-temperature resistivity measurements were performed in a PPMS using a Lake Shore Model 372 AC Resistance Bridge. We applied an AC current of 316 $\mu$A at a frequency of 13.7 Hz. The field was applied along the $c$-axis after degaussing the magnet at 1.1 K (above the transition) to reduce the effect of persistent magnetic fields. Measurements were taken upon cooling in various applied fields.

Hydrostatic pressure: Measurements under hydrostatic pressure were performed using a Bluefors dilution refrigerator down to 0.05 K and up to 2.0 GPa at the Institute for Solid State Physics (ISSP), The University of Tokyo. To track the variation of the superconducting transition temperature under hydrostatic pressure, the AC magnetic susceptibility of a sample was measured with a mutual induction method at a fixed frequency of 317 Hz with a modulation field of about 1 Oe. Measurements were performed on the optimized samples $N_{\mathrm{Cu}}=2$, cut from the same crystal used for magnetization, specific heat, and resistivity measurements. For applying pressure, a piston-cylinder cell made from nonmagnetic BeCu and NiCrAl alloys was used with Daphne 7373 as the pressure transmitting medium. The pressure was determined from the superconducting transition temperature of Pb.

Muon spin rotation: Zero-field and transverse-field muon spin rotation and relaxation measurements were performed at the TRIUMF facility in Vancouver, Canada. A spectrometer incorporating a dilution refrigerator was used on the M15 beamline, to allow for measurements down to 20 mK. The setup makes use of a superconducting magnet to allow for magnetic fields up to 5 T, and resistive coils for finer control and field-zeroing. The magnetic field was applied horizontally, parallel to the direction of the muon beam. In the ZF measurements, the muon spin was (anti)parallel to the beam direction, while in TF measurements the muon spin was perpendicular to the field and beam direction. Single crystals grown with the optimal ratio $N_{\mathrm{Cu}}=2$ and with the greatest thickness (average 1.0 mm) were cut along the $a$-axis. We used multiple co-aligned single crystals, totaling a mass of about 0.98g, to maximize the measured signal and reduce background from muons not implanted in \chLaCuSb2. The crystals were first mounted with the $a$-axis parallel to the applied field for the first measurement, and the same crystals were then individually rotated and remounted with the $c$-axis parallel to the applied field for the second measurement. The samples were mounted on a silver cold finger using a mixture of Apezion N grease and copper-loaded Cry-Con grease to ensure good thermal contact. For the field along the $a$-axis, we used $H_{\mathrm{app}}=40$ Oe. For the field along the $c$-axis, we used $H_{\mathrm{app}}=40$ Oe to access low critical fields and $H_{\mathrm{app}}=10$ Oe to access higher critical fields. All data at constant field were simultaneously refined and fit for various temperatures using the musrfit program. ³²

The flux-grown crystals were plate-like with average surface areas of 1 cm², often limited by the crucible diameter, and with varying thickness. The least Cu-rich fluxes produced crystals with the most appreciable thickness, compared to the most Cu-rich fluxes having large surfaces areas and decreased thickness. Grown crystals with copper-rich fluxes often came out with layers of solidified \chSb-\chCu2Sb flux on the surface, which could be removed with mechanical polishing. The presence of these phases is expected from the ternary phase diagram reported previously \citesecGschneidner. Inclusions, primarily Sb and less-so \chCu2Sb, were present in all samples, and Sb was used as a reference for the lattice constants for powder XRD. Importantly, \chSb is not superconducting at ambient pressure \citesecWittig and \chCu2Sb is only superconducting below 0.085 K \citesecAndres, and thus cannot account for the superconducting signals observed near 1 K in our samples. However, the nature of off-stoichiometry may effect the physical properties we observe. For example, removal of entire planes of Cu would in principle shrink the unit cell to an effective \chLaSb2 structure, which (in the orthorhombic structure) is known to be superconducting \citesecGuo, Ruszala

The fact that $N_{\mathrm{Cu}}$ can tune a wide range of compositions is likely due to the change of thermodynamic chemical potential and stable phases that can emerge in the \chLaCu_xSb2 solid-solution-type structure. If the main role of the Cu off-stoichiometry is interpreted as the introduction of non-magnetic impurities (as opposed to stacking defects, see TEM results below), in the form of vacancies or interstitial sites, then according to Anderson’s theorem on dirty s-wave superconductors, such defects should have no effect on $T_{c}$ \citesecAnderson. There may be several ways to reconcile this. One possibility is that the small effects of off-stoichiometry may be affecting the Fermi level and density of states, as suggested in resistivity measurements (see 4), thereby resulting in stronger dependence of $T_{c}$ as expected under conventional BCS theory ²⁰. Another possible explanation is that the gap function itself is anisotropic leading to strong dependence of $T_{c}$ on the scattering rate, a well-known result from Gor’kov \citesecPetrovic, Gorkov. Indeed, an anisotropic gap function is expected from our tight binding model for the pairing gap function. It is possible that one or both effects could explain the drastic suppression of $T_{c}$ found in the superconducting dome.

A Cu-deficient sample with $N_{\mathrm{Cu}}=1$, and a near-stoichiometric sample with $N_{\mathrm{Cu}}=2$, were ground into a fine powder for use in TEM. The samples used were from the same crystals cut for thermodynamic measurements. Microcrystallites nearly oriented along $a$ and $c$ were studied for their fringe patterns to estimate the lattice constants, and to gauge small variations in fringe width, which may give evidence of stacking defects in the unit cell. This high-resolution TEM data precludes variations in the lattice constant within $\pm 0.44$Å  in 27 unit cells. We used the GMS-3 software to integrate over regions shown in Extended Fig.  ‣ 2, and extract the intensity as a function of position over the length of the rectangles. Near each maxima, we fit the location of the central peak using a simple quadratic fit. This was used to deduce the distribution of distances between successive fringes. The resulting histogram shows spread within one pixel, i.e. a spread of 0.44 Å. The lack of clear bimodal distributions outside the resolution of the scan indicates there are no obvious stacking disorders that may produce \chLaSb2. Given the off-stoichiometry on the order of 0-5%, this suggests possibly statistical vacancies of Cu atoms as suggested in \citesecKlemenz. Interestingly, the samples with smaller copper flux ratios $N_{\mathrm{Cu}}\approx 1$ that were cut into the bulk of the sample appeared to have a copper-colored tint on the exposed cut surface after a period of several months. This indicates the copper ions may be mobile in the structure at room temperature, made possible by the presence of vacancies and a path by which ions can travel.

The volume magnetization was deduced from the total moment $\mu_{\mathrm{tot}}$ from SQUID measurements by taking into account the mass $m$ of the samples, and using theoretical density $\rho_{\mathrm{th}}=7.48$ g cm^-3, to compute $M=\mu_{\mathrm{tot}}/V=\mu_{\mathrm{tot}}\rho_{\mathrm{th}}/m$. In magnetization data, we computed the internal field to correct for the effects of the demagnetization factor,

where $H_{a,i}$ is the applied field in the direction $i$, $M_{i}$ is the (volume) magnetization, and $N_{i}$ is the demagnetization factor with $0\leq N_{i}\leq 1$. For diamagnetic samples in a rectangular prism geometry, the demagnetization factor could be estimated with \citesecProzorov

where $A,B$ ($C$) are the lengths of the sides perpendicular (parallel) to the applied field. Anisotropic magnetization measurements on an optimized sample $(N_{\mathrm{Cu}}=2$) were performed on a nearly-rectangular-prism sample with approximate dimensions 3.2 $\times$ 1.52 $\times$ 0.84 mm³ and mass $27.54$ mg. From this, we estimate a demagnetization factor $N_{a}\approx 0.18$ for fields along the $a$-axis and $N_{c}\approx 0.62$ for fields along the $c$-axis. In the $\mu$SR measurement, our collection of co-aligned samples formed the shape of a rectangular prism with effective dimensions 15 $\times$ 12 $\times$ 1.4 mm³. Here the magnetic field was always applied parallel to the thinnest dimension, yielding a demagnetization correction factor of about $N\approx 0.86$.

Quantities like the superconducting transition in susceptibility were estimated by fitting to simple models of a diamagnetic response assuming a Gaussian distribution. Supposing the susceptibility is exactly $4\pi\chi=-1$ below $T_{c}$ and exactly $4\pi\chi=0$ above $T_{c}$, we get a step-like transition at $T_{c}$. For $T_{c}^{\prime}$ distributed in a Gaussian distribution about some mean $T_{c}$, this becomes

For the ambient-pressure measurements, the susceptibility for all samples measured showed sharp transitions and no broad tails at low temperatures, with all samples having a saturated susceptibility by 0.4 K. In all cases $\sigma$, the width of the superconducting transition, was used to estimate the error bar on $T_{c}$, rather than the error reported by the fit routine. The susceptibility for various values of $N_{\mathrm{Cu}}$ is shown in Extended Fig. 2a.

For the high-pressure measurements above 1.7 GPa, the AC data were fit using a modified form due to the double-peak nature of the transition. We model this double-peak feature coming from the successive transition of two gaps opening at $T_{c}$, for the sake of extracting $T_{c}(p)$ and $T^{*}$:

where $A$ is the voltage amplitude, $f$ is the fraction of each component, $\sigma_{1}$ and $\sigma_{2}$ are the widths of the transitions, and $C$ is a constant voltage offset.

Ambient-pressure susceptibility data are shown in Extended Fig. 2b, highlighting the zero-field cooled (ZFCW) and field-cooled (FCW) curves for the optimized sample (measured upon warming), using $4\pi\chi_{0}=4\pi\chi_{v}(1-N_{i})$ where $\chi_{v}$ is the molar susceptibility calculated with the nominal applied field value. Note that $4\pi\chi$ is greater than $-1$, due to the large relative error suggesting the applied field might have been slightly larger than the reported 2 Oe. While the FC volume susceptibility $4\pi\chi_{v}$ is an indication of the Meissner fraction, errors in field calibration or density estimations resulted in inferred Meissner fractions greater than $100\%$. However, the demagnetization-corrected magnetization versus applied field demonstrate the $4\pi\chi_{v}=-1$ relation in Fig. 4c,d.

To determine the critical field from $c$-axis magnetization data $H_{\mathrm{int,c}}(T)$, we find the field that results in a discontinuous transition into the normal state ($M\neq 0$ to $M=0$). This was done by fitting the $H_{\mathrm{int}}$ data versus $4\pi M_{c}$ to a constant value, for each fixed temperature. For magnetization along the $a$-axis, we define the critical field $H_{c1}$ to be the inflection point of the magnetization. Numerically, this was determined by taking the derivative data $f=\differential(4\pi M)/\differential H$, and fitting it near $H_{c1}$ to the piecewise function

The critical field $H_{c1}$ was defined as the field where $4\pi M$ changes slope, or $f(H_{c1},H_{0})=0$, thus

Furthermore, $H_{c2}$ was determined by fitting $4\pi M$ as a linear function at high fields (for small values of $4\pi M$),

The critical field was defined such that $4\pi M(H_{c2})=0$, or $H_{c2}=-c/b$ from the fitted parameters. In all of these fits, the uncertainty was obtained through error propagation.

We can finally extract the zero-temperature critical field by fitting the data to conventional models for the temperature dependence:

where $t=T/T_{c}$ is the reduced temperature. This was done for the data in magnetization, specific heat, and $\mu$SR where we assumed only the $a$-axis oriented samples have a critical field $H_{c2}$, while the remainder of the data were fit for $H_{c1}$.

The reason why superconductivity exists in some samples and not in others with differing copper content is an important question. Given the complicated band structure with various bands crossing the Fermi level, it is not out of the question that tuning the stoichiometry will affect disorder scattering. Indeed, the resistivity of the optimized sample is about an order of magnitude lower than that reported for \chLaCuSb2 in the literature ^{6, 8}. As seen in Extended Fig.4a, the residual resistivity at low temperatures is $\rho_{0a}=1.883(15)\;\mu\Omega$-cm with a residual resistivity ratio (RRR) of about 14, whereas previous samples have been in the $~{}\sim 1$ m$\Omega$-cm range. Furthermore, the low-temperature residual resistivity is the lowest in the optimized sample $N_{\mathrm{Cu}}=2$, whereas it is larger for other samples near the endpoints of the superconducting dome. Importantly, the large linear magnetoresistance seen in Extended Fig. 4b appears resilient in \chLaCuSb2, as it is found in the normal state of the superconducting samples as well as in the samples from Chamorro et al. ⁶. This suggests that the Dirac electrons still play an important role at low temperatures in our samples.

Besides effects of disorder, varying the copper concentration can also affect the chemical potential, and thus the carrier density. To study this, Hall effect measurements were taken using a five-probe configuration with AC Transport in the PPMS, on a sample cut from the same crystal as used in other measurements. The lead separation was $0.38$ mm and the sample thickness was $0.15$ mm. The applied current was 40 mA with frequency 103 Hz. Symmetrization difficulties resulted in apparent non-linearity. Regardless, the Hall coefficient gave a carrier density of approximately $4.16(2)\times 10^{22}$ cm^-1. Interestingly, the sign of the charge carriers in the sample with $N_{\mathrm{Cu}}=2$ was positive, while the sign for samples from Chamorro et al. ⁶ was negative. This suggests the change in stoichiometry also results in a change in the Fermi level, which can result in different contributions from electron or hole charge carriers. In our work, there is also an apparent trend that the end members of the superconducting dome have a smaller carrier density than the optimized sample, for example seen in Extended Fig. 4c. Off-stoichiometry not only affects the sample quality but also the electronic structure, and may correlate with the presence of superconductivity at low temperatures. However, the carrier densities across the superconducting dome are within an order of magnitude of each other, so the density of states is likely to change in proportion to the small changes in the Copper content $x_{\mathrm{Cu}}$.

We also studied the resistivity in the superconducting state for the optimized $N_{\mathrm{Cu}}=2$ sample. The DR-temperature resistivity data in applied fields are shown in Fig. 4d. Like the susceptibility, the temperature dependence of the resistivity was fit assuming Gaussian broadening:

with $\rho_{0a}$ the low-temperature residual resistivity. The resistivity was measured for applied magnetic field along the $c$-axis and currents along the $a$-axis, varying temperature at fixed fields. In these measurements we did not reach a zero-resistance state, but the temperature of the midpoint $T_{c}$ found in the resistivity drop is consistent with $T_{c}$ from other thermodynamic measurements, in low fields $H<25$ Oe. At fields higher than $H_{c1}$, the width $\sigma$ became exceedingly large and the resistance did not saturate at low temperatures, so it was difficult to estimate $T_{c}$ and $\sigma$. This also suggests that the surviving superconducting state at these higher fields is due to percolation or filamentary superconductivity, and is not necessarily a bulk response. For this reason, we did not fit the resistivity data to extract critical field(s) in Extended Fig.14.

Overall, resistivity data for various samples highlights the importance of the copper stoichiometry. While the susceptibility data indicate changes in the superconductivity, we find concomitant changes in the resistivity and Hall effect that suggest these phenomena are related. However, despite changes in defect density, Fermi level, and carrier densities, we still find large linear magnetoresistance in all samples suggestive of the presence of Dirac fermions.

To extract the transition temperatures from specific heat, we first fit the data $c_{p}/T$ to a spline near the transition temperature. The derivative of this spliced data was taken, and the resulting derivative data were fit to a Gaussian function. The transition temperature was then reported as the midpoint and standard deviation of this fit function, at various fields, for use in the phase diagram.

In the normal state above 1 K, we extrapolate the Sommerfeld coefficient and phonon contribution $c_{p}=\gamma_{n}T+\beta_{3}T^{3}$ to zero temperature and to lowest order. This yields a Sommerfeld coefficient $\gamma_{n}=4.78(1)$ mJ/mol-K² and a phonon contribution $\beta_{3}=0.571(2)$ mJ/mol-K⁴. With $d=4$ atoms per formula unit \chLaCuSb2, the corresponding Debye temperature $\Theta_{D}=(12\pi^{4}Rd/5\beta_{3})^{1/3}=238.8(3)$ K. This is consistent with previously reported values ⁸.

In the isotropic free-electron model, the molar specific heat Sommerfeld coefficient $\gamma_{n}$ in the normal state can be used to determine the specific heat effective mass of the charge carriers \citesecAshcroft,

where $\gamma_{nV}$ is the volume specific heat Sommerfeld constant, and $V_{0}$ and $V_{\mathrm{fu}}=\frac{1}{2}V_{0}$ are the volumes of the unit cell and one formula unit, respectively. Using the measured Sommerfeld constant, along with the carrier density derived from the Hall effect, we find a specific heat effective mass of $m^{*}=1.44(1)m_{e}$. This is to be contrasted with the low effective in-plane masses of charge carriers in our previous report ⁶, likely due to the anisotropy of the Fermi surface.

From specific heat one can obtain the electron-phonon coupling parameter $\lambda_{e-p}$, given by

where $\mu^{*}\approx 0.13$ for intermetallic superconductors. We find a value of $\lambda_{e-p}\approx 0.466$, which is in contrast to the value assumed in a previous theoretical work finding $\lambda_{e-p}\approx 0$ \citesecRuszala. However, the free-electron value $\gamma_{b}\approx 2.85$ mJ/mol-K² and the experimental value differ by the enhancement factor $\lambda_{e-p}=\gamma_{n}/\gamma_{b}-1\approx 0.678$, in reasonable agreement with the calculated value in Eq. (14).

The main contribution to the total specific heat at lowest temperatures is from the nuclear Schottky anomaly. This is ascribed to the interaction of nuclear quadrupolar moments with local electric field gradients at the nucleus (see SI 5 for more details on the corresponding modeling). This contribution is dependent on both the spinful nuclear moments, the point group symmetry of the ions involved, and is to be expected for $T<100$ mK in low symmetry solids containing Cu \citesecCaspary and/or Sb \citesecOrtiz. The anomaly also grows with magnetic field, as seen in other Sb-based superconductors \citesecAoki. To extract the electronic specific heat, we model the total contribution to the specific heat for $T<0.35$ K by the equation

where $c(H)$ and $\Delta(H)$ are the phenomenological parameters used to model the activated behavior associated with the knee-like feature; $A(H)$ is related to the quadrupole coupling and increases in applied magnetic fields as the nuclear spin states undergo additional Zeeman splitting. In this way, we simultaneously fit the different contributions to separate out the Schottky anomaly (dominant at low temperature) and phonons (dominant at high temperature) from the electronic contributions.

The entropy of the electronic charge carriers is calculated as

For free electrons, $\Delta S_{\mathrm{el}}(T)/\gamma_{n}T=1$ at all temperatures. In a superconducting sample, $\Delta S_{\mathrm{el}}(T)/\gamma_{n}T=1$ for $T>T_{c}$ in the normal state, and there is an entropy balance at $T_{c}$ such that $\Delta S_{\mathrm{el}}(T_{c})/\gamma_{n}T_{c}=1$. For $T<T_{c}$ in the superconducting state, $\Delta S_{\mathrm{el}}(T_{c})/\gamma_{n}T_{c}<1$ and decreases to zero as $T\to 0$. Extended Fig. 6 shows the calculated entropy, which is close to 1 at $T_{c}$ and above. This indicates the sample is a bulk superconductor and that the sharp drop in the electronic specific heat at $T^{*}$ is intrinsic. Note that the sample coupling reported by the PPMS decreases at decreasing temperatures, due to the combination of the exponentially-activated thermal conductivity and the large sample mass required for a good measurement signal. However, this sudden drop near $T^{*}$ is necessary to obey entropy balance at $T_{c}$, where $\Delta S(T_{c})=\gamma_{n}T_{c}$. The slight overshoot of the entropy $\Delta S_{\mathrm{el}}/\gamma_{n}T_{c}>1$ is likely a result of the uncertainty associated with subtracting the large nuclear Schottky anomaly when isolating the electronic contribution to the specific heat capacity. Furthermore, around phase transitions the specific heat can be difficult to extract $C_{p}$ precisely, adding further error to calculating the entropy. These limit the rigor possible when fitting to multi-gap models, as these BCS and self-consistent models strictly obey entropy balance.

The existence of Dirac fermions in \chLaCuSb2, coupled with any unconventional time-reversal symmetry (TRS) breaking, would make \chLaCuSb2 a prime material candidate in the search for monopole superconductivity. TRS breaking is not a priori expected in \chLaCuSb2, since the nominal crystal structure P4/nmm is centrosymmetric, and all ions in \chLaCuSb2 are nonmagnetic. To determine whether TRS is broken below $T_{c}$ in \chLaCuSb2, we used $\mu$SR in zero field as a way to search for inhomogeneous magnetic fields. Extended Fig. 8 shows the asymmetry $A(t)$ as a function of time for temperatures above and below $T_{c}$. At all temperatures, we fit the asymmetry to the equation

where $G_{KT}(t)$ is the Kubo-Toyabe function due to random fields from nuclear moments; $\Lambda$ is the temperature-dependent relaxation rate; and $F$ is the fraction of muons that stop in \chLaCuSb2 as opposed to the silver mounting plate. In exotic superconductors that break TRS, $\Lambda(T)$ increases with decreasing temperature due to larger field inhomogeneity that often appears below $T_{c}$, such as from domains related by TRS. Simultaneous fits of our high- and low-temperature ZF spectrum reveal a relaxation rate $\Lambda(T)$ that is essentially constant with temperature (within error), where $\Lambda=9(4)\times 10^{-3}$ $\mu$s^-1 at 1.293(3) K and $\Lambda=11(4)\times 10^{-3}$ $\mu$s^-1 at 0.017(1) K. That is, the maximum field size that could exist consistent with these data is $\Delta\Lambda/\gamma_{\mu}=0.02(7)$ Oe.

For a Type-II superconductor, the relaxation rate $\sigma$ is related to the superconducting relaxation rate by $\sigma=\sqrt{\sigma_{\mathrm{SC}}^{2}+\sigma_{\mathrm{n}}^{2}}$, where $\sigma_{\textrm{n}}$ is the temperature independent nuclear spin induced relaxation rate. We estimate the London penetration depth from $\sigma_{\mathrm{SC}}$ as follows. Firstly, for $H\parallel a$, muons probe the London penetration depth $\lambda_{ac}=\sqrt{\lambda_{a}\lambda_{c}}$ due to superconducting currents flowing in the $ac$ plane \citesecLiarte. Assuming that the field ratio $b\equiv H_{\mathrm{app}}/H_{c2}$ satisfies $0.13/\kappa^{2}\ll b\ll 1$, where $\kappa$ is the GL parameter, the London penetration depth and superconducting relaxation rate are related by \citesecBrandt

where $\Phi_{0}\approx 2.067\times 10^{-15}$ Wb is the magnetic flux quantum. With an applied field of $H_{\mathrm{app}}=40$ Oe, we were above $H_{c1}=32(1)$ Oe to set us firmly in the mixed phase, and such that $b\approx 40$ Oe$/172$ Oe $\approx 0.23$ satisfies the condition that $0.13/\kappa^{2}\ll b\ll 1$, assuming $\kappa_{a}\sim 1$.

In this way, we were able to estimate the geometric mean of the zero-temperature anisotropic London penetration depth, $\lambda_{ac}(0)\approx 408(2)$ nm. We note that more accurate estimations of $\lambda_{ac}(0)$ will involve functions of $b$ and $\kappa$, which we do not extract directly from the $\mu$SR data. We can extract the superfluid density $\rho(T)=\lambda^{2}(0)/\lambda^{2}(T)$ independent of the assumptions made in Eq. (18) if we assume the carrier effective mass is constant.

LaCuSb2 comprises multiple Fermi surfaces. Besides several small pockets around the $\Gamma$-point, there are a pair of large quasi-$2$D diamond-shaped Fermi surfaces. Along with the pockets around the $X$-point, these arise from the topological bands with Dirac nodal lines protected by the nonsymmorphic symmetries of the space group. To better demonstrate the band topology and study the consequences for superconductivity, we construct an eight-band tight-binding model that captures essential features of the Fermi surfaces of topological bands.

According to the first principles calculations, the topological bands mainly consists of $p_{x,y}$-orbitals of $\ch{Sb}$, which allows us to focus on a $2$D $\ch{Sb}$ square net layer. The $\ch{Sb}$ square net is geometrically a square lattice where $\ch{Sb}$ atoms occupy the lattice sites. However, the unit cell of the square net (defined by the base vectors $\mathbf{a}$ and $\mathbf{b}$ of $\ch{LaCuSb2}$) is the doubled unit cell of the square lattice. Therefore, in each unit cell of $\ch{Sb}$ square net, there are two $\ch{Sb}$ atoms. These two sites per unit cell can be artificially distinguished by slightly displacing the atoms along $\pm{z}$-directions. It should be emphasized that, in reality, the $\ch{Sb}$ atoms forming the square net are geometrically all in a plane, and the difference comes from the chemical environment of the neighbouring layers above and below. Nevertheless, the artificial displacements reduce the accidental symmetries of a square lattice to the true symmetries of the layer group P4/nmm inherited from the space group bearing the same name.

Extended Fig. 10a,b show the high-symmetry points and paths in the Brillouin zone for space group $P4/nmm$, along with those at $k_{z}=0$, respectively. Extended Fig. 10c shows the unit cell of $\ch{Sb}$ square net. The orange circles represent the $\ch{Sb}$ atoms. The inversion center is indicated by the red circle in between the $\ch{Sb}$ sites, and the two-fold rotational symmetry along $z$-axis is centered on $\ch{Sb}$ sites.

After putting $p_{x,y}$-orbitals on each $\ch{Sb}$ site, and labeling the two sublattices by $r=\{A,B\}$, we now describe our tight-binding model. The nearest neighbour hoppings and second nearest neighbour hoppings are parameterized by Slater-Koster type parameters. Notice that, for nearest neighbour hoppings, the hoppings are along diagonals, so the orbital basis are rotated by $45^{\circ}$ as

The Hamiltonian for the nearest neighbour hoppings can be written as

where the sum of $\mathbf{R}$ runs through all the unit cells, and $\bm{\delta}$ points towards the nearest neighbouring sites. The matrices $\tau_{i}^{(r)}$ are Pauli matrices with superscript specifying the physical degrees of freedom the Pauli matrices act upon. Here $(r)$ indicates the sublattice index. We will use $\tau^{(\sigma)}_{i}$ for spins and $\tau^{(\Delta)}_{i}$ for particle-hole doubling in the Bogoliubov-de Gennes formalism of superconductivity.

The Hamiltonian for the second nearest neighbour hoppings can be written in a similar fashion as

where $\bm{\delta}^{\prime}$ points towards the second nearest neighbour sites.

We also include the orbital energies and the chemical potential for the sake of describing doping and discussing superconductivity. These are simply diagonal terms

where we assume the orbital energies $\epsilon_{r,p_{\alpha}}=\epsilon_{p}$ are the same for both sublattices required by the crystal symmetry (e.g., the inversion symmetry).

To determine the tight-binding parameters, we fit our model band structure to the first principles calculations. The full band structure is shown in Extended Fig. 10d. We find $t_{1}=3.375$, $t^{\prime}_{1}=0.875$, $t_{2}=0.125$, $t^{\prime}_{2}=0.125$, and $\epsilon_{p}-\mu=-0.625$ all in unit of eV. Extended Fig. 10e shows our model bands that give rise to the Dirac nodal lines.

The band structure in Extended Fig. 10e exhibits several band crossings near the $X$-point. Bands that are degenerate at the $X$-point remain doubly degenerate along $X$-$M$, which gives rise to a dispersive Dirac nodal line. We also note that the crossings between $\Gamma$-$X$ and $\Gamma$-$M$ have the same origin. They are, in fact, part of the diamond-shaped Dirac nodal line intersecting with the high symmetry planes. The Dirac nodal lines are protected by the non-symmorphic crystal symmetry, as we will demonstrate below.

The key symmetry element that protects the nodal lines is the glide mirror plane $g=\{R=M_{z}|\mathbf{t}=[1/2,1/2,0]\}$. For any $\mathbf{k}$ in $k_{z}=0$ plane, $g\mathbf{k}=\mathbf{k}$, therefore, we may choose the Bloch waves to be eigenstates of $g$ as well, i.e. $g|u_{\mathbf{n,k}}\rangle=\lambda_{n,\mathbf{k}}|u_{\mathbf{n,k}}\rangle,\lambda_{n,\mathbf{k}}=\lambda_{n}\exp(\mathrm{i}\mathbf{k}\cdot\mathbf{t})$. Since $R^{2}=1$, $\lambda_{n}=\pm 1$. Extended Fig. 10f shows band cuts along $\mathbf{G}_{x}=2\pi/a\hat{x}$. The solid curves are bands cut through the $\Gamma$-point, and the dashed curves are bands cut through $\Gamma+\pi/2\hat{y}$. Band eigenvalues are calculated according to our tight-binding Hamiltonian. Bands with eigenvalue $\lambda=-1$ are in red and $\lambda=+1$ in blue. Furthermore, the time-reversal and the inversion symmetries ensure two-fold degeneracy of each band due to spins. Therefore, the symmetry protected crossings are Dirac nodes.

It bears emphasizing that the nodes on BZ boundaries and those within the BZ are of different types. The TR symmetry $\Theta$ together with the non-symmorphic symmetry $g$ protects the nodal line along $X$-$M$, which can be understood as Kramers degeneracy with respect to the antiunitary operator $\tilde{\Theta}=g\Theta$, with $\tilde{\Theta}^{2}=-1$. Therefore, these crossings, located at BZ boundaries, are like type-II Dirac node.

Regarding four bands as two pairs of intertwined bands crossing at BZ boundaries, the band crossings between the pairs are also protected by $g$. However, these crossings are less robust compared to the previous case, in the sense that when the pairs of bands are deformed, the inter-pair crossings can move and even annihilate pairwise. Extended Fig. 10f shows two cuts along $\mathbf{G}_{x}=2\pi/a\hat{x}$ through the $\Gamma$-point and $\Gamma+\pi/2\hat{y}$. Note that the inter-pair crossing moves towards $k_{x}=0$. These crossings can be finally annihilated and give rise to a closed diamond-shaped nodal line inside BZ. Therefore, unlike the crossings at BZ boundaries, the inter-pair crossings are like type-I Dirac nodes.

In the presence of the spin-orbit couplings (SOC), the diamond-shaped nodal line will generically be gapped (also hybridized with other bands). However, as the SOC gaps are small and more than $200\text{meV}$ below the Fermi energy, our simple tight-binding model still gives a good description of the spin-orbital textures of the Fermi surface that are relevant for superconductivity and transport measurements. The SOC effect on the nodal line along $X$-$M$ also vanishes to leading order. Therefore we neglect SOC in our tight-binding model, and the model band Hamiltonian is

We now assume the mean field approximation applies to $\ch{LaCuSb_{2}}$ and discuss the superconducting pairings according to Bogoliubov-de Gennes effective Hamiltonian. We write the generic superconducting paring terms for Cooper pair with zero center-of-mass momentum as

where $\hat{\Delta}^{T}(-\mathbf{k})=-\hat{\Delta}(\mathbf{k})$ is required by the anticommutation relation for the fermionic operators. The transpose $T$ acts on all the indices including spins, sublattices, and orbitals. We choose the basis of matrices

Therefore, there are $28$ possible terms. From current experimental results, the superconductivity is consistent with singlet pairing, so we can fix the spin part to be $\mathrm{i}\tau_{y}^{(\sigma)}$. Transforming Eq. 24 back to real space, we assume that the leading pairing terms come from the local on-site attractive interactions, which further fixes pairing in the sublattice space to be $\tau_{0}^{(r)}$. There remain three possibilities for orbital degrees of freedom i.e., $\tau_{0}^{(\alpha)}$, $\tau_{1}^{(\alpha)}$, $\tau_{3}^{(\alpha)}$. Therefore, we have