Multiple Andreev reflections effect spectroscopy of LiFeAs single crystals: three superconducting order parameters and their temperature evolution

The structure of the superconducting order parameter of LiFeAs is studied by incoherent multiple Andreev reflections effect (IMARE) spectroscopy. The high transparent superconductor–thin normal metal–superconductor (SnS) contacts are created by a planar “break-junction” technique. Below Tc≈17.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c \approx 17.5$$\end{document} K, the obtained I(V) and dI(V)/dV characteristics of SnS junctions show a presence of at least three bulk superconducting order parameters in LiFeAs. We directly determine the magnitudes, characteristic ratios, and temperature dependences of the superconducting gaps and discuss their symmetry. Three-gap superconductivity with possible momentum dependence (anisotropy) in LiFeAs. Temperature dependences of the superconducting gaps indicate a moderate interband coupling. Superconducting energy parameters scale with Tc under minor lithium deficiency. Three-gap superconductivity with possible momentum dependence (anisotropy) in LiFeAs. Temperature dependences of the superconducting gaps indicate a moderate interband coupling. Superconducting energy parameters scale with Tc under minor lithium deficiency.


Introduction
Layered LiFeAs [1] belongs to one the most intriguing families of the pnictide superconductors, with the socalled 111 structure type. Its crystal structure contains superconducting FeAs blocks alternating with Li planes along the crystallographic c direction. Contrary to the majority of the iron-based pnictides (for a review, see [2,3], LiFeAs shows no magnetism [4], whereas its superconducting properties are optimal in the stoichiometric state with the maximum critical temperature T c ≈ 18 K. Any substitution, lithium deficiency, or applied pressure rapidly decrease the T c to zero [5][6][7][8], so the question of how to increase the T c in LiFeAs is still open. The Fermi level in LiFeAs is crossed by several bands forming two concentric hole barrels near the Γ point of the first Brillouin zone, and two electron barrels around the M point (in terms of 2-Fe unit cell), as shown theoretically [9,10] and experimentally using angle-resolved photoemission spectroscopy (ARPES) [4,8,[11][12][13]. The phase diagram, superconducting properties, and current advances of their studies are reviewed recently in [14].
Due to the presence of active alkali metal, LiFeAs appears extremely sensitive to even trace amounts of oxygen, water vapors, or nitrogen, which results in quick degradation of its superconducting properties and the T c dropping to zero in 5-10 min exposure in open air. Therefore, all the experimental steps should be done in a protective "dried" atmosphere. Unfortunately, this chemical activity strongly complicates any probing of LiFeAs and is responsible for a scarcity of the available experimental data.
ARPES probes [4,12,13] revealed that below T c the largest superconducting gap develops at the inner (shallow) hole barrel around the Γ point, the medium-size gap-at electron Fermi surface sheets, whereas the small gap develops at the outer hole barrel. Theoretically calculated superconducting gap distribution over the Fermi surface in the framework of a so-called s ++ pairing through orbital fluctuations additionally combined with s ± pairing through spin fluctuations [15,16] showed a good agreement with the ARPES data [4,12,13]. On the other hand, such gap structure could be reproduced within the s ± approach solely by accounting an orbital selectivity (different correlation strengths in the bands formed by different orbitals; in particular, Cooper pairing strength) [17]. Generally, the majority of theoretical [15][16][17][18][19] and experimental studies [4,[12][13][14][20][21][22] consider a substantial anisotropy of the superconducting properties of LiFeAs in the ab-plane.
Here we continue the IMARE studies [20,21] of LiFeAs single crystals. We consider the current-voltage characteristics (CVC) and the dynamic conductance dI(V)/dV spectra of planar SnS-Andreev contacts with incoherent transport ("long") and almost optimal local critical temperature T c ≈ 17.5 K. We show a reproducible I(V) and dI(V)/dV data that could be interpreted in the framework of a three gap structure, and directly determine the magnitudes, characteristic ratios, and temperature dependences of the superconducting order parameters, and discuss their symmetry. By comparing the current results with our earlier data [14,21], we unambiguously show the characteristic ratios of all the energy gap parameters remain almost constant within the range T c ≈ 15.5-17.5 K.
The paper is organized as follows. In the next section, we describe the IMARE spectroscopy and the planar breakjunction technique. Section 3 presents the experimental I(V) and dI(V)/dV of the obtained SnS junctions and details their features. In Sect. 4, we discuss the interpretation of the experimental data and summarize the determined parameters of the superconducting state of LiFeAs.

Experimental details
LiFeAs single crystals in the form of thin rectangular plates with dimension up to 7 mm were grown using a self-flux technique. The details of crystal growth and characterization are presented in [5]. The bulk critical temperature T c ≈ 18 K was confirmed by transport and magnetic measurements.
In order to directly determine the superconducting order parameters of LiFeAs and their temperature dependences, we used an incoherent multiple Andreev reflections effect (IMARE) spectroscopy. Generally, MAR effect [35][36][37] occurs in ballistic junctions of superconductor-thin normal metal-superconductor (SnS) type. Regardless to a coherence of the transport current, a series of dynamic conductance features called subharmonic gap structure (SGS) [35,37,38] appears at I(V) and dI(V)/dV at any temperatures up to T c . The position V n of SGS directly relates to the superconducting gap magnitude Δ(T ) as [35,37] where the subharmonic order n = 1, 2, … is the natural number. Unlike probing asymmetric NS and NIS junctions (I is insulator), where the finite temperature smears a shape of features located at eV < Δ , no fitting of dI(V)/ dV is needed in case of SnS contact till T c , which facilitates a precise measurement of temperature dependence of the gap [37,38]. In "long" SnS junction (with the characteristic dimension d exceeding the Josephson coupling length) and the barrier strength Z ≲ 0.5 1 at temperatures below T c , incoherent Andreev transport [35,37] causes no Josephson supercurrent branch, but an excess current at the whole bias voltage range, which drastically rises at eV → 0 , thus forming so-called foot area. The corresponding dynamic conductance spectrum shows an increase in zero-bias conductance G ZBC (strictly at eV = 0 ) as compared to the normal one G N (being the conductance at eV ≫ 2Δ ). For semiballistic SnS junction with l E ∕d ≃ 1 − 2 in a fully clean limit Z = 0 , expected are the fundamental ( n = 1 ) harmonic and a lower-amplitude subharmonic with n = 2 [37,38], whereas the features with the greater n's are suppressed due to the finite inelastic scattering length l E and finite corresponding scattering rate ℏ∕2 E . According to the CVC shape we assume Z = 0.2 − 0.4 for the breakjunctions obtained. In this case one can roughly estimate l E ∕d = 1.5 − 3 [38].
With temperature increase, the excess Andreev current and enhanced zero-bias conductance (ZBC) suppresses gradually; with it, the SGS features shift toward zero, and their amplitudes decrease in proportion to the concentration of Cooper pairs [38]. The local critical temperature T local c corresponds to the contact area transition to the normal state, accompanied with vanishing of all the features caused by IMARE transport. For a multiple-gap superconductor, several SGS series appear in the dI(V)/dV spectrum, corresponding to each gap.
A momentum-dependent (extended s-wave) order parameter would cause doublet-like SGS features, whereas the position of two dips forming the doublet corresponds to the maximum and minimum Cooper pair coupling energies [20,39,40]. A rough simulation of the Andreev feature in the framework of [41] was done in Fig. 4b in [39] and Fig. 5 in [40]. In the used break-junction configuration, the current always flows along the c-direction. In case of ballistic junction and almost cylindrical Fermi surface with k z ≪ k z , k y , the charge carriers would preferably keep their in-plane momentum components constant, whereas some mixing of the k z components could take place. As a result, the planar break junction provides information about namely the in-plane anisotropy of the superconducting gap.
The sample preparation and mounting was done in a "dry" argon atmosphere. The thin ab-plane oriented single crystal with dimensions about 3 × 1.5 × 0.1 mm 3 was attached to a U-shaped springy holder by four-contact pads made of liquid In-Ga solder (see Fig. 1 in review [39] for details). Then the sample was cooled down in helium atmosphere. At T = 4.2 K, under a gentle curving of the holder, the bulk sample splits along the ab-planes with a formation of two cryogenic clefts (with steps and terraces on the cryogenic surfaces) separated with a weak link, a kind of planar ScS contact (where c is constriction).
The resulting constriction forms in the bulk far from current and potential In-Ga leads, which provides a good heat sink from the constriction, and the true four-point probe. Note during the experiment the crack remains deep in the bulk, with tightly conjuncted (sliding) clefts. This prevents impurity penetration into the crack and maintains the purity of cryogenic surfaces that results in probing of bulk (not degraded) superconducting properties with a high spectroscopic resolution. In LiFeAs, similarly to the majority of Fe-based superconductors we studied, the formed constriction appears electrically equivalent to a thin layer of normal metal of high transparency (about 80%-98% ), thus providing an observation of MAR effect. This follows from the observed shape of the resulting I(V) and dI(V)/dV, which are typical for the high-transparent classical SnS-Andreev junction with incoherent transport [35,37,38].
Beside single ScS contacts, Andreev arrays with ScSc-… -S structure can be also formed in the MCBJ experiment with layered sample [39,43]. Such array resembles a natural stack of m equivalent elements (m is accidental but natural). Hence, the dI(V)/dV of the array shows the Andreev features at positions being scaled by a natural factor of m as compared to that of single SnS junction: n . The number m for each array could be unambiguously determined by comparing dI(V)/dV curves for various arrays: after scaling the bias voltage axis by m, if the equivalence of constrictions is realized, the dI(I)/dV spectrum turns to that of a single junction (for details, see the Appendixes in [43,44]). The CVC's and the dI(V)/dV spectra shown below are normalized to the single SnS junction ( V norm means V/m), whereas the vertical axis is kept original.
The arrays were typically observed by our group, with constrictions of low-transparency (Josephson regime) [39,45] or high-transparency (SnS-Andreev regime) [20,21,39,40,43,46], but the nature of the latter is still not understood. A well-known an intrinsic Josephson effect occurs in natural arrays of SIS junctions (SISI-… -S) with coherent transport and low transparency, and was observed in hightemperature cuprates and other layered superconductors [45,[47][48][49]. Similarly, one may suppose a development Summarizing the advantages of IMARE spectroscopy of mechanically-controlled planar break junctions and natural arrays, this technique provides a precise, local (within the contact area with d ≈ 10-90 nm) and high-resolution probe of the bulk superconducting order parameter, its temperature dependence and any fine structure. In our studies, the dynamic conductance spectra were measured directly by a standard low-excitation modulation technique [39]. The results obtained with this setup are insensitive to the presence of parallel ohmic conduction paths; if any path is present, the dynamic conductance curve shifts along the vertical axis, while the bias stays unchanged. Since both, the normal-state junction resistance R N = 20 − 200 Ohm and the capacitance value C n , are relatively small for the transparent classical SnS-Andreev contacts, the setup is generally insensitive to the highfrequency noise, since 1∕R N >> i C n . Figure 1 shows the typical CVC's of SnS-Andreev arrays (a) and the corresponding dynamic conductance spectra (b) shown in the same color. The I(V) characteristics have no Josephson supercurrent branch at V = 0 , but show an Andreev excess current at any bias voltages that rapidly rises at low bias voltages (see Fig. 1a): the current increases for the solid curve ♯1 at T = 4.2 K as compared to the normal-state CVC shown by dashed line). The corresponding dynamic conductance spectra in (b) panel show an enhanced ZBC peak and a set of dI(V)/dV dips. Therefore, we attribute these contacts to classical incoherent SnS-Andreev regime of high transparency with the barrier strength Z ≲ 0.5 [35,37]. Taking the inverse conductance at eV ≫ 2Δ(0) as an assessment for the normal resistances of these junctions R N1 ≈ 18 Ohm and R N2 ≈ 30 Ohm and following the rough estimation using Sharvin formula presented in [21], we get the junctions are semiballistic with about l∕d ≈ 1.7-2.1 (both l and d values are taken along the crystallographic ab-plane). Unfortunately, we cannot make similar l/d estimate along the c-direction since we do not know any studies in which the c l c value was determined for LiFeAs.

Experimental results
The arrays under consideration were formed sequentially in one and the same LiFeAs single crystal at T = 4.2 K. Under a gentle mechanical readjustment, the m = 4 junction-array (upper dI(V)/dV curve) transforms to the 5-junction array (lower spectrum). For the I(V) and dI(V)/dV curves shown in Fig. 1, the bias voltage axes were normalized by the above mentioned natural m numbers as V norm = V ∕m , thus turning to single SnS contacts. After such normalization, the positions of all the dI(V)/dV features and the width of the foot area well coincide. Note the Andreev dI(V)/dV structures are reproducible under R N and m variation, therefore they are caused by bulk properties of LiFeAs. The difference The corresponding dynamic conductance spectra at 4.2 K. Normal-state nonlinear dI(V)/dV background is suppressed. Vertical dashes mark the positions of Andreev features of three possible superconducting order parameters: the large gap Δ Γ ≈ 6.1 meV, the second, anisotropic gap with the edges Δ out L ≈ 3.8 meV and Δ in L ≈ 2.7 meV, and the small gap Δ S ≈ 1.25 meV. The inset shows the low-bias fragments of these dI(V)/dV spectra with additional background suppression in order to reveal some footprints of the second subharmonic of the small gap (arrows, Δ S labels). V norm ≡ V ∕m in the R N of the single contacts (see Fig. 1a) could be attributed to a minor change in the contact area (in the ab-plane) or in the barrier transparency during the readjustment.
Above T c the CVC and the dI(V)/dV of the SnS junction formed in LiFeAs remain nonlinear (see dashed line in Fig. 1a), as mentioned by us earlier in [21]. Since the nonlinearity is not related to the superconducting state, we suppress this monotonic background in the spectra shown in Figs. 1, 2, 4 in order to simplify the Andreev structures to that of the classical SnS junction.
At highest bias voltages eV ≈ ±12.2 meV, the Andreev dips caused by the largest superconducting order parameter are present (red bars, 2Δ Γ label in Fig. 1b). Turning to ARPES data [4,12,13], the largest superconducting gap is the most likely developed in the inner hole barrel of the Fermi surface (located tightly around the Γ point), therefore, hereafter we use the notation Δ Γ . Its second ( n = 2 ) subharmonic expected at eV ≈ ±6.1 meV is not distinguished due to a low amplitude even for the fundamental ( n = 1 ) dip. Possibly, it is caused by a small carrier mean free path or low Cooper pair concentration in the corresponding hole bands constituting the shallow barrel. Another reason is an influence of inelastic scattering usually characterized by the broadening parameter Γ E = ℏ∕(2 E ) (where E is the inelastic relaxation time) that smears the electron density of states (DOS) peaks at the gap edges and extends quasiparticle continuum inside the gap, mixing the real and the imaginary parts of DOS distribution, finally resulting in a smearing of Andreev dips.
At lower bias voltages, the spectra in Fig. 1 show a wide doublet at eV ≈ ±(5.4 − 7.6) meV (magenta and blue bars, 2Δ in L and 2Δ out L labels). Since the increase of their amplitude, as compared to the highest bias dip, and the distinct temperature evolution (see below), we attribute this doublet to a fundamental Andreev feature of another superconducting gap(s), hence the positions of the two dips forming the doublet directly determine the two characteristic energy values Δ out L ≈ 3.8 meV and Δ in L ≈ 2.7 meV. For possible interpretations of the origin of the doublet feature, see the discussion. The second subharmonic of Δ out L located at eV ≈ ±3.8 meV (thin blue bars, Δ out L label) is resolved in the dI(V)/dV spectra, whereas the n = 2 dip of Δ in L expected at about ±2.7 meV is overlapped with the pronounced small gap feature at eV ≈ ±2.5 meV (black bars, 2Δ S label). Being located at the drastic foot, the second Δ S subharmonic at eV ≈ ±1.25 meV is not easily visible but becomes distinguished after additional monotonic background suppression in this low-bias region, as shown in the inset of Fig. 1b (arrow, Δ S label).
Temperature evolution of the dI(V)/dV spectrum (upper curve in Fig. 1b) is shown in Fig. 2. The normal resistance R N of the contact that can be estimated at |eV | ≫ 2Δ is constant until T local c , which demonstrates a mechanical stability of this break-junction. In order to clarify the behaviour of the Andreev structures, the dI(V)/dV curves in Fig. 2 are manually shifted upward, their monotonic background is suppressed by normalizing to the normal state spectrum.
At the lower curve in Fig. 2, the fundamental harmonic of the Δ L gap (doublet) and the small gap are marked by black dashes. With temperature increase, the ZBC peak and Andreev dips become less intensive, the latter also shift toward low bias voltages. At T local c ≈ 17.5 K all the features caused by MAR effect totally vanish, thus indicating a contact area transition to the normal state. Tracking the positions of the Andreev dips in Fig. 2, one could directly determine the temperature behaviour of the superconducting energy parameters: the possible inner and outer edges of anisotropic gap Δ L (T ) and Δ S (T ) dependences are presented in Fig. 3 by magenta, blue, and open black circles, respectively. The anisotropy of the Δ L gap estimated as is about 29% and remains almost constant until T local c (see the lower panel of Fig. 3). Since that, it is possible to speculate on the attributing these both dips to one and the same superconducting order parameter. With it, the small gap tends to zero a bit , as demonstrated in the inset: the ratio between Δ in L (T ) and Δ S (T ) increases twice in the vicinity of the critical temperature. The latter observation indirectly forbids to assign all dI(V)/dV-minima marked by the vertical dashes in Fig. 2 to one and the same bulk order parameter.
Similar data for Andreev twin contact (SnSnS) array formed in another LiFeAs single crystal from the same batch is shown in Figs. 4, 5. The local critical temperature of this array is almost the same, about 17.5 K. For all the spectra, R N (T ) is constant until T c , whereas the curves in Fig. 4 are shifted vertically for clarity, and the monotonic background is suppressed. For the lower curve measured at 4.2 K, we have marked the fundamental Andreev harmonics of the superconducting order parameters Δ Γ (arrows), and Δ L and Δ S (vertical dashes). At high bias voltages about eV ≈ ±12.2 meV a feature caused by Δ Γ ≈ 6.1 meV is present, whereas at higher bias the spectra become flat and contain no features. At eV ≈ ±(5.4 − 7.6) meV, a doublet determining Δ in L ≈ 2.7 meV and Δ out L ≈ 3.8 meV is well distinguished. At lower bias voltages eV ≈ ±(2.0 − 3.2) meV another doublet, attributed to the small gap and possibly caused by Δ S anisotropy, is resolved. The positions of the edges of this doublet correspond to the energy parameters which we denote as Δ in S ≈ 1 meV and Δ out S ≈ 1.6 meV. The doublet structure of the small gap is not fully reproducible in our break-junction experiments, thus, we consider its anisotropy as a possibility, but do not want to speculate on this matter.
The temperature dependences of the energy gap parameters Δ Γ (T ) (red squares), Δ in,out L (T ) (magenta and blue circles), and Δ in,out S (T ) (black open circles) obtained using data from Fig. 4 are shown in Fig. 5. As becomes obvious from the lower panel of Fig. 5, the anisotropy of Δ L (violet circles) and Δ S possible splitting (black open circles) remains almost unchanged with temperature increase toward T c . For the Δ L gap, the anisotropy degree is approximately 34-35% which is close to that estimated for above considered contact (see Fig. 3, lower panel), thus demonstrating the reproducibility of the result for different m number of SnS contacts in an array. If we consider one of the superconducting gaps would be surface Δ * , then any coincidence (after the bias axis normalization by m) become hardly possible, since resulting Δ * ∕m obviously depends on random m.

Discussion
The data shown in Figs. 3, 5 are summarized in Fig. 6. For comparison, we also present our earlier data for SnS-Andreev array with lower critical temperature T local c ≈ 15.6 K, published in [21]. In order to compare the temperature trends of the superconducting gap parameters for the contacts with various local critical temperatures, Fig. 6 shows the normalized order parameter dependences 2Δ i (T )∕k B T local c versus arbitrary temperatures T ∕T local c . It is obvious, whereas all the trends are reproducible (especially for the Δ L ), the observed anisotropy ranges for the largest Γ-gap and the small superconducting gaps vary.
The shapes of the Δ(T ) temperature dependences (see Figs. 3, 5, 6) qualitatively resemble those typical for a moderate interband coupling, although it is hard to estimate it directly due to abundance of free parameters in a case of three gap superconductor. Now we briefly discuss alternative interpretations of the complex Andreev structure of the dI(V)/dV spectra.
More detailed consideration is presented in our previous study [21].
The 2Δ Γ , 2Δ out L and 2Δ in L dI(V)/dV features could not be considered as n = (1, 2, 3) SGS dips for one and the same superconducting condensate due to the relation between their amplitudes. If supposing this, one should expect a lowering of the amplitude of the Andreev dips with n number increase in case of semiballistic SnS contact (inelastic l E < 3d , d is the contact dimension) [38]. Contrary, the set of dI(V)/dV spectra of SnS-Andreev contacts obtained by us in LiFeAs single crystals demonstrates the opposite tendency: the 2Δ Γ feature has the lowest amplitude, whereas 2Δ in L dip is reproducibly the most intensive. Therefore, the dI(V)/dV features labelled as 2Δ Γ and 2Δ out,in L are interpreted as fundamental ( n = 1 ) dips determining three distinct energy parameters.
Despite the low amplitude of its fundamental Andreev harmonic, we attribute Δ Γ to a distinct superconducting order parameter rather than to a fine structure. Generally, the fine structure dip could be caused by a boson resonance (or monochromatic photons/phonons emission) during MAR process [43,44,50], where the energy of the bosonic mode is the difference between the positions of the fine structure dip and the "parent" dip (in our case, more intensive 2Δ in L or 2Δ out L doublet features could be parent). Supposing any bosonic mode even existing in the superconducting state only, the temperature dependence of the boson energy is expected to be rather weaker than Δ(T ) ; in case of the above mentioned monochromatic emission of ℏ quantum, their energy could be considered constant in the range 0 < T < 20 K. However, the experimental data in Fig. 6 show that both, Δ Γ (T ) and the difference (Δ Γ (T ) − Δ in,out L (T )) tends to decrease rapidly towards T local c . Therefore, we attribute the highest-bias dI(V)/dV features as the fundamental n = 1 Andreev harmonic of the distinct, largest superconducting order parameter Δ Γ .
As mentioned above, the low intensity of its dI(V)/dV Andreev structure could result from a large broadening parameter Γ > Δ(0) and/or low partial carrier concentration in the related bands. Actually, recent ARPES probes indicated an unexpectedly strong inelastic scattering could be typical for LiFeAs [51]. Noteworthily, in the spectra of the most qualitative SnS junctions (with Γ ≪ Δ(0) ) a narrow doublet at eV = ±2Δ Γ was resolved by us [21] that could be caused by a minor Δ Γ anisotropy in the momentum space. For the single-dip peculiarity observed in this study, the 2Δ Γ (T )∕k B T c data determined from Fig. 2 (see red squares in Fig. 6) lays in the middle of the narrow doublet shown by gray rhombs (taken from [21]).
The doublet Andreev feature caused by the Δ L gap is reproducibly observed in the dynamic conductance spectra of various SnS break-junctions formed in LiFeAs single crystals, regardless to any random parameters: the area, normal resistance R N , and the number of junctions m in the array. Therefore, this doublet represents the intrinsic (bulk) properties of the related bands in LiFeAs characterized by two gap-edge energy parameters Δ in L and Δ out L . Note their temperature behaviour is quite similar, with the anisotropy degree A(T) keeping almost constant until T c (see lower panels of Figs. 3 and 5). With it, the intensity of the inner dip (at |eV | = 2Δ in L ) is reproducibly stronger than that of the outer dip. Earlier, the doublets for the large gap with similar temperature evolution were typically observed by us in (Ba,K)Fe 2 As 2 and Ba(Fe,Ni) 2 As 2 [40,44].
The energy parameters Δ out L and Δ in L could be either two isotropic superconducting order parameters (developing at different Fermi surface sheets) or the maximum and minimum magnitudes of the anisotropic gap Δ L ( ) related to one and the same superconducting condensate (where is the angle in the k x k y -plane). Supposing the first case, a moderate coupling between the related bands would be realized (note in a two-band approach [52,53,54] the Δ 1 (T )∕Δ 2 (T ) = const within the whole temperature range if the determinant of the coupling constant 2 × 2 matrix is zero, and the geometric mean of intraband coupling potentials √ V 1 V 2 is nearly equal to the interband coupling V 12 ).
Turning to the simple numerical simulation of the dI(V)/ dV Andreev feature caused by the anisotropic gap with an extended s-wave symmetry [20,40,41], and accounting the anisotropy of superconducting properties expected in LiFeAs [4,12,13,15,16,18,19], one may suppose the doublet is caused by anisotropic superconducting gap Δ L . In this case, the Δ L gap would be momentum-dependent having about 30% anisotropy, whereas 2Δ out L and 2Δ in L are the maximum and the minimum absolute values of Cooper pair coupling energies in the k x k y -plane of the momentum space. In order to distinguish between these two cases, further detailed studies of the doublet shape are necessary. Unfortunately, our technique cannot directly judge, which of the cases is being implemented.
It is generally possible to interpret alternatively the lower eV features of the small gap Δ S as subharmonics of Δ L . This seems not the case here due to the following reasoning: (i) the amplitude of 2Δ S feature is much larger than that of the 2Δ L feature, (ii) their positions evolve differently with temperature, while a temperature trend of all features, constituting one SGS, should be the same. Contrary, as shown in the inset of Fig. 3, the ratio Δ in L (T )∕Δ S (T ) increases from 2.2 to 4 in the vicinity of T c , therefore the Δ S features have another temperature trend. Moreover, in general, a substantial lowering of the small gap in the vicinity of T c (as compared to the large gap temperature trend) is typical for a "driven" superconducting order parameter, whereas its existence near T c is induced in the k-space by a stronger superconducting condensate, with larger Δ(0) . Hence, we consider Δ S as a distinct superconducting gap, but its symmetry is still ambiguous.
On the one hand, a doublet Δ S feature was resolved in Fig. 4, showing Δ in S (T )∕Δ out S (T ) ≈ const similarly to the Δ L doublet behaviour. On the other hand, such Δ S doublet is poorly reproducible, whereas the majority of the obtained dI(V)/dV spectra show a pronounced single feature of the small gap, resembling that in Fig. 2 (see also stars in Fig. 7, as well as the "break-junction" data in Fig. 5 in [14]). Accounting that the position of this single Δ S -feature (black circles in Fig. 6) is about an average between the doublet (black squares, different SnS-Andreev contact), one may suppose a common merging of the Δ S doublet due to any scattering processes. Of course, one should not exclude any parasitic effects (an interference, shape resonance, etc.) causing the observed Δ S doublet, therefore, this issue also needs in further detailed studies.
Summarizing the above discussion, we assume that have revealed at least three bulk superconducting order parameters in LiFeAs. The temperature dependences and the corresponding characteristic ratios (see Fig. 6) are reproducible, thus demonstrating the self-consistency of our IMARE spectroscopy data. All the directly determined Δ i (T ) dependences lay below the single-band BCS-like curves (dash-dot lines) which is typical for a moderate interband interaction. The largest superconducting gap has the characteristic ratio 2Δ Γ (0)∕k B T c ≈ 8 and possibly shows a minor anisotropy about 10% (as mentioned in [14,21]). For the "large" gap, we observe a moderate anisotropy about 30% (or two distinct superconducting gap values Δ in,out L differing by ≈ 30% ) with the characteristic ratios 2Δ out L (0)∕k B T c ≈ 5 and 2Δ in L (0)∕k B T c ≈ 3.5 . For the small gap, the average 2Δ S (0)∕k B T c is as small as 1.7 which is less than the weak coupling BCS limit, thus indicating a "driven" superconductivity character in the related bands. Figure 7 shows the characteristic ratios of the superconducting order parameters obtained in our IMARE studies and in the literature. Here, the current data are shown by red stars, whereas our earlier data [14,21] are presented by gray stars. For comparison, we also present the r i ≡ 2Δ i (0)∕k B T c values versus critical temperature T c estimated by STM (triangles) [22][23][24], bulk techniques (specific heat and H c1 measurements, squares) [27][28][29][30][31][32], and surface probes (rhombs) [33,34]. Dashed areas cover the ranges of the characteristic ratios of three superconducting order parameters.
Consider the range of the largest r(T c ) values in Fig. 7. Note this superconducting order parameter ( Δ Γ in our notation) is detected in STM and IMARE probes only. Nonetheless, the STM data show a valuable data scattering: the lowest values r Γ ≈ 8 obtained by STM [22][23][24] agree well with our IMARE data, whereas a bit higher values r Γ ≈ 10 were estimated in some STM probes [25,26].
Summarizing our IMARE data (stars in Fig. 7), at least three distinct superconducting order parameters with possible anisotropy are reproducibly observed in Li 1− FeAs. Within the T c ≈ 15.5-17.5 K, the characteristic ratios remain almost constant.

Conclusion
Using incoherent multiple Andreev reflections effect (IMARE) spectroscopy, we have studied the I(V) and dI(V)/ dV of planar break junctions with almost optimal local critical temperature T local c ≈ 17.5 K formed in LiFeAs single crystals. The complex Andreev structure reproducibly observed in the dI(V)/dV spectra could be described in terms of at least three distinct, possibly anisotropic superconducting order parameters coexisting below T c . The directly determined characteristic energy gap parameters at T ≪ T c are Δ Γ ≈ 6.1 meV (likely having a minor anisotropy of about 10%), the middle gap(s) Δ out L ≈ 3.8 meV and Δ in L ≈ 2.6 − 2.7 meV, and the average small gap Δ S ≈ 1.3 meV (possibly also showing about 35% anisotropy). The directly measured temperature dependences of the superconducting gaps are typical for a substantial Fig. 7 The characteristic ratios 2Δ i (0)∕k B T c versus critical temperature T c obtained in Li 1− FeAs single crystals. The current data obtained using IMARE spectroscopy (Figs. 2, 4) are shown by red stars, our earlier data detailed in [14,21]-by gray stars. The characteristic ratios estimated for the large superconducting gap (blue symbols), middle gap (cyan symbols), and the small gap (black symbols) using STM (triangles) [22][23][24][25][26], bulk techniques (squares) [27][28][29][30][31][32], and surface probes (rhombs) [33,34] are shown for comparison. Dashed areas cover the 2Δ i (0)∕k B T c ranges for clarity, the vertically connected symbols illustrate the gap anisotropy. Dashdot line shows the weak-coupling BCS limit 3. interband coupling, cannot be fitted by the single-band BCS-like trends, and coincide well with our previous results for the LiFeAs with T c ≈ 15.6 K. The directly determined characteristic ratios of all the superconducting energy gap parameters remain almost constant in Li 1− FeAs within the range of critical temperatures T c ≈ 15.5 − 17.5 K.
Acknowledgements Some journals require declarations to be submitted in a standardised format. Please check the Instructions for Authors of the journal to which you are submitting to see if you need to complete this section. If yes, your manuscript must contain the following sections under the heading 'Declarations': Author contributions SK and TK contributed equally to the sample mounting and the IMARE spectroscopy measurements. The data analyzis and interpretation were provided by TK. Single crystal growth and characterization were made by IM, AB, and AS. The first draft of the manuscript was written by SK, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations
Competing interests The authors have no competing interests to declare that are relevant to the content of this article.

Ethics approval
The submitted work is original and not has been submitted to another journal for simultaneous consideration. The manuscript is not published elsewhere in any form or language.

Consent for publication Not applicable
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.