Nodal multigap superconductivity in the anisotropic iron-based compound RbCa2Fe4As4F2

The 12442 compounds are a recently discovered family of iron-based superconductors, that share several features with the cuprates due to their strongly anisotropic structure, but are so far poorly understood. Here, we report on the gap structure and anisotropy of RbCa2(Fe1-xNix)4As4F2 single crystals, investigated by a combination of directional point-contact Andreev-reflection spectroscopy and coplanar waveguide resonator measurements. Two gaps were identified, with clear signatures of d-wave-like nodal structures which persist upon Ni doping, well described by a two-band d-d state with symmetry-imposed nodes. A large London penetration depth anisotropy was revealed, weakly dependent on temperature and fully compatible with the d-d model.


INTRODUCTION
The so-called 12442 family of iron-based superconductors (IBSs) is the youngest addition to this class of materials [1]. It belongs to the generalized 122 structure, in which superconducting (SC) Fe 2 As 2 planes are alternated by spacer layers. In the proper 122 family (such as the AFe 2 As 2 compounds) the spacer layer is made only of A atoms, while in the 1144 systems two 122 layers are alternately stacked. The 12442 family is similar to the 1144 but one of the alternate spacers is not composed of one atomic species but by a larger insulating Ca 2 F 2 layer [2] (Fig. 1a). The 12442 structure is therefore less trivial and clearly more anisotropic than those of the rest of the IBSs and, in this respect, more resembling that of cuprates such as Bi 2 Sr 2 CaCu 2 O 8+δ and YBa 2 Cu 3 O 7−δ (YBCO) [3][4][5]. This intuitive and intriguing similarity acquires even more charm in light of the several (although not uncontroversial) experimental observations suggesting the presence of a d-wave-like gap opening on the complex, multiband, Fermi surface [6][7][8]. It seems therefore that the 12442 family might be a host of several coexisting features that are key for the understanding of the physics of unconventional superconductors. The Fermi surface of these materials is composed of six to eight sheets (observed by angle-resolved photoemission spectroscopy, ARPES, [9] or proposed by ab-initio calculations [10]), with large hole pockets at the center * These authors contributed equally. † Corresponding author: gianluca.ghigo@polito.it of the Brillouin zone and a very small electron pocket at the M point. This situation of almost vanishing electron pockets is close to what is found at the end of the doping series in 122 IBSs (KFe 2 As 2 and KFe 2 Se 2 at the opposite extrema) in which the typical s ± symmetry of the order parameter seems not to hold, with the emergence of nodal behavior [11][12][13][14] due to the proximity, in terms of energy, of the two instabilities.
Features compatible with the presence of a d-wave gap on one of the bands, possibly coexisting with a fully gapped state on other bands, have in fact been observed by µSR superfluid density measurements (on K-and Cs-12442 polycrystals [6,7]) and low temperature specific heat (on K-12442 single crystals [8]). Recently, schemes that explain the coexistence of gapless and gapped features were proposed and shown to be suitable candidates for multiband systems [15,16]. Conversely, no nodal gaps were observed by heat transport [17] and optical spectroscopy [18] in Cs-12442 samples, and by ARPES [9], lower critical field measurements [19] and tunneling spectroscopy (STS) [20] in K-12442 single crystals. STS measurements revealed a single (contrasting all other reports) s-wave gap, with an intraband coupling mechanism that would be incompatible with spin-fluctuation-mediated superconductivity, for which clear signatures were observed by neutron spin resonance measurements [21,22]. This large set of contrasting reports highlights the complexity of these systems, and the need of joint investigations of the same high-quality samples by means of different techniques in order to shed light on this complicated matter.
In this work, motivated by these considerations, we report on a combined study of the gap symmetry and superfluid density of RbCa 2 (Fe 1−x Ni x ) 4 As 4 F 2 single crystals employing directional point-contact Andreevreflection spectroscopy (PCARS) and a coplanar waveguide resonator (CPWR) technique on samples with increasing Ni doping. The study of doped samples gives the possibility to investigate the persistence of nodal features in the presence of scattering provided by structural disorder, allowing us to speculate and possibly to discern whether they are symmetry-imposed or accidental [24,25].

Synthesis and basic characterization
Single crystals of RbCa 2 (Fe 1−x Ni x ) 4 As 4 F 2 with x = 0, 0.03 and 0.05 were grown by the self-flux method using RbAs [3,4,26]. The high quality of the investigated single crystals was assessed through X-ray diffraction, energy dispersive X-ray (EDX) spectroscopy and transport measurements, as shown in Fig. 1b-c and in Supplementary Fig. 1. The EDX analysis is performed routinely after crystal growth to determine the exact composition; a typical spectrum and the compositional study are reported in Ref. 3.
Particularly important is to notice that the resistivity at T c increases with increasing Ni content (Fig. 1c), indicating an enhanced scattering in the normal state. However, the superconducting transitions remain very sharp despite the strong decrease of critical temperature (inset of Fig. 1c). This highlights the primary role of Ni substitution in shifting the chemical potential and changing the electron doping level, rather than enhancing the pairbreaking scattering in the superconducting state (that would broaden the transition). This effect was pointed out also by Hall coefficient measurements in Ref. 3, although the T c dependency on the Ni content was there interpreted as an indirect indication of the s ± symmetry of the order parameter. Actually, due to this double role of Ni doping, direct measurements are necessary to investigate fundamental details of the gaps, as we will show in the following.

Directional point-contact Andreev-reflection spectroscopy
The number and symmetry of the SC energy gaps were determined via directional PCARS [27][28][29][30][31] in crystals with x = 0 and x = 0.05. Fig. 2a-d reports some examples of normalized differential conductance spectra (dI/dV , symbols) acquired by injecting the current along either the ab plane or the c axis (see Methods, Supplementary Methods, and Supplementary Fig. 4 for details). As in all the Fe-based compounds, in addition to the structures associated to the energy gap(s), the spectra present higher-energy shoulders (falling in the grey regions) that can instead be associated to the strong coupling between carriers and the bosonic mode responsible for the pairing [32][33][34]. In both the undoped and 5% Ni-doped crystals, 100% of the spectra display zero-bias peaks or cusps that are strongly suggestive of a gap with vertical node lines and a change of sign on the same Fermi surface [35]. This agrees with all the evidences of nodal gaps in other compounds of the same family (Cs-12442 and K-12442), but is completely incompatible with a single s-wave gap picture (as observed in K-12442 by STM [20]). The suppression of the critical temperature T c by ∼ 10 K in the 5% Nidoped samples with respect to the undoped compound results in a decrease in the gap amplitudes but, also, in a reduction of the energy at which the boson structures set in; as a result, the latter fall very close to the edge of the large gap, which makes it more difficult to disentangle the relevant contribution. To extract quantitative information on the gaps, it is necessary to fit the spectra with a suitable model for Andreev reflection in the presence of a nodal gap. We used the one reported in Ref. 36, which perfectly fits the case of ab-plane contacts, and its generalization to the case of c-axis contacts in a highly anisotropic material [37]. Some details about this model and its parameters are given in the Methods and Supplementary Methods. Since it is based on the BCS, weak-coupling assumption  of energy-independent order parameters, this model is unable to reproduce the bosonic structures; therefore, the fit has been made by disregarding them and focusing on the central part of the curves, that instead contains all the information about the SC energy gaps; the validity of this approach has been demonstrated in Refs. 33 and 34. As shown in Supplementary Fig. 7c-d, if one tries to fit these structures as if they were due to a gap, very high (and probably non-physical) values of the gap ratio κ = 2∆/k B T c are obtained.

Gap structure
As for the gap symmetry, a guess has to be made before trying to fit the spectra. Since specific-heat [38] and µSR [6,7] measurements in similar compounds were interpreted as suggesting the coexistence of a nodal gap together with a nodeless one, we initially tried to fit the spectra by using one s-wave gap and one nodal gap, for simplicity assumed to have a d-wave symmetry, reasonably residing on different Fermi surface sheets. We found that all curves for undoped and 5% Ni doped crystals, both in ab-plane and c-axis contacts, can be satisfactorily fitted by this s−d model, with a larger s−wave gap ∆ 1 and a smaller d−wave gap ∆ 2 , but the contribution of the s gap is barely visible and the gap values extracted from different contacts are rather scattered (for x = 0, ∆ 1 = 4.3 − 6.6 meV and ∆ 2 = 1.9 − 3.5 meV) and, within the experimental uncertainty, look to be directly proportional to T c (see Supplementary Fig. 8).
Despite this ostensible success, the existence of a nodal and a nodeless gap in the same system poses some questions from the theoretical point of view, and may not be the most reasonable explanation of the observed results. In particular, doubts might arise regarding the compatibility of this gap structure with the underlying I4/mmm point group symmetry and Fermi surface structure. Moreover, the electron-electron interaction in the IBSs is mainly provided by repulsive antiferromagnetic spin fluctuations (AFM-SF), that provide coupling over regions of the Fermi surface with opposite sign of the SC gaps. Thus, a single s-wave gap cannot exist within this picture, but two isotropic gaps with opposite sign (s ± state) are possible. Therefore, in an effective two-band model with s-and d-wave gaps, one needs to admit that two s-wave gaps with opposite sign open up on either band (∆ s,1 and ∆ s,2 ); as for the d-wave gap, it can be considered to exist in both the bands for the sake of generality: ∆ d,1 = ∆ d,2 . In other words, each band should host a s + d gap. This gap displays nodes only if the dwave component is larger than the s-wave one, in contrast with the experimental results of the fit (see also Supplementary Fig. 8). Therefore, a two-band s−d model is likely to be, at least, oversimplified for this system.
Since the PCARS spectra (as well as the CPWR measurements, as discussed below) unambiguously point to the existence of nodes in the gap, it makes sense to fit them by assuming either a single nodal gap (which however does not allow properly fitting all the curves) or two nodal gaps. For simplicity, we chose a d x 2 −y 2 -wave symmetry for both of them, with the nodal direction at an angle π/4 with respect to the k x axis in the reciprocal plane (see Methods). Such a two-band d−d picture automatically: i) ensures the existence of nodes, ii) respects the overall symmetry of the system, and iii) does not pose problems of coexistence of the different gaps even in the simplest two-band model. In the following, we will show that it is also fully compatible with our PCARS and CPWR results. Fig. 2g depicts a possible realization of this d−d picture in the case where the Fermi surface consists of a holelike sheet at the centre of the Brillouin zone, and of an electronlike pocket at the corner: a schematic simplification of the 12442 case. In general, the quality of the fits of the PCARS spectra within the d−d model is comparable to that obtained with the s−d model. However, in the cases where the zero-bias peak is the more relevant feature, the d−d fit is by far the most natural choice and reproduces very easily the experimental data. In the case of ab-plane contacts, for example, the change of sign in the gap is expected to give rise to a constructive interference of electronlike and holelike quasiparticles (ELQ and HLQ in the following) whenever the current is injected at an angle α = 0 with respect to the k x axis, with a maximum probability for α = π/4. This phenomenon explains both the peakshoulders feature of the spectrum in Fig. 2a (α π/8) and the high zero-bias peak of Fig. 2c (α = π/4). In c-axis contacts this interference cannot take place since ELQ and HLQ always feel gaps with the same sign, but the existence of nodes makes unpaired quasiparticles to be available for conduction even at very low energy. This naturally explains the zero-bias maxima in Figs. 2b,d; the latter looks narrower mainly because of the smaller amplitude of the larger gap.
As a further check of the correctness of our approach, we report in Fig. 2e the temperature dependence of the spectrum shown in Fig. 2d . The experimental conductance curves (symbols) are compared to the relevant fitting curves (solid lines). We chose to keep all the parameters fixed to their low-temperature values, and to tune only the gap amplitudes. Despite this very strict constraint, the theoretical curve follow very well the evolution of the central part of the spectra. Moreover, the gap amplitudes approximately follow a BCS trend, as shown in Fig. 2f.
When the whole series of measurements is considered, it turns out that the gap amplitudes extracted from the d − d fit are less scattered than in the s−d case: for x = 0, ∆ 1 ranges between 4.9 and 6.6 meV and ∆ 2 from 1.5 to 2.5 meV. Fig. 4a reports the gap amplitudes as a function of the T c of the contacts; the shaded regions are delimited by lines of constants gap ratio. For the larger gap, 2∆ 2 /k B T c ranges from 3.3 and 5.0 and is centered at about 4.2, the value predicted by the BCS theory in the d-wave case. The smaller gap has instead a gap ratio between 1 and 2, much smaller than the BCS value. Interestingly, the c-axis and ab-plane gap amplitudes follow the same trend, with no evidence of a dependence on the direction of current injection.

Coplanar waveguide resonator analysis
As already mentioned, a combination of multiple bands and unconventional gap symmetry in 12442 compounds was also suggested on the basis of the low-temperature linear feature in the superfluid density observed by µSR in K-and Cs-12442 polycrystals [6,7]. Fig. 3 shows the temperature dependence of the normalized superfluid density ρ s in our undoped (x = 0) and Ni-doped (x = 0.03 and 0.05) single crystals, determined by starting from the CPWR measurement of the London penetration depth λ, i.e. ρ s = [λ(T = 0)/λ] 2 (see Supplementary Fig. 3 and Refs. [39][40][41][42]. A linear temperature dependence of ρ s is, indeed, evidenced also in our samples, in the experimentally accessible temperature range. The absence of kinks in the superfluid density curve is typical of single band superconductors and of multigap systems with mainly interband coupling, as the IBSs [40]. The temperature dependence of the superconducting gaps measured by PCARS (Fig. 2f) indeed shows that both gaps slowly decrease and close -together-at T c , therefore a smooth ρ s (T ) curve is expected. With respect to Refs. 6 and 7, our measurements display a much less pronounced curvature at high T , and a wider T range where the behavior is quite linear, particularly evident in Ni-doped samples. These differences might be due to differences among the systems (Cs-, K-and Rb-based) or to the single-vs. poly-crystalline type of samples. Notably, our ρ s data can be fitted, for all doping levels, with a BCS models involving two d-wave gaps, as done in Refs. 6 Fig.4a-b) that are in very good agreement with those found experimentally by PCARS. This is a further support to the validity of the d−d picture for Rb-12442, but is not in contrast with the µSR experiments on Csand K-12442 [6,7], since in those systems the s−d and the d−d fits are also very similar, possibly due to the contribution of several bands to the overall density. It should be noted that, although also a single d-wave gap fit would give a reasonable result, the agreement with the PCARS data (unambiguously pointing to the presence of at least two gaps) would be lost. The fact that the same nodal features are observed, both in PCARS and CPWR measurements, up to 5% Ni doping is a strong indication that Ni doping, up to this level, neither affects the number and the symmetry of the energy gaps, nor lifts the nodal lines, suggesting that they might be symmetryimposed and not accidental [24]. In principle, however, this could also mean that the amount of scattering provided by this doping level is just not sufficient to lift the nodes. In this respect, we recall the above discussion, regarding the weak role of Ni atoms as pair-breaking scattering center. This is further evidenced by the fact that there is no clear change in the ρ s (T ) behavior at low T, that is expected to change from linear to quadratic for a d-wave superconductor in the dirty limit. Additional studies, with disorder introduced via ion or electron irradiation [43][44][45][46] in undoped samples are desirable and ongoing.

Anisotropy
Starting from the CPWR measurements of the London penetration depth of crystals with different shape factor, the penetration depth anisotropy parameter γ λ = λ c /λ ab can be extracted [49,51] (see Methods and Supplementary Methods). The measured values, shown in the inset to Fig. 3, are much larger than what is generally found for other IBS of the generalized 122 family, but in line with those reported for the upper critical field in the same system [3] and in K-12442 [52] and close to the typical values found in cuprates [50], as shown in Fig. 4c. Moreover, γ λ is found to decrease with increasing Ni substitution (Fig. 4b), as expected due to the presence of Ni impurities, that act also as weak isotropic scattering centers [53]. However, at 5% Ni content the anisotropy parameter value is still about twice what is found for undoped CaKFe 4 As 4 [49]. The temperature dependence of γ λ is quite flat, with a slight upturn at higher temperatures.  Fig. 3 (left scale). c, Anisotropy of the London penetration depth (measured at T 5 − 6 K) in the undoped Rb-12442 single crystal compared to that of MgB2 [47], other IBSs of the generalized 122 family [48,49], and YBCO [50].
For long, a temperature dependence of this quantity was interpreted as a signature of the multiband nature of the system. However, recent calculations [54] have shown that in multiband systems with spheroidal Fermi surfaces any temperature dependence (increasing, decreasing or constant) is achievable depending on several microscopic parameters. In short, these calculations indicate that: i) for a single d-wave gap with vertical node lines, of the form ∆ cos(2φ) (φ being the azimuthal angle) γ λ should be temperature independent, as in the single-band s-wave case; ii) for the multiband case, even two s-wave gaps can give rise to any temperature dependence of γ λ , depending on the band anisotropies. Following the same approach as in Ref. 54, it is possible to generalize the calculations to a two d-wave gap system on spheroidal Fermi surface (details in the Supplementary Methods). One finds that the results discussed in Ref. 54 for the two-band s-wave case are strictly valid also for the two d-wave gaps case. Therefore, the weak temperature dependence found for γ λ is fully compatible with the d−d gap structure considered. It also follows, from the constancy of the gap ratios measured along the ab and c directions for undoped and doped samples (compared to the strong decrease of γ λ with doping) that the main factors influencing the value of γ λ could be the anisotropy of the Fermi velocity and the 2D (cylindrical as opposed to spheroidal) topology of the Fermi surfaces. The decrease of γ λ with doping could be determined by isotropic scattering provided by the impurity atoms as in the 1144 system [49], or by a change of values of the Fermi velocities that seems reasonable when the Fermi level is shifted by doping.

Concluding remarks
In conclusion, we have brought evidence for the existence of two nodal (possibly d-wave) gaps in Rb-12442 by using two techniques, PCARS and CPWR, that allow investigating the symmetry of the order parameter from different perspectives and yielded concording results. The shape of the PCARS spectra, as well as the low-temperature linear behaviour of the superfluid density measured by CPWR, clearly indicate the presence of nodes in the gap that persist upon Ni doping, suggesting that the nodes might be symmetry imposed. Moreover, both the PCARS spectra and the ρ s (T ) curve can be fitted very well within a d−d model (compatible with the Fermiology and symmetry of the system) yielding gap amplitudes in very good agreement. These facts, together with the unusually high values of the London penetration depth anisotropy, seem to tighten the similarity between the 12442 family of IBSs and the copper oxides, which makes the former an ideal platform to investigate the nature of superconductivity in both these classes of unconventional superconductors.

PCARS measurements
The superconducting gap structure was investigated via pointcontact Andreev reflection spectroscopy (PCARS), which entails the measurement of the differential conductance (dI/dV ) of a pointlike contact between the superconductor (S) and a normal metal (N). Here, the contacts were fabricated via the so-called soft technique [27,[29][30][31], in which a thin Au wire is stretched over the crystal, until it touches the surface in a single point, is (optionally) mechanically stabilized by drop-casting a small (≈ 50 µm) droplet of Ag glue, and acts as a parallel of nanoscopic N/S contacts. Thanks to the platelet-like shape of the crystals, it was possible to control the direction of (main) current injection, by making the point contacts either on the topmost surface (thus injecting the current along the c axis) or on the side (thus injecting the current along the ab planes). We refer to these contacts as c-axis and ab-plane contacts both in the following and in the Results and Discussion section. Actually, to ensure a good directionality of the current injections, c-axis contact were made in flat, mirrorlike regions of the top surface, excluding growth terraces; as for ab-plane contacts, the very thin, but smooth and flat side surfaces of the crystals allowed us to assume a small, if any, contribution of conduction along the c axis. A picture of contacts along either direction is reported in Supplementary Fig. 4.
Proper PCARS measurement conditions require the contact to be smaller than the electronic mean free path, resulting in ballistic conduction with no Joule dissipation in the contact region, and a contact resistance largely exceeding the resistance of the normal bank [27,33]. In these conditions -and in the absence of a surface insulating layer -electronic conduction across the S/N interface is dominated by Andreev reflection [55], and a quantitative analysis can be carried out by fitting the experimental dI/dV vs V spectra with suitable models.
Prior to fitting, the experimental dI/dV curves need to be normalized, i.e. divided by the normal-state differential conductance curve [27,30,31]. In principle, the dI/dV spectrum recorded at a given T < Tc should be divided by the normal-state conductance at the same temperature, obtained by suppressing superconductivity with, e.g., a magnetic field. However, when the upper critical field is very high (as in most IBS) a common practice is to use the normal-state conductance recorded just above Tc. The small thickness of the crystals (a few micrometers) however leads to a rather large normal-state resistance that, in the pseudo-four-probe configuration used for PCARS, is in series with the contact resistance and gives rise to a downward shift and a horizontal stretching of the conductance curves with respect to those at low temperature [34,56]. To recover the presumable normal-state conductance at low temperature, we rescaled the dI/dV at T > Tc by subtracting the contribution of this so-called spreading resistance Rs. The value of Rs for each contact was typically selected as the one that allowed to match the tails of the superconducting and normal dI/dV spectra, but the choice remains somewhat arbitrary. Therefore, multiple reasonable values of the spreading resistance were used to obtain the normalized spectra, and the fitting procedure repeated each times. The resulting values of the fit parameters were then averaged and their maximum differences taken as the uncertainty of the procedure.
The normalized PCARS spectra were fitted to a two-band model, using either two d-wave gaps (as discussed in the Results and Discussion section), or one s-wave gap and a d-wave gap (summarized in Supplementary Figs. 7,8). The quasi-cylindrical shape of the Fermi surfaces is taken into account while performing the integration over all possible directions of incident electrons. For c-axis contacts, we integrated over a cylindrical belt about the basal plane kz = 0, as explained in Ref. 37. For ab-plane contacts, we integrated in the whole half-plane kx > 0 as explained in Ref. 36. In either case, the model is actually the same and takes into account the interference effects that occur when HLQ and ELQ (injected with different wavevectors) feel order parameters of opposite sign. This actually occurs only in ab-plane contacts, when the normal to the N/S interface makes a non-zero angle α with the (antinodal) kx direction, and can give rise to zero-energy peaks in the spectra (see Supplementary Methods for more details). In both cases the model contains as adjustable parameters the gap amplitudes ∆ 1 and ∆ 2 , the broadening parameters Γ 1 and Γ 2 , the barrier parameters Z 1 and Z 2 , and the relative weight of the large-gap band in the conductance w 1 . In the d−d case, α was supposed to be the same for both gaps. Despite the large number of parameters, the values they acquire in the fits is not arbitrary because each parameter has different and specific effects on the shape of the spectrum, as extensively discussed in Refs. 27, 29-31. In spectra with a low signal, as in Fig.2a, there is some interplay between Γ and ∆ which gives rise to an increased uncertainty on the gap amplitude, here included in the error bars. In all cases, the high-energy shoulders highlighted by the shaded gray bands in Fig. 2a-d and Supplementary Fig. 7 must be excluded by the fitting procedure, since they are strong-coupling features that are not captured by models based on the weak-coupling BCS theory [32][33][34]. Indeed, if these structures are fitted as if they were due to an energy gap, the resulting gap amplitude is unphysically large (see Supplementary Methods for details).

CPWR measurements
The superfluid density, surface impedance and London penetration depth anisotropy were measured by means of a coplanar waveguide resonator (CPWR) technique that has proven suitable to study IBS crystals [41,43,57]. The measurement device consists of an YBa 2 Cu 3 O 7−δ resonator in the shape of a coplanar waveguide, to which the sample is coupled. Resonance curves are recorded with a vector network analyzer as a function of temperature, making it possible to track resonant frequency shifts and variations of the unloaded quality factor across the superconducting transition of the investigated samples. After a calibration procedure [41,58] is performed, this gives access to the absolute value of the penetration depth and its temperature dependence, as well as to the surface impedance temperature dependence [42,59]. Due to the fields distribution at the sample position (vanishing electric field and rf magnetic field parallel to the ab-planes of the sample) and to the finite size of the crystal, the measured penetration depth λ is a combination of the main anisotropic components λ ab and λc. This combination depends on the geometry of the sample under consideration, and for this reason both components can be retrieved by measuring samples with significantly different aspect ratios fs = c · (1/a + 1/b) [49]. Accordingly, and for samples with λ ab < c and λc < a, b (where 2c, 2a, 2b are respectively the thickness, width and length of the samples), the fraction of penetrated volume can be estimated as λ ab /c + λc/a + λc/b. Thus, the measured penetration depth can be expressed as λ = λ ab + fs · λc.

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Methods. Addi-tional data related to this paper may be requested from the authors.