Vortex-lattice melting and paramagnetic depairing in the nematic superconductor FeSe

The full H-T phase diagram in the nematic superconductor FeSe is mapped out using specific-heat and thermal-expansion measurements down to 0.7 K and up to 30 T for both field directions. A clear thermodynamic signal of an underlying vortex-melting transition is found in both datasets and could be followed down to low temperatures. The existence of significant Gaussian thermal superconducting fluctuations is demonstrated by a scaling analysis, which also yields the mean-field upper critical field Hc2(T). For both field orientations, Hc2(T) shows Pauli-limiting behavior. Whereas the temperature dependence of the vortex-melting line is well described by the model of Houghton et al., Phys. Rev. B 40, 6763 (1989) down to the lowest temperatures for H $\perp$ FeSe layers, the vortex-melting line exhibits an unusual behavior for fields parallel to the planes, where the Pauli limitation is much stronger. Here, the vortex-melting anomaly is only observed down to T*= 2-3 K, and then merges with the Hc2(T) line as predicted by Adachi and Ikeda, Phys. Rev. B 68 184510 (2003). Below T*, Hc2(T) also exhibits a slight upturn possibly related to the occurence of a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state.


I. INTRODUCTION
In 1957 Abrikosov 1 predicted that a magnetic field can penetrate a superconductor as an array of vortices, each carrying a magnetic flux quantum Φ 0 = h/2e.This occurs in type-II superconductors in which the normal-superconducting surface energy is negative, i.e. when the Ginzburg-Landau parameter κ = λ/ξ, the ratio of the London penetration depth λ to the coherence length ξ, exceeds the threshold value 1/ √ 2. 2 Vortices repel each other and typically crystallize into a hexagonal lattice.In the presence of weak and randomly distributed disorder, e.g. point defects, [3][4][5] this long-range periodicity is lost and a new (dislocation-free) state of matter, still displaying well defined diffraction peaks -the so-called 'Bragg glass'is formed 6,7 .
In increasing magnetic field, the density of vortices increases until they overlap at the upper critical field H c2 (T) where superconductivity disappears at a second-order phase transition 2,8 .However, thermally induced and/or static disorder can lead to a melting of the vortex solid well below the upper critical field H c2 (T) (for reviews see Refs 9-11).
This possibility has been first considered by Eilenberger 12 , but attempts to observe it in low-dimension geometries, to enhance thermal fluctuations, remained unsuccessful 13,14 .Unequivocal thermodynamic evidence of a genuine vortex-lattice melting transition finally came out soon after the discovery of the cuprate superconductors, in which thermal fluctuations are greatly enhanced due to their high T c , very short coherence length and large anisotropy.In a very limited number of exceptionally high-quality single crystals, vortex melting manifests as a tiny discontinuity in the reversible magnetization [15][16][17][18] , while specific-heat [19][20][21][22][23][24][25][26] and thermal-expansion 27 measurements exhibit a peak superimposed on a step in-dicative of the additional degrees of freedom in the hightemperature vortex-liquid phase.This first-order transition represents the only genuine phase transition for superconductivity in a magnetic field, since H c2 (T ) becomes a broad crossover.Since then calorimetric features related to melting were also reported for conventional low-T c superconductors e.g.Nb 3 Sn, SnMo 6 S 8 and for Fe-based superconductors [28][29][30][31] .However, the fate of the melting transition for T → 0 remains unclear since the low-temperature/high-field region is usually inaccessible, as in cuprates due to the high values of H c2 (0), or because residual disorder disrupts the melting transition 26 .Theoretically, the melting line for T → 0 may i) be suppressed due to quantum fluctuations [32][33][34] , or ii) merge with H c2 (T) or iii) even disappear at finite temperature in strongly Pauli-limited superconductors 35 .
The Zeeman effect represents another mechanism which can affect the high-field superconducting phase transitions, and intense research efforts in low-T c unconventional superconductors including organics, ruthenates and heavy fermion, have focused on the emergence of high-field phases where this effect is dominant (also referred to as strongly Paulilimited superconductors).Prominent examples are κ-(BEDT-TTF) 2 Cu(NCS) 2 36 and CeCoIn 5 37,38 , which exhibit thermodynamic evidence of a modulated phase having Cooper pairs with nonzero total momentum and a spatially non-uniform order parameter ∆(r) 39,40 .While for the former the highfield phase appears to be a physical realization of the original Fulde-Ferrell-Larkin-Ovchinnikov state (FFLO) 41,42 , the modulated phase in CeCoIn 5 is believed to result from a particular coupling between d-wave superconductivity and a field-induced incommensurate spin-density wave (SDW) [43][44][45] .
The recently discovered Fe-based superconductors offer another interesting platform for the study of vortex matter.As anticipated theoretically 46 , thermal fluctuations of intermediate magnitude between cuprates and conventional materials, accompanied by a clear vortex-melting anomaly, were highlighted in the 122 and 1144 families via high-resolution thermodynamic measurements 30,31,47 .In parallel, a first-order superconducting transition detected in the magnetostriction of KFe 2 As 2 stressed the relevance of Zeeman depairing in these materials and raised the possibility of observing a FFLO phase [48][49][50] .
Among the Fe-based materials, FeSe has attracted considerable interest as superconductivity emerges deep inside a nonmagnetic but electronic nematic phase that breaks 4-fold rotational symmetry below T s = 90 K 51,52 .Superconductivity is argued to arise from a spin-nematic pairing driven by orbital-selective spin fluctuations [53][54][55][56][57][58] .Despite its low T c ≈ 9 K, FeSe can be considered as a high-T c superconductor because of its very low carrier density 59,60 and Kasahara et al. 61 have argued that it lies deep inside the Bardeen-Cooper-Schrieffer/Bose-Einstein-Condensate (BCS/BEC) crossover.In this context, the same authors have claimed i) that a fieldinduced phase transition of the Fermi liquid with strong spin imbalance occurs for H c within the superconducting state at a field H * at which the Zeeman energy becomes comparable to F 61,62 and ii) that a genuine FFLO phase is observed for H ⊥ c 63 .Thus, both thermal fluctuations and paramagnetic effects are expected to be large in this high-κ superconductor.The moderate value of H c2 (0) < 30 T offers a unique opportunity to study the H − T phase diagram down to the lowest temperatures in clean single crystals.
In this Article, using thermodynamic probes on highquality FeSe single crystals, we demonstrate the existence of sizable field-induced Gaussian superconducting fluctuations using a scaling approach and provide compelling thermodynamic evidence of the existence of an underlying vortexmelting transition.Our analysis of these data also clearly reveals that Pauli depairing exerts a large influence on the vortex-melting properties in high magnetic fields, in particular for H FeSe layers.Here we find that the vortex-liquid phase disappears for T 3 K, i.e. below which the vortexmelting line merges with the H c2 (T) line.Interestingly, such a merging is predicted by mean-field theory 35 .Here, it occurs near the expected tricritical point from which the FFLO phase could emerge in very clean single crystals 64 .Finally, our results exclude that FeSe lies within the BCS/BEC crossover, and we find no thermodynamic signature of the reported highfield phase for H ⊥ FeSe layers [61][62][63] .
This Article is organized as follows.In Sec.II, the experimental methods (crystal growth and specific-heat and thermalexpansion measurements) are explained in detail.In Sec.III, we present our raw specific-heat and thermal-expansion data, which already provide clear evidence for the existence of both large superconducting fluctuations and a vortex-melting transition.Scaling analysis of our thermodynamic data is presented in Sec.IV and the resulting H − T phase diagram is analyzed thoroughly using existing models of the mixed state.The possible occurence of the FFLO state is discussed and a consistent check of our analysis is provided.Conclusions are given in Sec.V.

II. EXPERIMENTAL METHODS
Stoichiometric single crystals of FeSe were synthesized by chemical vapor transport using a eutectic mixture of KCl and AlCl 3 and characterized using single-crystal x-ray diffraction.Samples 1 and 2, used in this work, were taken from batches 2 and 5 of Ref. 65, respectively.Specific-heat measurements were performed on Sample 1 up to 30 T and down to 0.6 K using a home-made miniature AC calorimeter.It consists of a bare Cernox chip from Lakeshore Cryotronics Inc. split in two parts and suspended to a small copper ring with PtW(7 %) wires.One part is used as an electrical heater while the other part is employed to record the temperature oscillations (a few % of the average sample temperature in the range 1-10 Hz , see Ref. 66 for further details).A precise in-situ calibration and corrections for the thermometer magnetoresistance were accounted for in the data analysis.This setup allows to measure the specific heat of minute samples with an accuracy better than ≈ 5 %, inferred from measurements on 6N copper, with a signal/noise ratio of about 10 4 .High-resolution thermal-expansion measurements were carried out on Sample 2 using a home-built capacitance dilatometer 67 with a typical relative resolution ∆L/L ≈ 10 −8 -10 −10 in fields up to 10 T.

III. EXPERIMENTAL RESULTS
Figs 1(a)-(b) display the temperature dependence of δC e (T, H) = C s (T, H) − C n (T, H), the difference between the superconducting-C s (T, H) and the nematic-state C n (T, H) specific heats, for magnetic fields applied perpendicular (⊥) and parallel ( ) to the FeSe layers, respectively.Here, C n (T, H) = γ n T +B 3 T 3 was determined by fitting the 18 T-data where superconductivity is fully suppressed down to 0.5 K.The inferred values of the Sommerfeld coefficient and the Debye term amount to γ n = 6.5 mJ mol −1 K −2 and B 3 = 0.4 mJ mol −1 K −4 , respectively, in good agreement with previous reports 65,68 .Figs 2(a)-(b) display the corresponding inplane thermal expansion δα e (T, H) = α s (T, H) − α n (T, H) measured upon heating (solid line) after cooling in a magnetic field (dashed line).Here, the normal-state contribution α n (T, H) was determined by fitting the 10 T-data for T ≥ 10 K.
A well-defined discontinuity is observed in the zero-field specific heat at T c = 8.9 K, with a width of about 1 K related to disorder (twin boundaries and/or a very small number of Fe vacancies 69,70 ) indicating the transition from the nematic to the superconducting state.A similar anomaly of comparable width, but with a slightly higher T c = 9.1 K, is found in our thermal-expansion data.We note that the anomaly is very mean-field-like, in the weak coupling limit, i.e. at odds with the λ transition expected in the 3d-XY universality class for an interacting Bose-Einstein condensate 71 , which superfluid 4 He belongs to.This rules out that FeSe actually lies deep within the BCS/BEC crossover for which a cusp-like anomaly is expected.For an exhaustive discussion of the BCS/BEC crossover in a two-band model, we refer to Ref. 72, and in the context of high-temperature superconductivity to Ref. 71 and 73.

A. Large superconducting fluctuations
The strength of thermal fluctuations is usually quantified by the Ginzburg number 9,10 given by where ξ⊥ (0), Hc (0) and Γ = ξ⊥ (0) ξ (0) are the respective Ginzburg-Landau values of the in-plane coherence length, thermodynamic critical field and H c2 anisotropy expressed in SI units.These can all be inferred from our thermodynamic data.Here Hc (0) = 0.21 T and ξ⊥ (0) = 3.5 nm are calculated from the zero-field specific heat and the initial slope of H ⊥ c2 (T ) (see Sec. IV B 1), respectively.Γ ≈ 4.5 is inferred from Fig. 6 where we find that the specific-heat curve for H ⊥ = 1 − 2 T matches that of H = 7 T. We obtain Gi = 5 × 10 −4 for FeSe, which is several orders of magnitude larger than in classical superconductors, e.g.Nb (Gi ≈ 10 −11 ), but slightly lower than in cuprate superconductors (10 −1 < Gi < 10 −3 ) 31 .This large value of Gi strongly suggests that thermal fluctuations cannot be neglected in FeSe.
A telltale sign of thermal fluctuations is a broadening of the superconducting transition in magnetic field 28,74,75 .As shown in Figs 1 and 2(a)-(b), a significant broadening of the superconducting transition for both field directions in both measurements is observed.This becomes particularly evident for H ⊥ ≥ 2 T and H ≥ 7 T, where the broadening clearly exceeds the intrinsic transition width.Evidence for large fluctuations are also quite prominent in the field-sweep measurements displayed in Fig 3(a), e.g. for T < 6.8 K where the transition to the normal state extends over several Teslas for both field orientations.A quantitative analysis of this broadening is obtained by the scaling analysis presented in Sec.IV A. As expected for thermally induced fluctuations, this broadening finally reduces progressively with decreasing temperatures for

B. Evidence for an underlying vortex-melting transition
Beside this large broadening of the transition, a small anomaly is clearly resolved on the upward side of the heat capacity anomaly, particularly for H ⊥ ≥ 2 T and H ≥ 12 T (see shaded area in the insets of Figs 1 and 2(a)).To make this feature more visible and to facilitate comparison with other superconductors, we have subtracted the H ⊥ = 0 T (resp.H ⊥ = 7 T) data from those obtained for H ⊥ ≤ 8 T (resp.H ⊥ > 8 T), as illustrated in Figs 2(c)-(d).The broad remaining discontinuity is very reminiscent of the vortex melting transition T m (H) initially reported by Roulin et al. 23 on twinned YBa 2 Cu 3 O 6.94 single crystals and more recently in Ba 0.5 K 0.5 Fe 2 As 2 30 and RbEuFe 4 As 4 31 .It was interpreted as a second-order melting transition between a vortex glass and a vortex liquid 76 .A clear signature of vortex melting is also visible in the H-sweep measurements where specific heat in excess (shaded areas in Fig. 3) is detected for both field orientations.Interestingly, this melting anomaly, which represents the only genuine phase transition in a hard type II superconductor with strong fluctuations, is clearly less pronounced for H FeSe layers.Indeed it is only observed in the range 7 < H < 21 T (see Figs 2(d) and 3(c) and no melting anomaly could be detected below ≈ 3 K (see Fig. 3(c)) where only a broadened mean-field feature persists which in turn vanishes for T ≤ 0.7 K.
Additional evidence for the existence of an underlying vortex-melting transition is obtained from thermal-expansion measurements, which have proven to be a very sensitive and complimentary probe of the vortex matter 27,30 .As shown in the heating curves (after field-cooling) of Figs 2(a)-(b), extra peaks that grow in magnitude with increasing field, are detected slightly below the broadened superconducting transition, at positions which coincide rather well with T m (H) defined as the mid-point of the broadened discontinuity in specific-heat measurements.Their absence in the cooling curves reveals that they are not electronic in origin, but rather related to irreversible magnetostrictive effects at the melting/irreversibility line due to flux pinning.Similar irreversible peaks were already observed at the melting transition of both YBa 2 Cu 3 O 7−δ 75 and Ba 0.5 K 0.5 Fe 2 As 2 30 crystals with weak flux pinning.Lortz et al. 30,41 argued that they are caused by flux gradients that yield non-equilibrium screening currents which form during cooling and suddenly vanish at the melting temperature upon heating.The applied field exerts a force on these currents which is transferred to the crystal lattice via the pinning centers.These types of anomaly in the cuprates were found to exhibit a behavior comparable to a kinetic glass transition and are most likely related to some glassy vortex phase 41,77 , whose phase line however corresponds rather well to the first-order melting line in fully reversible crystals. )

A. Scaling of the specific heat and thermal expansion
A natural way to study thermal fluctuations is to examine the predicted scaling behavior of the specific heat.The same scaling relations apply to the reversible thermal expansion which is closely related to the specific heat through the Ehrenfest or Pippard relations 78 .
In optimally-doped YBCO, thermal-expansion measurements have demonstrated the existence of 3d-XY fluctuations over a ± 10 K interval around T c 79,80 and their persistence in field up to 11 T without phonon-background subtraction 75 .For comparison, analysis of calorimetric data were inconclusive mainly because the fluctuation component represents at most 5 % of the large phonon contribution [81][82][83][84][85][86][87][88][89][90][91] .In FeSe, the phonon subtraction is straightforward because T c is low and superconductivity is fully suppressed for H ⊥ >16 T.However, an analysis in zero field is impossible since the intrinsic transition width ≈ 1 K clearly exceeds the width of the critical region |T − T c | ≤ GiT c ≈ 4 mK.On the other hand, in this critical XY regime, the field introduces an additional length scale H = φ 0 /πH which reduces the effective dimension-Table I. Superconducting parameters of FeSe.ξ(0) and Γ are the Ginzburg-Landau values of the coherence length and Hc2(T ) anisotropy inferred from the initial slope of Hc2(T) (see Sec. IV B 1) and Fig. 6, respectively.κ is the Ginzburg-Landau parameter and Gi the Ginzburg number derived in Sec.III A. H orb (0), Hp(0) and αM are the orbital and Pauli critical fields and the Maki parameter determined in Sec.IV B.  ality 71,75 leading to a broadening of the transition.According to Jeandupeux et al. 90 , the magnetic field breaks the XY symmetry if the correlation length ξ XY ≈ ξ 0 t −2ν exceeds H (ν XY = 0.669 is the XY critical exponent and t = T /T c − 1).Thus, fluctuations are expected to be observed only for ) we obtain H XY ≈ 2 mT clearly indicating that this regime can be neglected.Thus we restrict our analysis to the 3d Lowest Landau Level (3d-LLL) framework 92 .
In this model, the broadening of the transition is enhanced in field by Gi(H) = Gi 1/3 (H/H c2 T c ) 2/3 92 .It is applicable if the field is high enough to confine the Cooper pairs in their lowest Landau level, i.e for H > H LLL = GiH c2 T c ≈ 10 mT for H ⊥ FeSe layers.In the vicinity of the meanfield transition temperature T c (H), Thouless 98 has shown that δC e (T, H) normalized by δC mf (T, H), i.e. the difference in heat capacities expected from mean-field theory, is a universal function of the single scaling variable . ( It measures the temperature shift with respect to T c (H) normalized by the fluctuation broadening 31 and is a temperature-and field-independent constant.In Ginzburg-Landau theory, δC mf (T, H) is temperature independent.However, similar to Nb 74 , δC mf (T, H) has a sizable temperature variation in the transition region as illustrated in Fig. 1.Therefore, we have normalized our measurements to δC mf (T, H) rather than the mean-field discontinuity since we are only concerned with that part of the temperature dependence ascribed to fluctuations 74 .Here, we have determined C mf (T, H) for each field by fitting the data of Fig. 1 to a second-order polynomials for T << T m (H) and have extrapolated it through the transition region.The same procedure is employed to evaluate α mf (T, H) for the cooling curves of Sample 2.
In Fig. 4, we compare our scaled specific-heat data to the calculations of Li and Rosenstein [93][94][95][96][97] (thick solid line) who successfully derived an analytical expression of the 3d-LLL scaling function for −25 < a T < 8, which includes the expected contributions from vortex melting.This expression was found to describe the broadening of the calorimetric transitions and the melting discontinuity in RbEuFe 4 As 4 31 extremely well.An excellent agreement with the Li-Rosenstein calculation is also achieved in both our specific-heat and thermal-expansion data for a large range of field with the constants r T ⊥ = 60 K −1/3 T 2/3 and r T = 160 K −1/3 T 2/3 which lead to Γ = (r T /r T ⊥ )  3) and the values given in Table I, demonstrating the pertinence of our scaling analysis.We note that the 3d-LLL scaling breaks down for large field values because H c2 in Eq.( 3) is no longer T-independent at high fields because higher-order gradient terms in the Ginzburg-Landau functional, neglected in the 3d-LLL approximation, become important.
A similar scaling approach should also work at very low temperatures for field curves at constant temperature.We find that a similar scaling function can be employed to analyze the mixed-state specific heat shown in Fig. 3 ) the scaling function is now with the constant Here, we estimated C mf (T,H) for each temperature by fitting our data to H 2 away from H m (T) (see dotted lines in Figs 5(c)-(d)), which is characteristic of Pauli-limited superconductors.Our scaled specific-heat data are compared to the Li-Rosenstein [93][94][95][96][97] calculation in Figs 5(a)-(b).For H ⊥ FeSe layers, we find that our scaled data precisely collapse on the theoretical curve obtained with r H⊥ = 23 K 2 3 T − 1 3 calculated using Eq.( 5) and the values given in Table I.This agreement confirms the robustness of our scaling analysis and the existence of Gaussian thermal fluctuations in FeSe.
However, for T < 6 K, the 3d-LLL scaling breaks down for H FeSe layers as illustrated in Fig. 5(b).We show in Sec.IV B 2 that it is related to strong paramagnetic effects, which are not accounted for in the 3d-LLL scaling approach.
Further, we note that the mid-point of our broad melting discontinuity lies around a T ≈ −11 i.e. below the Li-Rosenstein value [93][94][95][96][97] (a T = −9.5)calculated for an ideal vortex lattice.We ascribe this difference to the weak flux pinning observed in our thermal-expansion measurements.
As explained in Ref. 99, the influence of disorder on the locus of the melting line can be quantified by the parameter D/c L with D ≈ (j c /j 0 ) 1/2 (j c and j 0 are the zero-field critical-current and depairing-current densities, respectively) and c L the Lindemann number.Using j c ≈ 3 × 10 4 A cm −2 inferred from Ref. 100, j 0 ≈ 10 7 A cm −2 and c L = 0.2 (see Sec. IV B 1), we obtain D/c L ≈ 0.3 << 1 indicating that the observed melting line in FeSe lies very close to the genuine first-order transition line of the defect-free sample 99 .
Hereabove, we have employed expressions derived for a single-band system whereas FeSe is a two-band superconductor.We believe that it is a fair approximation since the kaveraged energy gaps on the electron and hole bands are found almost equal, i.e. ∆ h (k) k ≈ ∆ e (k) k ≈ 1.3 meV. 53,65n the opposite case, e.g.∆ h >> ∆ e , 3d-LLL scaling breaks down because of the existence of two distinct energy modes 101,102 .

B. H − T phase diagram
The above scaling approach yields a mean-field H c2 (T) and, together with the position of the melting anomaly, allows us to construct (H − T ) phase diagrams which are displayed in measurements and triangles are inferred from thermal expansion.The locus of the mean-field H c2 (T ), derived from our scaling analysis, is represented by the red symbols and corresponds to a T , a H = 0 as indicated by the dotted lines in Figs 4 and 5.In the following, we critically analyze these thermodynamically-derived phase diagrams, which are more representative of the real phase transitions than transportderived phase diagrams 103 , since the zero-resistance criteria of the latter only mark the vortex melting line.

H ⊥ FeSe layers
The large value of the Ginzburg-Landau parameter κ ⊥ ≈ 100 clearly indicates that FeSe is a strong type II superconductor and thus represents a good candidate to study the influence of Pauli-depairing effects on the vortex state.The importance of spin paramagnetism is typically quantified by the Maki parameter 104 where gµ B and H orb (0) = −0.727Hc2 T c are the zero-temperature Pauli and orbital critical fields, 105,106 respectively (g is the gyromagnetic factor).In Fig. 7(a), we show the temperature dependence of H orb (T ) computed within the clean-limit Helfand-Werthamer framework 106 using the measured initial slope H c2⊥ = −3 T K −1 (black dotted lines).We obtain H orb⊥ (0) = 20 T which clearly exceeds the experimental value of H c2⊥ (0) ≈ 15 T by a significant amount, strongly suggesting that Pauli depairing is already significant in FeSe for H ⊥ FeSe layers.To account for this effect, we analyze our data within the Werthamer-Helfand-Hohenberg (WHH) formalism 107,108 including Pauli effects.An excellent fit to our data (red dashed line) is obtained for H p (0) = 26.5 T (left as a free parameter) which leads to a moderate value of the Maki parameter α M ⊥ ≈ 1.1.For completeness, we also plot in Fig. 7 the temperature dependence of H p (T) (magenta curve) 105 .
The melting line H m⊥ (T) (blue symbols) exhibits a characteristic upward curvature near T c , similar to that of YBCO, 26 and then crosses over to a quasi-linear dependence at lower temperature, where it finally merges with H c2 (T) at T ≈ 0. The observation of the vortex-melting transition down to T /T c ≈ 0.1 in FeSe allows us to examine the numerous theoretical models put forward to describe this phenomenon.Here, we compare our results to the semiquantitative approach of Houghton et al. 109 based on the Lindemann criterion.In this approach, the flux-line lattice melts when the mean-square amplitude u 2 th of the fluctuations where  4) for H = H m 99 .Using our WHH calculation for H c2 (T) and Gi from Table I, we solve Eq.( 7) and the resulting curve is depicted in Fig. 7(a) (solid blue line).For c L = 0.15, our calculation fit our data remarkably well down to T = 0. Thus, the melting transition in the presence of weak/moderate paramagnetic depairing remains the only genuine thermodynamic phase transition for H ⊥ FeSe layers, as emphasized by Adachi and Ikeda 35 , and the vortex liquid smoothly crosses over to the nematic state near H c2 (T).Finally, we note that we find no thermodynamic evidence of the high-field modulated phase reported by Kasahara et al. 61 and Watashige et al. 62 inferred from heat-and electrical-transport and torque measurements.

H FeSe layers
The situation is very different for H FeSe layers since paramagnetic effects are much stronger in this direction due to the larger value of α M = Γα M ⊥ ≈ 4.5.In such a case, Adachi and Ikeda 35 predicted that H m (T) and H c2 (T) could already merge at finite temperature.This appears to be the case realized in our data, as illustrated in Fig. 7(b) where the two lines merge around T * ≈ 3 K.Below this temperature, the normal state is recovered before the vortex solid had a chance to melt due to strong paramagnetic pair-breaking, which strongly suppresses H c2 (T).The vortex liquid phase thus no longer exists in this part of the phase diagram, and H c2 (T) again turns into a genuine phase transition for T < T * .This reveals that paramagnetic effects tend to suppress superconducting fluctuations, as can be expected as the transition becomes first-order like.We note that the progressive disappearance of the vortexmelting anomaly nicely correlates with the drastic reduction of the resistive transition width reported by Kasahara et al 63 .
In Fig. 7(b), we show our calculation of H c2 (T) (red dashed curve) using the WHH formula with the inferred values of α M = 4.5 and H orb = ΓH c2⊥ = -13.8T K −1 .We find that it accurately reproduces our H c2 (T) data only for T > T * .For T < T * , the H c2 (T) values lie above this line, and the behavior is reminiscent of FFLO-type physics.The vortex melting line (solid blue curve), obtained within the Lindemann approximation with the same value of Gi (see Sec. IV B 1) and a slightly higher c L = 0.2, describes our experimental data quite well only for T > 6.5 K. Below 6.5 K, the real melting (blue symbols) deviates strongly from this line, which is not unexpected, since strong paramagnetic effects are not accounted for in this model.We note that the 'jump' observed by Kasahara et al 63 in heat-transport measurements appears to coincide with the dashed red line at low T. Theoretically, this line no longer represents a genuine phase transition for T < T * , as explained in Refs 8 and 105, and the correlation with our experiments is unclear.
In purely Pauli-limited superconductors, a spatiallymodulated phase is predicted to appear at high field for T < T cp = 0.56 T c , as shown independently by Fulde and Ferrell FF, ∆(q)e iqr and Larkin and Ovchinnikov (LO, ∆ (r) cos (qr)) 39,40 .Two effects are already expected to mark the emergence of the FFLO phase below T cp : (i) a firstorder transition between the uniform BCS and the modulated states (dotted magenta line) and (ii) an enhancement of H c2 (T) which now defines the transition between the FFLO-and the polarized normal-state phases (solid magenta line).However, these lines are expected to lie very close to each other in purely Pauli-limited 3D systems.Here, they are only 2 T apart at T = 0 (see shaded area in Fig. 7(b)).Accounting for the orbital effect, we expect a shift of T cp towards lower temperature and the two lines to lie even closer to each other [111][112][113][114] .For α M ≈ 4.5, the FFLO state should emerge theoretically below T cp = 0.33 T c ≈ 3 K in FeSe 64 .Interestingly, this corresponds to T * where H m (T) and H c2 (T) are found to merge.Therefore, we have recalculated H c2 (T) (without additional parameters) allowing for a finite modulation of the order parameter, Q = v F q/2T c , assuming a second-order FFLO/N transition 115 .The results for H c2 (T ) and the related Q(T) are displayed in Figs 7(b)-(c) (green solid lines), respectively.They reproduce the enhancement of H c2 (T) observed for T < T * remarkably well.Unfortunately, we found no experimental evidence for the required BCS/FFLO transition in our specificheat measurements.However, recent heat-transport data from Kasahara et al. 63 provide some evidence for such a modulated phase for T < 2 K.
It appears that the transition to the normal state possibly changes its character below T*.While it appears secondorder like at T = 3 K, the transition exhibits a rather broad specific-heat discontinuity at T = 0.7 K (see Fig. 3(d)).This feature is reminiscent of the broadened discontinuity observed around ≈ 27 T in heat-transport measurements 63 .Overall our phase diagram is found to be in rough agreement with that of Kasahara et al. 63 .However, additional measurements (e.g.magnetocaloric measurements) would be very useful to further establish the firm existence of an FFLO phase in FeSe.

Consistency check of our analysis
It is worth noting that the analysis presented in Sec.IV B 1 and IV B 2, for both field orientations, were carried out with only two free parameters i.e.H p (0) = 26-29 T and c L = 0.15-0.20 for determining H c2 (T) and H m (T), respectively.The other quantities reported in Table I are directly inferred from our measurements using standard thermodynamic relations.Similarly, the agreement between the experimental values of r T and r H in Sec.IV A with these calculated directly from Eqs (3) and ( 5) is better than 15 %.
Further, the transition at T = 0 in the purely paramagnetic case occurs when the polarization energy equals the condensation energy, i.e. for 8 : (χ n − χ s )H 2 p (0) = H 2 c (0) where χ n = (µ 0 /2)(gµ B ) 2 N (0) represents the normal-state Pauli susceptibility and χ s = 0 for a singlet superconductor.Using our values of H p (0), H c (0) = 0.12 T from our δC e (T,H=0) data and ∆(0) = 1.3 meV from BQPI experiments, we obtain χ n ≈ 1.7×10 −5 and N (0) ≈ 2.2×10 47 J −1 m −3 spin −1 .These values lead to a value of the Sommerfeld coefficient γ n = 2π 2 3 k 2 B N (0) ≈ 6.7 mJ mol −1 K −2 in excellent agreement with the value inferred from direct specific-heat measurements.These consistency checks provide us with great confidence concerning the relevance of our scaling analysis, the accuracy of the inferred mean-field H c2 (T), and the validity of the presented phase diagrams.

V. CONCLUSIONS
We have determined the full H −T phase diagram of the nematic superconductor FeSe for both field orientations.Compelling evidence of an underlying vortex-melting transition is found in both specific-heat and thermal-expansion measurements down to low temperature and high magnetic fields.We demonstrate the existence of significant Gaussian thermal fluctuations via a scaling analysis of our thermodynamic data which yields the temperature dependence of the mean-field upper critical field.The antagonist interplay between superconducting fluctuations and Pauli depairing effects is studied.We argue that the predominance of the paramagnetic limitation at low temperature is responsible for the unusual disappearance of the melting transition at finite temperature, around T * ≈ 2 K, for H FeSe planes, as anticipated theoretically.A slight upturn of H c2 (T) for T < T * , possibly related to the occurence of the Fulde-Ferrell-Larkin-Ovchinnikov phase, is observed.Additional thermodynamic measurements e.g. of the magnetocaloric effect or magnetostriction, with accurate in-plane field alignment, are necessary to firmly establish the existence of this modulated phase.

Figure 1 .
Figure 1.(Color online) (a)-(b) Temperature dependence of δCe/T, the difference between the superconducting-and normal-state specific heats of Sample 1 for H ⊥ and to the FeSe layers, respectively.The insets highlight the excess specific heat (shaded area) ascribed to the vortex melting transition Tm(H) (see arrows).
l .h e a t .

Figure 2 .
Figure 2. (Color online) (a)-(b) Temperature dependence of δαe, the difference between the superconducting-and normal-state thermal expansions of Sample 2 for H ⊥ and to the FeSe layers, respectively.Measurements were carried out upon heating (solid line) after cooling in a magnetic field (dashed line).The inset highlights the excess thermal expansion (shaded area) related to the vortex-melting transition Tm(H) (see arrows).(c)-(d) Specific heat of Sample 1 as a function of temperature for H ⊥ and to the FeSe layers, respectively.To emphasize the discontinuity at Tm(H), the 0 T-data (7 Tdata) are subtracted from those obtained for H ⊥ ≤ 8 T (for H ⊥ > 8 T).In (b), only the 0 T-data are subtracted.

2 Figure 3 .
Figure 3. (Color online) (a) Field dependence of δCe(H, T )/T , the difference between the mixed-and the normal-state specific heats of Sample 1 for both H and ⊥ FeSe layers at the indicated temperatures.The shaded areas indicate the excess specific heat related to the vortex-melting transition at Hm(T).(b)-(c) same as in (a) but at lower temperatures.No melting anomaly is observed for H layers (c) and only a broadened mean field jump, vanishing below ≈ 1 K, is observed.(d) Comparison of δCe(H, T )/T for the two orientations, plotted as a function of H/Hc2, at T = 0.6-0.7 K.

Figure 4 .
Figure 4. (Color online) (a)-(b) 3d-LLL scaling of the T-dependent specific heat of Sample 1 for H ⊥ and H to the FeSe layers, respectively.(c)-(d) Same for the thermal expansion of Sample 2 (cooling curve).The thick line represents the scaling function calculated by Li and Rosenstein 93-97 for the given values of rT (see text for details).The discontinuity at aT = -9.5 corresponds to the vortex-lattice melting transition.The dashed line indicates the mean-field Tc(H).

3 2 = 4 . 3 .
These values are in very good agreement with the respective values 68, 185 and 4.5 calculated using Eq.( (a).The argument of

Figure 5 .
Figure 5. (Color online) (a)-(b) 3d-LLL scaling of the H-dependent specific heat of Sample 1 for H ⊥ and H to the FeSe layers, respectively.The thick line represents the scaling function calculated by Li and Rosenstein 93-97 for the given values of rH (see text for details).The discontinuity at aH = -9.5 corresponds to the vortexlattice melting transition.The dashed line indicates the mean-field Hc2(T).(c)-(d) Field dependence of δCe(H, T )/T , the difference between the mixed-and the normal-state specific heats of Sample 1 plotted as a function H 2 for both field orientations.

Figs 7 (Figure 6 .
Figure 6.(Color online) Comparison of the specific-heat curves measured for H perpendicular (red) and parallel (black) to the FeSe layers.The inset shows an estimate of the temperature dependence of the Hc2(T ) anisotropy.The line represents the ratio H c2 (T )/H c2⊥ (T ) obtained in Sec.IV B .

Figure 7 .
Figure 7. (Color online) (a)-(b) H − T phase diagram of FeSe for H ⊥ and FeSe layers, respectively.Blue (red) symbols represent the melting (mean-field upper critical field Hc2(T)) data.Black dotted lines are calculations of the orbital critical field H orb (T).Dashed red lines represent fits to the Hc2(T) data including Pauli limitation (assuming a second-order transition).The solid blue lines represent calculation of the melting within the Lindemann approximation (see text).T * marks the temperature where Hm(T) and Hc2(T) merge.The green solid line is the same as the dashed red line allowing for a FFLO modulation Q(T) (shown in (c)), assuming a second-order phase transition between the FFLO and the normal states.For comparison, the solid and dotted magenta line are second-and first-order transitions calculated for the pure Pauli-limited case.The shaded area denotes the extent of the FFLO region in the H − T plane.Tcp stands for the tricritical point.
0) and t = T /T c .Here, c L ≈ 0.1 − 0.2 is the Lindemann number and the analytical expression of F T is given in Ref.110.We note that, for b m ≈ 1, Eq.(7) leads to r H [H c2 (t) − H m (t)]/[H m T ] 2/3 ≈ cste which coincides with Eq.(