Thermodynamics of transition to BCS-BEC crossover superconductivity in FeSe$_{1-x}$S$_x$

The BCS-BEC crossover from strongly overlapping Cooper pairs to non-overlapping composite bosons in the strong coupling limit has been a long-standing issue of interacting many-body fermion systems. Recently, FeSe semimetal with hole and electron bands emerged as a high-$T_{\rm c}$ superconductor located in the BCS-BEC crossover regime, owing to its very small Fermi energies. In FeSe, however, an ordinary BCS-like heat-capacity jump is observed at $T_{\rm c}$, posing a fundamental question on the characteristics of the BCS-BEC crossover. Here we report on high-resolution heat capacity, magnetic torque, and scanning tunneling spectroscopy measurements in FeSe$_{1-x}$S$_x$. Upon entering the tetragonal phase at $x>0.17$, where nematic order is suppressed, $T_{\rm c}$ discontinuously decreases. In this phase, highly non-mean-field behaviors consistent with BEC-like pairing are found in the thermodynamic quantities with giant superconducting fluctuations extending far above $T_{\rm c}$, implying the change of pairing nature. Moreover, the pseudogap formation, which is expected in BCS-BEC crossover of single-band superconductors, is not observed in the tunneling spectra. These results illuminate highly unusual features of the superconducting states in the crossover regime with multiband electronic structure and competing electronic instabilities.


I. INTRODUCTION
The crossover between weak-coupling BCS theory of superconductivity and strong-coupling BEC describes the fundamentals of quantum bound states of particles, including information on the superconducting critical temperature (T c ) and the possible pseudogap formation [1,2]. In the strong-coupling regime, the pairformation temperature T pair and the actual superconducting transition temperature T c are distinctly separated ( Fig. 1(a)). Between T pair and T c , the so-called preformed Cooper pairs exist, which give rise to large superconducting fluctuations and possibly a depletion of the low-energy density of states (DOS), namely the pseudogap. Almost all experimental studies on the crossover physics have been performed in ultracold atomic systems for the past decades [2], because of the difficulty in tuning the attractive interaction between electrons in solids. Therefore, it is extremely challenging to realize such strongly-coupled pairs in electron systems. In underdoped high-T c cuprates, the pseudogap formation has been discussed in terms of the BCS-BEC crossover [3], but recent experiments support a phase transition of different kinds at the onset of pseudogap and the relevance of BCS-BEC crossover in cuprates remains unresolved [4].
Indeed, in high-quality single crystals with long mean free paths of electrons [8], large superconducting fluctuations have been reported from the torque magnetometry, providing support for the presence of fluctuated Cooper pairs created above T c in the crossover regime [11]. These observations underline the importance of the BCS-BEC crossover physics both in the superconducting and normal states of FeSe. On the other hand, scanning tunneling spectroscopy (STS) measurements do not observe the pseudogap formation above T c [12].
These results imply that the superconductivity in FeSe is not simply described by a picture of the conventional BCS-BEC crossover. An important aspect which has not been taken into account in the previous studies is the effects of multiband electronic structure, which may have other tuning factors of the crossover in addition to ∆/ε F . In fact, Fermi surface of FeSe is composed of separated hole and electron pockets with different values of ∆/ε F [8,9,13]. It has been pointed out that in the multiband systems the developments of superconducting fluctuations and pseudogap formation are sensitively affected by the interband coupling strength [14], indicating that the normal and superconducting properties associated with the crossover in multiband systems can be different from those in the single band systems [15]. It is quite important, therefore, to clarify how the superconducting properties evolve when the multiband electronic structure is tuned.
Here, we focus on FeSe 1−x S x , where the orthorhombic (nematic) phase transition at 90 K in FeSe can be tuned  [1,2]. The temperature T is normalized by the Fermi temperature TF and the strength of attraction is given by dimensionless coupling constant 1/(kFas). Below the pairing temperature Tpair, bosonic pairs form, while the superconducting coherence is acquired when the condensation occurs at Tc. In the weak-coupling BCS limit, ∆/εF is much smaller than unity, and the coherence length ξ is much longer than the interparticle distance (∼ 1/kF) and Tpair is very close to the actual Tc. In the BCS-BEC crossover regime ∆/εF ∼ 1, Tc/TF reaches a maximum, and preformed Cooper pairs are expected to exist in an extended temperature region between Tpair and Tc. (b) Experimentally determined T -x phase diagram of FeSe1−xSx. The superconducting transition temperature Tc (red squares) is determined by the present heat capacity measurements using small crystals. The nematic transition temperature (blue diamonds) determined by transport measurements [16] is also plotted. The abrupt change in Tc indicates a significantly different superconducting ground states (SC1 and SC2) divided by the nematic end point at xc ≈ 0.17.
by isovalent sulphur (S) substitution. The nematic transition that distorts the Fermi surface is completely suppressed at a critical concentration x c ≈ 0.17 [16,17]. The superconductivity persists in the tetragonal (nonnematic) phase above x c . In the case of applying hydrostatic pressure, the suppression of nematic order is accompanied by the pressure-induced antiferromagnetic order [18]. In contrast, S substitution suppresses the nematic order without inducing the antiferromagnetism that significantly changes the Fermi surface through band folding even when crossing x c . This enables us to investigate the evolution of superconducting state in the BCS-BEC crossover with little influence of antiferromagnetic fluctuations. Moreover, the observations of quantum oscillations up to x ≈ 0.19 [19] demonstrate that the S ions do not act as strong scattering centers.

II. HEAT CAPACITY
Through the high-resolution heat capacity, magnetic torque, and STS measurements on high-quality single crystals, we investigate the evolution of superconduct-ing and normal state properties with x in FeSe 1−x S x . In Fig. 1(b), we show the phase diagram obtained by the heat capacity measurements. The drastic change in T c and superconducting fluctuations are found across x c . To ensure the sample homogeneity, we used very small crystals (10-190 µg) for the thermodynamic measurements, and employed the long relaxation method with a homemade cell ( Fig. 2(a)) designed for small single crystals. Figures 2(b)-(e) depict the T -dependence of electronic heat capacity divided by T , C e /T , for x = 0, 0.05, 0.10 and 0.13. Here we subtract the phonon contribution, which shows T 3 -dependence, as confirmed by the measurements in the normal state at high fields (Supplementary Materials). In Fig. 2(b), C e /T (T ) in FeSe (x = 0) exhibits an ordinary BCS-like jump at T c = 9 K, consistent with the previous studies on larger samples [20][21][22]. As the temperature is lowered below T c , C e /T decreases with decreasing slope. Below ∼ 3 K, C e /T shows Tlinear dependence, which is consistent with the strongly anisotropic superconducting gap [8,22,23]. Similar C e /T (T ) is observed in orthorhombic phase x =0.05, 0.10, and 0.13 with T c of 9.7, 10 and 9.5 K, respectively (Figs. 2(c), (d) and (e)). Upon entering the tetragonal phase at x > x c , C e /T exhibits a dramatic change not only in the superconducting state but also in the normal state. As represented by x = 0.20 in Fig. 2(h), C e /T increases with increasing slope in the normal state when approaching T c from high temperatures. In addition, C e /T exhibits a kink at T c , without showing a discontinuous jump. Such a continuous T -dependence at T c is reminiscent of the BEC transition in the free Bose gas systems [24] although it is too simple to employ the free Bose gas model to our system. We note that similar enhancement of C e /T in the normal state is obtained in the recent calculations for the strongly interacting Fermi systems near the unitary limit in the BCS-BEC crossover regime [25]. As the temperature is lowered below T c , C e /T (T ) for x = 0.20 decreases with increasing slope, in stark contrast to C e /T (T ) in the orthorhombic crystals as shown in Figs. 2(b)-(e). Below T c , C e /T decreases nearly in proportion to T 1/2 with large residual C e /T in the limit T → 0, which is indicative for a large DOS at zero energy. The large residual DOS is consistent with the recent thermal conductivity and STS measurements [26,27]. Similar C e /T (T ) is also observed for x = 0.18, and 0.21 in the tetragonal phase as shown in Figs. 2(g), and (i). These results together with the sharp change in T c across x c demonstrate that the superconducting properties in the tetragonal phase are drastically different from those in the orthorhombic phase.
We stress that the observed peculiar C e /T (T ) in tetragonal phase does not stem from chemical inhomogeneity within the small sample for the following reasons. First, topographic image in STM of FeSe 1−x S x crystals reveal no evidence for the segregation of S atoms, indicating an excellent homogeneity [27] (see also Figs. 3(a) and (b)). Second, the observation of the quantum oscillations of the crystals in the same batch [19] demonstrates a long mean free path of the carriers, ensuring the high quality of our crystals. Third, according to the elemental mapping of the energy dispersive X-ray spectroscopy (EDX) at different scales of distance (see Fig.S1), the typical spatial variation of the composition ∆x is ∼ 0.01, and there are no discernible segregations or large inhomogeneity of the chemical compositions which can explain our observation in C e /T (T ). Fourth, we measured the heat capacity on two different crystals for x = 0.20(see Fig. S2), and two other concentrations x = 0.18 and 0.21(Figs. 2(g), and (i)), and confirmed the reproducibility of the behavior. We also note that our results are consistent with the previous data of the crystal with much larger volume within the error bar (see Fig. S2). Fifth, most importantly, there is no concentration of the sample, which shows the superconducting transition between ≈ 4 K and ≈ 7 K in the phase diagram of FeSe 1−x S x . This is supported by C e /T (T ) of the crystal at x ≈ 0.16, which is located in the very vicinity of x c . As shown in Fig. 2(f), C e /T at x ≈ 0.16 exhibits two well separated peaks at T c1 = 4 K and T c2 = 7 K, indicating a phase separation (due to the small spatial variation of ∆x ∼ 0.01) at x c . This T -dependence is well reproduced by the sum of C e /T of x = 0.13 (T c = 9.5 K) and 0.18 (T c = 4.0 K), assuming that both contribute equally. Here, for the fitting, T c values are shifted slightly. These results indicate a discontinuous change of T c at x c . Such a discontinuous change in the superconducting properties is consistent with the jump in ∆ values at x c recently reported by STS studies [27]. The observed enhancement of C e /T above T c for the sample in tetragonal phase, therefore, can be attributed to an intrinsic electronic property.

III. SCANNING TUNNELING MICROSCOPY/SPECTROSCOPY
To elucidate whether the depletion of DOS associated with the pseudogap formation occurs, the STS measurements on the crystal of x = 0.25 are performed. By counting the number of S atoms in the STM topographic image, we determined x value accurately ( Fig. 3(a)).  Fig. 3(a). Spectra show little variation within 10 nm scale, demonstrating the uniform spatial distribution of the superconducting gap. We also note that the gap value is consistent with the previous report [27]. Figure 3(c) depicts the STS spectra normalized by the conductance above the superconducting gap in a wide temperature range. Large residual DOS, which is more than half of the normal state value outside the superconducting gap (∆ ≈ 1 meV) is observed even at T = 0.4 K (T /T c = 0.1). This large residual DOS is consistent with the large C e /T at low temperatures. Remarkably, STS spectrum becomes nearly energy independent at high temperatures above ∼ 5 K, showing no signature of the gap formation. This is also seen by the T -dependence of the gap depth plotted in the upper panel of Fig. 3(d). The superconducting gap closes at ≈ 5 K below which the T -dependence of magnetization shows a rapid decrease as shown in the lower panel of Fig. 3

(d).
Nearly energy-independent DOS at ε F observed by STS above T c indicates that the observed enhancement of C e /T above T c in tetragonal phase is not caused by the DOS effect. Therefore, it is natural to consider that the enhancement of C e /T stems from the superconducting fluctuations. Indeed, similar enhancement of C e /T toward T c due to the strong fluctuations has been previously reported in other iron-based superconductors [28][29][30]. In order to discuss the effect of the superconducting fluctuations in FeSe 1−x S x quantitatively, we extract the superconducting contribution C sc , which is obtained by subtracting the normal state electronic contribution γT  from C e for x=0 (blue circles in Fig. 4(a)) and x=0.20 (red circles in Fig. 4(b)). We compare C sc with the conventional term of the mean field Gaussian fluctuations [31](Supplementary Materials). Obviously, the heat capacity contribution that significantly exceeds C gauss can be seen for x = 0.20. This extra heat capacity is observable up to t ∼ 0.8. In Fig. 4(c), we display the x > x c is reduced far below the mean-field BCS value. As a consequence of unusual suppression of ∆C/γT c , the low-T C e /T is largely enhanced through the entropy balance. Figure 4(d) depicts the x-dependence of the C e /T at T /T c = 0.1 taken from Fig. 2. Remarkable enhancement of C e /T at x c indicates the fundamental difference of the low-energy excitations across x c . We note that recent NMR measurements for x > x c report no splitting of the spin-echo signal in the crystal from the same batch [32,33], which excludes the microscopic phase separation of superconducting and non-superconducting re-gions.

IV. MAGNETIC TORQUE
The anomalous C e /T in the tetragonal phase is reflected by the entropies in the normal and superconducting states, S n and S s , respectively, as shown in Fig. 4(e) and (f). For x = 0, the T -dependence of S s shows a kink at T c , which is typical for the second-order superconducting transition. For x = 0.20, on the other hand, owing to the absence of the jump in C e /T at T c , no kink anomaly is observed in the T -dependence of S s . Moreover, the entropy is significantly suppressed from S n even above T c . In Fig. 4(f), we show the x dependence of the lost entropy normalized by S n , (S n − S s )/S n , just above T c , T = 1.03T c . Upon entering the tetragonal phase, an abrupt enhancement of (S n −S s )/S n can be seen, indicating a drastic change of the superconducting fluctuations associated with the preformed Cooper pairs.
To obtain further insight into the superconducting fluctuations in the tetragonal phase, the magnetic torque was measured by using micro-cantilevers. The torque can sensitively detect the diamagnetic response through the anisotropy of magnetic susceptibility [11]. By measuring the field angular variation of torque as a function of polar angle θ from the c axis, the anisotropy ∆χ = χ c − χ ab is determined, where χ c and χ ab are the magnetic susceptibilities along the c axis and in the ab plane, respectively. As shown in the inset of Fig. 4(h), τ (θ) shows an almost perfect sinusoidal curve, τ (θ) ∝ ∆χ sin θ. Although ∆χ(T ) at high fields is almost independent of T and H, it exhibits strong T and H-dependence at low fields as shown in Fig. 4(h). It is well settled that the fluctuation induced diamagnetic susceptibility of most superconductors including multiband systems can be well described by the standard Gaussian-type (Aslamasoz-Larkin, AL) fluctuation susceptibility χ AL (Supplementary Materials). Here we focus on the H-dependent diamagnetic susceptibility at low fields [11]. Then, ∆χ at low fields is much larger than the |χ AL | which is added to high field 7 T data. This implies that the superconducting fluctuations in FeSe 1−x S x are distinctly different from those in conventional superconductors, supporting the fluctuation effects observed in heat capacity measurements. We stress that although the multi-gap superconductivity may lead to a small value of ∆C/γT , it cannot give rise to the non-mean-field behaviors observed in heat capacity and magnetic susceptibility above T c .

V. DISCUSSION AND CONCLUSION
In FeSe, while large superconducting fluctuations that well exceed the Gaussian fluctuations are observed, the C e /T around T c does not exhibit a significant deviation from the mean-field behavior. This suggests that tetragonal FeSe 1−x S x are closer to the BEC regime. However, with increasing S composition, the volume of the Fermi surface increases, as reported by quantum oscillation experiments [19], and ∆ is reduced, according to STM measurements [27]. Then ∆/ε F is expected to become smaller with increasing S concentration, which is opposite to the tendency of approaching the BEC regime [34]. This suggests that there is another important factor, which has not been duly taken into consideration.
A key factor that controls the BCS-BEC crossover in FeSe 1−x S x is the multiband character. As shown previously [12,15], in FeSe the strong interband interaction between electron and hole pockets prevents the splitting of the T pair and the T c , determined by the superfluid stiffness. We note that this strong interband interaction is also an important ingredient in the orbital-selective scenario [23,[35][36][37] and therefore we expect that even taking orbital selectivity into account, the main result, i.e. mean-field like jump of the specific heat will remain the same. At the same time, much less is known on the evolution of the ratio of the interband versus intraband interactions upon approaching x c in FeSe 1−x S x . Furthermore, recent experimental [38] and theoretical [39,40] works pointed out an interesting possibility of the substantial non-local nematic order on the d xy orbital and/or inter-orbital-nematicity between d xy and d yz orbitals in FeSe. The most important consequence of this phenomenon would be an additional Lifshitz transition in FeSe 1−x S x upon S substitution, which occurs close to x c . This in turn implies the appearance of the incipient band near x c (see for example Fig.6 in Ref. [39]). On general grounds, the presence of the incipient band would enhance the tendency towards the BCS-BEC crossover. In addition, given nearly equal mixture of the s-and dwave symmetry components of the superconducting gap in FeSe [41] to express such anisotropic gap, it is natural to assume that they are competing in the tetragonal phase of FeSe 1−x S x . In particular, the ratio of the pairing interactions in the d-wave and s-wave channels |g d ee |/g s eh may, in addition to the incipient band, be another key parameter of the BCS-BEC crossover as follows from our model calculations (see Fig. S4). The pair-formation temperature T pair and condensation temperature T c can be split more with increasing the ratio of intraband and interband interactions |g d ee |/g s eh (see Fig. S4(b)) without changing ε F . This quantitatively explains the appearance of strong coupling superconductivity for x > x c with no apparent enhancement of ∆/ε F . A remarkable and unexpected feature is the absence of pseudogap above T c in the STS spectra, despite the presence of giant superconducting fluctuations in thermodynamic quantities. It should be noted that the relationship between the preformed pairs and the pseudogap formation is still elusive even in the ultra-cold atoms [2]. In addition, the onset temperature of the pseudogap may dif- fer from that of fluctuations [25]. Moreover, it has been pointed out that the pseudogap phenomena in multiband systems are markedly altered from those in the single band systems. For instance, in multiband systems, the pseudogap formation is less prominent by the reduction of its onset temperature [42]. Quite recently, the BEClike dispersion of the Bogoliubov quasiparticles has been reported in the tetragonal FeSe 1−x S x by using ARPES measurements [34]. Although the direct comparison of the ARPES data with the present thermodynamic results is not straightforward, because ARPES measures only the portion of hole Fermi pocket, it provides spectroscopic evidence of the BCS-BEC crossover. The ARPES measurements also report the presence of distinct pseudogap in the hole pocket. However, we point out that this is not an apparent contradiction to the present STM results. It has been suggested that the DOS of electron pocket is much larger than that of hole pockets. As STM measures the DOS integrated over all bands, STM is insensitive to the pseudogap formation in the hole pocket. Therefore, the STM spectra, combined with the ARPES results, suggest that the pseudogap formation is also orbital dependent.
It still, however, cannot fully explain the reduction of the gap magnitudes found experimentally and further mechanisms may be needed such as influence of magnetic fluctuations which would further reduce the superfluid stiffness [43], a Lifishitz transition near the orthorhombic to tetragonal transition [39], or nematic fluctuations. In FeSe 1−x S x , the nematic fluctuations are peaked at q nem ≈ 0 because its ordered state breaks the rotational symmetry while preserving translational symmetry. Enhanced nematic fluctuations [16] play an important role for the normal state properties [17,44], particularly at x ∼ x c [17]. It has also been pointed out theoretically that nematic fluctuations strongly influence the superconductivity [45,46]. Therefore, it is tempting to consider that the nematic fluctuations and incipient bands in FeSe 1−x S x upon approaching a tetragonal phase can enhance the pairing interaction, leading to the system to approach the BEC regime. Further theoretical and experimental studies are required to uncover whether the nematic fluctuations and multi-band nature affect the crossover physics.
Our thermodynamic studies have revealed a highly non-mean-field behavior of superconducting transition in FeSe 1−x S x when the nematic order of FeSe is completely suppressed by S substitution. This is consistent with the recent ARPES measurements which report the BEC-like dispersion of the Bogoliubov quasiparticles in the tetragonal regime. Our findings demonstrate that FeSe 1−x S x offers a unique playground to study the superconducting properties in the BCS-BEC crossover regime of multiband systems.

APPENDIX A: MATERIALS AND METHODS
Crystal growth and characterization The single crystals of FeSe 1−x S x were grown by the chemical vapor transport technique as described in Refs. [6,16]. The composition x is determined by the energy dispersive Xray spectroscopy.
Heat capacity The heat capacity of the crystals is measured using long relaxation method [47,48], whose experimental setup is illustrated in the inset of Fig. 1(a). A single bare chip of Cernox resistor is used as the thermometer, heater and sample stage, which is suspended from the cold stage by silver-coated glass fibers in order that the bare chip has weak thermal link to the cold stage, and electrical connection for the sensor reading. The mass of the samples measured in this study is 81, 41,171,137,179,50,107,192, and 10 µg for x = 0, 0.05, 0.1, 0.13, 0.16, 0.18, 0.20#1, 0.20#2, and 0.21, respectively. The samples are mounted on the bare chip using Apiezon N grease. The heat capacity of the crystals are typically obtained by subtracting the heat capacity of bare chip and grease from the raw data.
Torque The torque measurements are performed using micro-cantilevers. The small single crystals with typical size of 150 µm×100 µm×10 µm are mounted on the tip of micro cantilevers with small amount of Apiezon N grease. The magnetic field is applied inside the ac(bc) plane by vector magnet, and the field-angle dependence of torque is obtained by rotating the whole cryostat within ac(bc) plane. Magnetic torque τ can be expressed as τ = µ 0 V M × H, where µ 0 is vacuum permeability, V is sample volume, M is magnetization, and H is external magnetic field. Here the H is rotated within the ac (bc) plane, and θ is polar angle from the c axis. In this configuration, the τ is expressed as τ 2θ (T, H, θ) = 1 2 µ 0 H 2 V ∆χsin(θ). Here ∆χ = χ c -χ ab is anisotropy of magnetic susceptibility between ac(ab) plane and c axis.
Scanning tunneling microscopy/spectroscopy The scanning tunneling microscopy/spectroscopy measurements have been performed with a commercial low temperature STM system. Fresh and atomically flat surfaces are obtained by in-situ cleavage at liquid nitrogen temperature in ultra high vacuum. All conductance spectra and topographic images are obtained by a PtIr tip. Conductance spectra are measured by a lock-in technique with a modulation frequency 997 Hz and modulation voltage 200 µV. The backgrounds of raw conductance data are normalized by curves which are derived by 2nd order polynomial fitting using outside superconducting gaps.
Theoretical calculation The theoretical calculations were done within a simplified interacting two-band model in two dimensions, employed previously [15] with hole and electron pockets with small Fermi energies [ Fig. 8]. We assume superconductivity due to repulsive interactions. While the dominant interband interaction favors Cooper-pairing in the s-wave symmetry channel, the dwave projected intraband interaction yields attraction in the d-wave channel. We obtained T pair from the solution of the linearized mean-field gap equations, whereas T c in each case is estimated from superfluid stiffness. Details of the calculations are given in the Supplemental Information.

APPENDIX B: ENERGY DISPERSIVE X-RAY SPECTROSCOPY (EDX) ANALYSIS
Energy Dispersive X-ray spectroscopy (EDX) analysis is conducted on tetragonal x = 0.20 samples in comparison with x = 0 to investigate the spatial homogeneity of the sulfur content. For x = 0, sulfur intensity is almost same as the background outside the sample position in Fig. 5(a), showing that there is no sulfur substituted into the sample as expected. Figures 5(b)-(f) show the elemental distribution of x = 0.20 at the scale of 500 µm, 50 µm, 5 µm, 500 nm, and 50 nm, respectively. Each data shows S intensity with the much thicker color than x = 0 although there is some small darker areas in Fig. 5(b), and (c) which originate from the large impurity attached on the crystal surface and surface roughness in a macroscopic scale and are not from the inhomogeneity of chemical composition. The data at every length scale show the uniform distribution of sulfur which is almost comparable to iron and selenium, indicating that there is no discernible segregation and inhomogeneity down to the mesoscopic scale 10 nm. We quantitatively investigate the sulfur content x at the several area inside the sample and the typical variation of x is 0.01, demonstrating that such variation of composition cannot explain the results of specific heat and torque measurements.

APPENDIX C: COMPARISON WITH PREVIOUS HEAT CAPACITY DATA
The heat capacity in the tetragonal FeSe 1−x S x has been reported previously by Y. Sato et al., using quasiadiabatic method [26]. We plot our C e /T data for two x = 0.20 samples in Figs. 6(a), (b) together with the data of Ref. [26]. The two samples in this work exhibit almost identical temperature dependence and absolute values, indicating the reproducibility of the heat capacity data in x = 0.20 with better resolution in our studies. Although the absolute value of our data is slightly smaller than the data in Ref. [26], our data is within the error range of previous data, and both data are basically consistent up to 6 K which is the highest temperature in the previous reports. The heat capacity jump due to the superconducting transition becomes broader with increasing H for both S concentrations. This behavior is possibly due to the enhanced fluctuations in the magnetic field as discussed in YBa 2 Cu 3 O 7−δ [49]. In zero field, we observe the large superconducting fluctuations in heat capacity in x = 0.20. The additional heat capacity due to the superconducting fluctuations just above T c is suppressed by applying magnetic field, and then the C e /T approaches to the normal state value. This fact indicates that the additional heat capacity responds sensitively to the magnetic field, and does not stem from the incorrect estimation of electronic heat capacity in the normal state when we subtract the lattice contribution. From the T c defined as the peak temperature in C e /T under field, we obtain the H-T phase diagram of x = 0.0, 0.10 and 0.20 as shown in Fig. 7(c). In x = 0.0 and 0.10 the T -dependence of upper critical field H c2 is consistent with previous reports [50,51], while the H c2 for µ 0 H > 12 T cannot be determined from our measurements. The H c2 in x = 0.20 shows small value compared to x = 0.0, reflecting the lower T c . From the H c2 (0) estimated by the WHH relation H c2 (0) = 0.69T c |dH c2 /dT | T =Tc [52], we obtain the coherence length ξ ab = 13.5 nm, ξ c = 4.1 nm for x = 0.20 through H c2 (0) = Φ 0 /(2πξ 2 ab ), and Φ 0 /(2πξ ab ξ c ) for H//c, and H//ab, respectively.
The contribution of the mean field Gaussian fluctuations to the heat capacity [31] is given by C gauss = C + t −0.5 , where C + = k B /(8πξ 2 ab ξ c ) and t ≡ T −Tc Tc is the reduced temperature, and ξ ab and ξ c are in-plane and out-of-plane coherence lengths at T = 0, respectively. The dashed lines in Figs. 4(a) and (b) represent the contribution of Gaussian fluctuations obtained by using ξ ab = 5.5 nm, ξ c = 1.5 nm for x = 0 [11] and ξ ab = 13.5 nm, ξ c = 4.1 nm for x = 0.20 (Fig. 7).
The Gaussian-type (Aslamasoz-Larkin, AL) fluctuation contribution in susceptibility is given by in the zero-field limit. Here Φ 0 is the flux quantum. In the multiband case, the behavior of χ AL is determined by the shortest coherence length of the main band, which governs the orbital upper critical field. As the diamagnetic contribution χ AL is expected to become smaller in magnitude at higher magnetic fields, |χ AL | yields an upper bound for the conventional Gaussian-type amplitude fluctuations.

APPENDIX E: CALCULATIONS OF PAIRING INSTABILITY AND PAIRING CHANNEL BASED ON THE TWO BAND MODEL
To make an estimate of the splitting between the Cooper-pair formation temperature, T pair , and the actual superconducting transition temperature, T c , we followed Ref. [15] and consider a simplified interacting twoband model in two dimensions with hole and electron pockets with small Fermi energies (E h F = 20meV and E e F = 10meV) as depicted in the Fig. 8(a). As we are interested in the evolution of the FeSe system as a function of the Sulfur doping we assumed the system is orthorhombic. However, once the intraband interaction in the d-wave interaction dominates and the s−wave component of the gap is small the results can be easily extrapolated to the purely tetragonal system. The Hamiltonian readsĤ where α ∈ {e, h}, and ξ ke , ξ kh are the electron and hole energy dispersions separated by the large momentum, respectively as shown in the Fig. 8(a). Assuming superconductivity due to repulsive interaction in the A 1g and B 2g symmetry channels we write the interaction terms as follows Here, we further assume that the inter-band repulsion drives s ± -wave symmetry of the superconducting order parameter, while interaction  in the d-wave channel is mostly intra-band. We define superconducting order parameters as and perform a mean-field decoupling to find the meanfield gap equations The dimensionless coupling constants are now given by dimensionless g αα = N 0 U αα with the density of states in two dimensions N 0 = m 2π , Λ = 1eV E h F is the high energy cut-off and E α (k) = ξ 2 α + [∆ s α + ∆ d α cos(ϕ)] 2 are the energy dispersion of the Bogoliubov quasiparticles. Note that while the inter-band term between electron and hole pockets is assumed to be repulsive, the intraband g d ee < 0 and g d hh < 0 are attractive in the d-wave channels. In case of ∆ α /E Fα ∼ 1, where E Fα is the Fermi energy of each band, we need to renormalize the chemical potential assuming the total number of particles is conserved. The equation determining the particle number is given by and it has to be solved self-consistently with Eqs. ( In the usual BCS case the phase fluctuations are costly and T c ≈ T pair . In our case due to smallness of the Fermi energies, the condensation of pairs may happen at lower temperature, T c ≤ T pair . We estimate T c ≈ π 2 ρ s where ρ s (T = 0) is the superfluid-stiffness following Ref. [15]. Note that in two dimensions, the superconducting transition temperature T c ∼ ρ s (T c ) (see, e.g., Refs. [53,54]). The interplay between T c and T pair depends on the ratio ρ s (T = 0)/T pair . If this ratio is large, the superfluid stiffness rapidly increases below T pair . In this situation, T c = T pair minus a small correction, i.e., the phases of bound pairs order almost immediately after the pairs develop (phase fluctuations cost too much energy). If ρ s (T = 0)/T pair is small, ρ s (T ) increases slowly below T pair and T c is of order ρ s (T = 0). In Ref. [15] its expression was extended to the multiband case where ρ s is given by ρ s ≈ ρ e + ρ h , We solved Eqs. (3) together with Eq. (3) for the mixed s + d-wave gaps in the orthorhombic state at T = 0 and used them as an input parameter to calculate T c . The Cooper-pair formation temperature, T pair is found from the condition that the determinant in Eq. (4) vanishes and the results are shown in the Fig. 8 The pair-formation temperature T pair and condensation temperature T c are calculated as a function of the ratio of intraband and interband interactions |g d ee |/g s eh shown in Fig. 8(b). The two temperatures split and the difference grows with increasing |g d ee |/g s eh . Fig. 8(c) shows the evolution of the gap symmetry as a function of |g d ee |/g s eh also in units of E F h , and we can see the dominant d-wave character for each band with large |g d ee |/g s eh regime. Observe that for smaller ratio of |g d ee |/g s eh the pairing symmetry is a pure s +− driven by inter-band interactions. In this case both gaps are equal in magnitude, but with the larger gap at the band with smaller E F . With increasing intra-band coupling in the d-wave gaps ∆ d e and ∆ d h grow while the s-wave gaps become smaller. Since we assume |g d ee | g d eh , ∆ d e and ∆ d h differ in magnitude with the larger gap at the band with larger E F .