Inhomogeneous Superconductivity Onset in FeSe Studied by Transport Properties

Heterogeneous superconductivity onset is a common phenomenon in high-Tc superconductors of both the cuprate and iron-based families. It is manifested by a fairly wide transition from the metallic to zero-resistance states. Usually, in these strongly anisotropic materials, superconductivity (SC) first appears as isolated domains. This leads to anisotropic excess conductivity above Tc, and the transport measurements provide valuable information about the SC domain structure deep within the sample. In bulk samples, this anisotropic SC onset gives an approximate average shape of SC grains, while in thin samples, it also indicates the average size of SC grains. In this work, both interlayer and intralayer resistivity were measured as a function of temperature in FeSe samples of various thicknesses. To measure the interlayer resistivity, FeSe mesa structures oriented across the layers were fabricated using FIB. As the sample thickness decreases, a significant increase in superconducting transition temperature Tc is observed: Tc raises from 8 K in bulk material to 12 K in microbridges of thickness ∼40 nm. We applied analytical and numerical calculations to analyze these and earlier data and find the aspect ratio and size of the SC domains in FeSe consistent with our resistivity and diamagnetic response measurements. We propose a simple and fairly accurate method for estimating the aspect ratio of SC domains from Tc anisotropy in samples of various small thicknesses. The relationship between nematic and superconducting domains in FeSe is discussed. We also generalize the analytical formulas for conductivity in heterogeneous anisotropic superconductors to the case of elongated SC domains of two perpendicular orientations with equal volume fractions, corresponding to the nematic domain structure in various Fe-based superconductors.


Introduction
The superconducting transition often occurs in a non-uniform way, when superconductivity (SC) initially appears in the form of isolated domains, which then acquire phase coherence leading to zero resistance. Such a heterogeneous SC onset attracts a great research activity and takes place in most high-temperature superconductors, including copper-oxide and iron-based families [1][2][3][4], where the spatial inhomogeneity of the SC energy gap has been directly observed in numerous scanning tunneling microscopy (STM) and scanning tunneling spectroscopy (STS) experiments [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. The heterogeneous SC onset is also common to various organic superconductors [15][16][17][18][19][20][21], homogeneously disordered conventional superconductors [22], polycrystals [23], and many other materials. alyze the experimental data on diamagnetic response, one needs to know the approximate size d of SC domains, at least if it is smaller than the penetration depth λ of the magnetic field into the superconductor. In the analysis of experimental data in Refs. [32,33], it was assumed that d λ, but this condition is violated if the SC domains are not larger than the nematic domains. The typical width of nematic domains is d n ∼ 100 nm, and their length exceeds λ, as observed by STM [14]. While the nematic domains in FeSe have been directly measured by STM/STS [11,13,14] and nano-ARPES [27], there is no any evidence that the SC and nematic domains in FeSe coincide. The size of SC domains in the conducting x-y plane can be measured by STS or point-contact spectroscopy, but as far as we know, no such experiments have been performed in FeSe, except for Ref. [12], where the spatial resolution was insufficient to study the SC domain shape and size.
The SC domain size along the z axis can be estimated from the SC percolation in finite-size samples, which corresponds to the onset of zero resistance in thin FeSe samples. A similar mechanism has recently been proposed to explain the anisotropic occurrence of zero resistance in organic superconductors [35]. On the other hand, the SC domain shape can be estimated from the anisotropy of excess conductivity above T c in combination with diamagnetic response data, similarly to Refs. [32][33][34]. Thus, the combination of these data helps to estimate the SC domain size in all directions and to answer the question of whether the SC and nematic domains coincide. This is important for understanding the mechanism of superconductivity in FeSe.
In Section 2, we describe the methods of both the experiment and the theoretical analysis. In Section 3, we generalize the Maxwell-Garnett approximation for elongated SC islands with two perpendicular orientations in anisotropic media. In Section 4, we present the results of our measurements of electronic transport in finite-size FeSe samples, used to estimate the interlayer size of SC domains, and the results of our numerical calculations of the percolation probability as a function of SC volume fraction φ and of sample thickness L z for the preliminary parameters of the domain aspect ratio. In Section 5, we desccribe a detailed theoretical analysis of the obtained and previous experimental data in order to extract useful information about the SC domain structure in FeSe. In particular, we numerically calculated the SC percolation threshold for the sample of relevant finite size and shape and compared it with our experimental data on resistivity to estimate the interlayer size of SC domains. In order to study the shape of SC domains, we also reanalyzed previously-obtained combined experimental data on excess conductivity and on diamagnetic response in FeSe above T c under the assumption that the SC domain width is d ∼ 100 nm < λ. In Section 6, we present the main conclusions.

Experimental
Our experimental method was similar to that described in Ref. [32]. We used high quality platelike single crystals (flakes) of FeSe, grown in evacuated quartz ampoules at a permanent temperature gradient using AlCl 3 /KCl flux technique, as in Ref. [38]. The FeSe mesa structures (microbridges), as shown in Figure 1, were made using the focused ion beam (FIB) technology described in Ref. [39,40] from selected single-crystal samples of thickness of 2-4 µm (see Figure 1 and Figure 2a,b in Ref. [32]). Prior to FIB processing, a gold film was deposited by laser ablation to prepare the electrical contacts to the crystal. The electric resistance was measured in the conventional 4-probe configuration. In order to improve heat exchange, most of the structures were coated with collodium.
It is known that FIB may cause damage to the samples. The typical thickness of amorphous layer damaged by Ga ions in FeSe is about 50 nm. The minimum cross section of our mesa structures is 500 × 500 nm, the size presented in the article is 2 × 2 µm which is much larger than the expected depth of the damaged layer. We also evaluated the resistivity of all our structures and did not notice any strong change in the transport properties of thinner samples caused by FIB exposure. The obtained thin microbridges may crack during cooling. Such "defect" samples are clearly visible under an FIB or an SEM [40]. Additionally, the damage to the sample during its cooling or measurement is easily detected by a sharp jump in resistance and by the analysis of its transport properties. Such "defect" samples are excluded. Since FeSe mesas are rather fragile, we cooled them at a slow uniform rate of about 2 K/min. The cooling rate in the uncracked FeSe samples does not affect their transport properties. Figure 1. Photos of the microbridges, used in our experiment, from two different angles. (a) illustrates the 3D shape of our mesa structure, while (b) shows its dimensions.

The Origin of Anisotropic Resistivity Drop above T c and the Maxwell-Garnett Approximation
A new method, based on combined transport and diamagnetic-response measurements in bulk material, was recently suggested and applied to study the SC heterogeneity in several strongly anisotropic materials, including FeSe [32,33], YBa 2 Cu 4 O 8 [34], and several organic superconductors [34][35][36]. These materials have a layered crystal structure and, hence, a strong anisotropy of electronic properties, which is typical of all ambientpressure high-T c superconductors. The resistivity drop above T c in all these compounds is anisotropic and strongest along the least conducting axis [16][17][18][19]32,33,[41][42][43], which contradicts the standard theory [44,45] of superconducting fluctuations in homogeneous superconductors. This anisotropic effect of nascent superconductivity has been explained and analytically described [32][33][34] using a classical effective-medium model, namely, the wellknown Maxwell-Garnett approximation (MGA) [46], generalized for strongly anisotropic heterogeneous metal with elliptical superconducting inclusions [32][33][34]. This simple model indeed predicts that an incipient superconductivity in the form of isolated domains in anisotropic conductors reduces the electrical resistivity anisotropically with a maximal effect along the least-conducting direction [32][33][34].
The qualitative picture behind this model [32][33][34] is very simple. In strongly anisotropic conductors with interlayer conductivity σ zz σ yy σ xx , the direct interlayer current perpendicular to the conducting layers is small, and the ratio is expressed by the parameter η ≡ σ zz /σ xx 1. However, if SC emerges in a form of isolated domains, there is a second way of interlayer electric current via superconducting islands. If there are few of these domains, the major part of the current path goes in the normal phase. Instead of going along the weakly-conducting z-axis in the non-SC phase, this second path between the SC domains goes along the highly conducting layers, until it reaches another superconducting island, providing next lift across the layers (see Figure 2 for illustration). Then, there is no local current density along the poorly-conducting z direction in the non-SC phase. Hence, the contribution from this channel to interlayer conductivity does not contain the small anisotropy factor σ zz /σ xx 1. Instead, this channel gives another small factor-the volume fraction φ of superconducting phase. The second way makes the main contribution to interlayer conductivity if φ/η 1. For the case of in-plane isotropy, σ yy = σ xx , and co-directional isolated spheroidal SC islands of volume fraction φ 1, the analytical formulas for conductivity is rather simple: where γ = a z /a x is the aspect ratio of main axes of spheroidal SC domains, and σ xx0 is the in-plane conductivity in the absence of SC domains. Note that in Ref. [33] γ denoted the square of this aspect ratio. Equation (1) confirms the above qualitative picture: the interlayer conductivity σ zz indeed consists of two terms. The first term (1 − φ) −1 coincides with that in σ xx , while the second term contains a factor γ 2 φ/η and at γ 2 /η > 1 determines the excess conductivity due to SC domains. Note that the domain size d does not enter Equation (1), which is valid for arbitrary distribution of d, provided their aspect ratio a z /a x ≡ γ remains fixed. Equation (1) was generalized for fully anisotropic case a z = a y = a x and σ xx0 = σ yy0 = σ zz0 in Ref. [34]. The averaged aspect ratios a z : a y : a x of SC islands may be extracted by comparing the measured conductivity with Equation (1), provided the SC volume fraction φ is known, for example, from diamagnetic response [32][33][34]. Alternatively, if the averaged aspect ratios a z : a y : a x of SC islands is known from the anisotropic diamagnetic response, the transport measurements can be used to extract the SC volume fraction φ as a function of some driving parameter, such as temperature, pressure, doping level, cooling rate, etc.

Numerical Calculations of Percolation Threshold
Another method [35], based on the transport measurements in finite-size sample, can be used to extract the averaged domain size. It was initially proposed to explain the anisotropic zero-resistance onset observed in organic superconductors, where zero resistance or a sharp resistance drop several times corresponds to a current percolation along the SC domains. Using this approach, we can calculate the SC volume fraction required for this current percolation for a given shape and size of the sample and of SC domains. However, for this method, the sample dimensions should be comparable or only several time greater than the SC domain dimensions. Otherwise, at the limit of an infinitely large sample, the percolation threshold will be isotropic, and we will not obtain information about the geometry of the domains.
The percolation probability p(φ) was calculated numerically using the Monte-Carlo algorithm. At each step, a random state with the proper number of spheroidal domains with a fixed size d and a fixed aspect ratio γ = a z /a x was generated in a box of dimensions L x × L y × L z , matching to our experiment. The required number of SC domains was determined by the volume fraction φ of the SC phase, and was selected in advance before the main simulation cycle. Each state corresponds to a graph whose vertices are SC domains. The vertices of the graph are connected by edges if the corresponding domains intersect. Thus, the problem of checking the presence of percolation along the axis is algorithmically reduced to finding the connected components of the graph containing vertices corresponding to the SC domains on the opposite edges of the sample, i.e. to the search for a percolation cluster. For each state, corresponding to one realization of SC islands, the percolation along each axis, i.e., the existence of a continuous path via intersecting SC domains, was checked, and the averaging over random realizations was made. About 10 4 -10 5 steps are usually sufficient to estimate the average percolation probability with an acceptable accuracy.

Generalization of Maxwell-Garnett Approximation
To consider elongated SC domains aligned in two perpendicular directions with equal probabilities, we have to generalize the MGA described above in Section 2.2. We start from the general equation for the effective conductivity tensorσ e of a heterogeneous media with M − 1 types of unidirectionally aligned isotropic similar ellipsoidal inclusions inside an isotropic media with conductivityσ 1 (see Section 18 of Ref. [46]): whereσ j = Iσ j is the effective conductivity tensor of the inclusions of type j, I is the 3 × 3 unity matrix, the so-called electric field concentration tensor: andÃ (j) is the diagonal depolarization tensor. For an ellipsoid with main semiaxes a i the depolarization tensorÃ (j) has only the diagonal components A (j) i expressed via the integral (see Equation (17.25) of Ref. [46]): where i = 1, 2, 3 corresponds to x, y, z axes. The depolarization tensor has the property that its trace is unity, i.e., ∑ i A i = 1. The Equation (4) can be expressed via the elliptic integrals, as given by Equations (B1)-(B3) of Ref. [34]. In our case of elongated SC ellipses equally distributed along two perpendicular in-plane axes x and y, we have inclusions of two types j = 2, 3 with equal conductivity, σ 2 = σ 3 = ∞, and with equal volume fractions, φ 2 = φ 3 = φ/2, but with different depolarization tensorsÃ (2) =Ã (3) , because the elongated SC domains are differently aligned. Evidently, A A similar result appears if one takes the SC ellipsoids randomly oriented in the x-y plane. Due to the layered crystal structure, in the normal-metal phase FeSe is strongly anisotropic: the conductivity ratio η ≡ σ zz /σ xx ≈ 0.0025. To describe such compounds with highly anisotropic conductivityσ 1 with diagonal components σ m ii , following the method used in Refs. [32][33][34], we applied the coordinate mapping: where with the simultaneous change of conductivity toσ 1 = σ m I = σ m xx I in Equations (2) and (3).
This mapping does not change the electrostatic continuity equation for the electric potential distribution inside matrix phase 1 of the heterogeneous medium: Hence, the voltage distribution in the original and mapped spaces are given by the same function: V(r). As a result of this mapping, the main semiaxes of SC inclusions change according to the rule: and the tensorsR (j1) andÃ (j) change toR (j1) * andÃ (j) * expressed by Equations (3) and (4) with the replacement in Equation (9). If initially the SC domains are not spherical but have ellipsoidal shape with the principal semiaxes a = a 1 , b = βa 1 and c = γa 1 , then after the mapping to an isotropic media these domains keep an ellipsoidal shape but change the principal semiaxes to: In our case of FeSe µ = 1, because σ m yy = σ m xx , and 1/ √ η ≈ 20 1. Hence, after the mapping, we assumed a z * a x * , a y * , but a x * = a y * for elongated SC domains of the shape resembling that of nematic domains in FeSe. Then one may use the simplified formulas (B4)-(B6) of Ref. [34]: Substituting Equations (11)-(13) for Equation (5), applying the mapping given by Equation (10) and using we obtain: and where e = 2.71828. From Equations (15) and (16) we see that the relative excess conductivity is anisotropic, which can be used to determine the aspect ratio γ = a z /a x from transport measurement, provided another aspect ratio β = a y /a x of SC domains is known.
Equations (15) and (16) are obtained for elongated SC islands, equally distributed along two main in-plane directions. This result is similar to the case of randomly oriented elongated SC islands, which is described by taking the trace of the matrix R (j1) [46]. Evidently, Equations (15)- (17) are invariant under the in-plane coordinate permutation x ←→ y, which changes β → 1/β and γ → γ/β 2 . Below we take β = a y /a x ≥ 1, corresponding to SC domains oriented along y.
Let us compare Equations (15) and (16) with Equation (1) derived for spheroid SC islands at φ 1, when the MGA is valid. For spheroid SC islands, when β = 1, Equations (15)- (17) and Equation (1) give the same result: ∆σ x /σ x ≈ 2φ, ∆σ z /σ z ≈ γ 2 φ/η ln 2γ/e √ η . For β = 1 the relative excess conductivity ∆σ x /σ x in Equation (15) is greater than in Equation (1) at the same φ by a factor of (1 + β) 2 /4β, which considerably exceeds unity at β − 1 ∼ β. For β 1 the increase ∆σ x /σ x ≈ φβ/2 ∝ β, which has an evident physical interpretation: thin elongated SC inclusions of random orientation give the excess conductivity almost as if they were spheroid with the largest dimension a y = βa x , but the volume fraction φ in this case is smaller by a factor of a x /a y = β −1 . However, ∆σ z /σ z decreases at β 1 by the factor β −1 . This is also evident: the increase of a y does not affect ∆σ z but increases the SC volume fraction φ ∝ β. Hence, at the same φ, ∆σ z ∝ β −1 .

Experimental
The experimental results for the excess conductivity above T c and for diamagnetic response in bulk samples are given in Figures 2-4 of Ref. [32], and we do not show these data here. Nevertheless, we re-analyze these data below, taking the expected size of SC islands into account. Here we show the measured R zz (T) curves for thin samples, which may help to estimate the size of SC islands.
In Figure 2c of Ref. [32] the results for R xx (T) and R zz (T) measurements in the microbridge of thickness L (0) z ≈ 200 nm is shown. One sees that the SC transition temperature T c itself is higher when determined from R zz (T) than from R xx (T). Similar T c anisotropy was reported in Ref. [37]. Below we explain this effect, analyze how this T c depends on sample thickness, and how this dependence can be used to extract information about the SC domains.
In Figure 3 we show the measured temperature dependence of normalized resistance R zz (T)/R zz (T = 15 K) in several microbridge samples of the same in-plane size 2 × 2 µm 2 but of different thickness, indicated in the figure legend for each curve. The normalization temperature T = 15 K was chosen because (i) we expect negligible volume fraction and the corresponding effect of SC domains at T > 15 K, and (ii) the R zz (T)/R zz (T = 15 K) curves at T > 15 K indeed coincide, as evidenced from Figure 3. The microbridge thickness for thicker samples is estimated visually from the SIM image of FeSe microbridge (overlap structure) oriented along the interlayer c axis, as shown in Figure 2b of Ref. [32] or in Figure 1a above. Therefore, we take the first microbridge thickness L (1) z ≈ 300 nm with an error about 10%. The black curve in Figure 3 shows R zz (T)/R zz (T = 15 K) in this microbridge. The SC transition temperature, corresponding to a 50% drop of resistance R zz (T), is about T c (50%) ≈ 8.5 K for this sample, while a 90% drop of R zz (T) happens at T c (90%) ≈ 8 K, which are only very slightly higher than T c determined from the in-plane resistance R xx (T) or from R zz (T) in bulk samples. This indicates that L (1) z ≈ 300 nm d z , and this sample almost behaves as a bulk one for the SC onset.
The in-plane resistance R xx (T)/R xx (T = 15 K) (not shown here) was measured only for larger samples, ∼ 1µm thick, as in Figure 2a of Ref. [32]. It is quite close to that in Figure 2c of Ref. [32] and, more importantly, to the R zz (T)/R zz (T = 15 K) curve in this microbridge of thickness L (1) z ≈ 300 nm. If normalized resistivity along two axes is similar, R xx (T)/R xx (T = 15 K) ≈ R zz (T)/R zz (T = 15 K), from symmetry arguments one may conclude that the mean aspect ratio of SC domains γ ≡ d z /d x ≈ L z /L x . This symmetry insight is confirmed by our percolation calculations for the spheroid SC domains below. For our FeSe sample it would give γ ≈ 0.15.  For thinner samples in our experiment T c determined from R zz (T) is higher, while T c determined from R xx (T) in large samples almost does not change. For each of the thinner microbridges m the thickness L (m) z was estimated according to: because for microbridges of the same in-plane area 2 × 2 µm 2 the measured interlayer resistance is proportional to the microbridge thickness L z . Unfortunately, this method of measuring microbridge thickness L z has an error, increasing with the decrease of L z , because L z may slightly vary along the microbridge area 2 × 2 µm. This approach may underestimate L z by 10-20%, especially, for thinnest microbridges. Therefore, we take L z = 50 nm for our preliminary percolation calculations in the next subsection (see Figure 4a).  illustration showing that the current percolation along the sample thickness is easier than along the sample length. Circular SC islands (blue) with diameter d = 0.4 are randomly distributed inside a rectangular sample (yellow) of dimensions 7 × 2, forming SC channels between contact electrodes (black).
From Figure 3, we see that the SC transition temperature T c strongly increases, when the sample thickness L z decreases: from T c ≈ 8 K at L (1) z ≈ 300 nm to T c ≈ 12 K at L z ≈ 40 nm. Note that L z ≈ 40 nm is still much larger than the in-plane SC coherence length ξ 0x ≈ 5 nm ξ 0z , so that the surface effects should not be important. We attribute this T c anisotropy to a heterogeneous SC onset and different percolation thresholds via SC domains in different directions for very thin samples, as described in the next subsection. Figure 4a shows the calculated probability p of current percolation along the in-plane x and out-of-plane z axes via spheroid SC domains as a function of SC volume fraction φ for two different domain heights, d z = 20 nm and d z = 5 nm, in a sample of dimensions 2 × 2 × 0.2 µm 3 , close to our experiment. The aspect ratio a z /a x = 0.62 of spheroid SC domains was chosen in agreement with Ref. [33]. Although we do not know exact domain shape and size, and the aspect ratio a z /a x = 0.62 is corrected in the next section, we make several important conclusions from this calculation. First, (i) the percolation probability along the shortest sample dimension z is indeed much higher than along the other two directions, which explains the observed anisotropic SC transition temperature T c in thin FeSe microbridges. This result has a simple explanation: the percolation along the shortest sample dimension (thickness) requires a much smaller number of SC domains than along the longest dimension (length), as illustrated in Figure 4b. Second, (ii) the effect of T c anisotropy depends strongly on the SC domain size d z as compared to sample thickness L z . The agreement with experiment is better for a larger domain size d z = 20 nm than for a smaller d z = 5 nm, suggesting an approximate average size of SC domains d z ∼ 20 nm. Third, (iii) the volume fraction φ c of SC domains, required for current percolation in the thinnest sample, is still rather high: φ c ∼ 0.2. From Figure 4a we find the sample-averaged percolation threshold φ c as corresponding to percolation probability p = 1/2.

Theoretical Analysis and Discussion
The SC volume fraction φ c ≈ 0.2, corresponding the percolation threshold along z for thinnest sample in Figure 4, gives an estimate of the SC volume fraction at the SC transition temperature T c ≈ 12 K in this thinnest sample. The φ(T c ≈ 12 K) ≈ 0.2 found in this way is much larger than the SC volume fraction φ(T = 12 K) < 10 −2 proposed in Refs. [32,33] basing on diamagnetic response data (see Figure 4d of Ref. [32]). A thinner sample, a larger size of the SC domains, or their elongated shape with a random orientation along x or y reduces φ c , but still keeps it large enough. Note that SC fluctuations can only enhance the diamagnetic response, further enhancing this discrepancy. This discrepancy is probably related to the assumption that the size d x of SC domains is larger than the SC penetration depth λ, which was made in the analysis of experimental data on diamagnetic response in Refs. [32,33]. In FeSe, the in-plane λ(T = 0) ≈ 400 nm and increases to ∼ 650 nm at T ≈ T c = 8 K, as observed from H c1 measurements (see Figure 6d of Ref. [47] At d λ, this simplifies to: If the SC domain length d y < λ, one can estimate the diamagnetic response as a contribution from small SC spheres of volume fraction φ (see Equation (8.22) of Ref. [44] and Equation (17) of Ref. [34]): where the demagnetizing factor of SC islands n 1 is small because of their oblate shape, and the factor (1 − n) can be omitted. We see that Equations (20) and (21) produce similar results, differing only in a numerical coefficient ∼ 1.
The experimental data on diamagnetic response in Figure 4c of Ref. [32] give ∆χ(T = 12 K) ≈ 0.7 · 10 −2 /4π. The assumption d λ, implicitly made in Ref. [32], means that instead of Equations (19)-(21) the SC volume fraction φ and diamagnetic susceptibility ∆χ are related by ∆χ ≈ −φ/4π(1 − n), which gives a strongly underestimated SC volume fraction φ 1 (T), as shown in Figure 4d of Ref. [32]. In particular, it gives φ 1 (T = 12 K) ≈ 0.7 · 10 −2 , which is much smaller than the value φ(T = 12 K) ≈ 0.2, expected from the SC percolation threshold along z axis for the thinnest sample in Figures 3 and 4. The difference between φ(T) and φ 1 (T) can be used to estimate the SC domain size. According to Equation (19), φ(T = 12K) ≈ 0.2 and 4π∆χ(T = 12 K) ≈ 0.7 · 10 −2 give an estimate d x /λ ≈ 0.65. Spheroid domain shape, according to Equation (21), gives a smaller domain diameter d x /λ ≈ 0.3, which better agrees with nematic domain width d n ∼ 200 nm. Note that, according to the BCS theory [44], λ(T) diverges at T → T c , but the H c1 measurements give finite λ ≈ 800 nm even at T = 9 K [47]. The penetration depth λ(T) is, therefore, poorly defined for SC domains at T > T c . If for our estimates we take λ ≈ 800 nm, corresponding to H c1 measurements at T = 9 K [47], we than obtain d x ∼ 0.3λ ∼ 240 nm. Taking λ = λ(T c ) ≈ 650 nm [47] gives d x ∼ 200 nm. This SC domain size d x slightly exceeds the average nematic domain width d n ∼ 100 − 200 nm but is much less than the length of nematic domains. The inequality d x > d n for averaged domain width is not too surprising. It can be explained by: (i) a significant fraction of wide SC and nematic domains with a width of d n 200 nm, which, due to their large size, make the main contribution to the diamagnetic response; (ii) the contribution of rare SC clusters consisting of several Josephson-coupled domains; (iii) additional diamagnetic response from SC fluctuations.
A weaker diamagnetic response ∆χ(T) of small SC domains, given by Equations (19)-(21), corrects to a higher value the SC volume fraction φ(T) extracted from ∆χ(T). It also corrects the estimated aspect ratio γ = a z /a x of SC domains to a smaller value than a z /a x ≈ 0.62 proposed in Ref. [33]. From Equation (16), neglecting the weakly dependent logarithmic factor, one obtains ∆σ z /σ z ∝ φγ 2 /β. At a fixed measured excess conductivity ∆σ z /σ z , this gives γ ∝ β/φ. Hence, the ∼ 28 times increase of estimated φ, from φ(T = 12 K) ≈ 0.7 · 10 −2 to 0.2 decreases γ ∝ 1/ √ φ about 5.3 times to γ = a z /a x ≈ 0.12 as compared to the earlier proposed [33] value. The corresponding percolation calculations for the sample geometry as in our experiment and the SC domain size d z ≈ 20 nm with γ = 0.12 are shown in Figure 5. These calculations suggest the percolation threshold φ c ≈ 0.12 rather than φ c ≈ 0.2 for the thinnest sample of L z ≈ 40 nm, where T c ≈ 12 K is determined from R zz (T) (see Figure 3). Hence, it corresponds to φ(T = 12 K) ≈ 0.12 in FeSe. This slightly modifies the estimate of γ ∝ 1/ √ φ to γ ≡ a z /a x ≈ 0.15. The corresponding percolation calculations for γ ≈ 0.15, d z ≈ 20 nm and for samples of different thicknesses, as in our experiment, are shown in Figure 6. The results of these calculations are in good agreement with the experimental data on anisotropic T c for FeSe microbridges, shown in Figure 3. In particular, in Figure 6, T c is almost isotropic for L z = L (1) z ≈ 300 nm, while for smaller L z T c is anisotropic and significantly higher if taken from R zz (T) rather than from R xx (T) curve, in agreement with Figure 3. Thus, our model of SC domain shape γ = a z /a x ≈ 0.15 and size d z ≈ 20 nm now agrees with the available combined experimental data on anisotropic SC excess conductivity above T c in bulk samples [32,33], on the anisotropic SC transition temperature in thin FeSe microbridges (see Figure 3), and on diamagnetic response in FeSe above T c [32]. . Probability p of current percolation along the in-plane x (red dashed curves) and the out-ofplane z axes (dotted blue curves) via SC domains of spheroid shape with aspect ratio a z /a x = 0.15 and height d z = 20 nm as a function of SC volume fraction φ in a sample of x-y area 2 × 2 µm 2 , calculated for four different sample thicknesses L z = 40 nm (a), L z = 125 nm (b), L z = 235 nm (c), and L z = 300 nm (d). The percolation probability is isotropic for L z = 300 nm.
The current percolation and T c anisotropy in thin samples are very sensitive to the aspect ratio γ = a z /a x of SC domains, which allows its accurate measurement. Indeed, for the aspect ratio γ ≡ a z /a x = 0.12, as in Figure 5, the percolation threshold φ c is isotropic for the sample thickness L (2) z ≈ 235 nm, but not for L (1) z ≈ 300 nm, as in Figure 6 for slightly different γ = 0.15 and in our experiment. This suggests a new precise method for measuring the averaged aspect ratios d x : d y : d z of SC domains deep inside the sample, which is not accessible by STM or other surface measurements. Indeed, if small samples are fabricated, only several times larger than the expected domain size, then an anisotropic SC transition temperature T c to almost zero resistance should be observed, as in Figure 3. However, if the sample aspect ratios L x : L y : L z match the average aspect ratios d x : d y : d z of the SC domains, this T c anisotropy should disappear even if the sample size remains small, L i 10d i , as for the sample with thickness L z = 300 nm and L x = L y = 2 µm in our experiment. The aspect ratios of this sample give the aspect ratios of the SC domains.
Our experimental data in Figure 3 and percolation calculations in Figure 6 suggest the aspect ratio value γ ≡ d z /d x ≈ 0.15 of SC domains in FeSe. This value agrees with the one obtained from the comparison of measured [32] excess conductivity and diamagnetic response at T > T c in bulk samples if the SC domain width d x ∼ 200 nm is comparable to the nematic domain width d n in FeSe. Note that the corresponding domain height d z = γd x ∼ 30 nm falls within the required interval 10 nm< d z < 40 nm, where, according to our percolation calculations, the T c anisotropy is significant for the sample thickness 40 nm< L z < 300 nm in our experiment. If the SC domains are not spheroid and have elongated shape with β = a y /a x > 1, the estimated ratio a z /a x grows ∝ β. For β = 5 we obtain γ = a z /a x ≈ 0.25. Percolation calculations can also be performed for this case. However, a direct experimental study of the shape and size of SC domains using STS measurements would be very useful in confirming our semiphenomenological predictions about the geometry of SC domains.
The methods employed above are also useful for many other compounds with heterogeneous superconductivity onset. For example, in FeS the spatial inhomogeneity has much larger length scale than in FeSe; the domain size d ≈ 35 µm, far exceeding the SC penetration depth λ, was observed with submicrometer resolution spatially resolved ARPES (µ-ARPES) in FeS [48]. It is noteworthy that in FeS there is neither a "nematic" phase transition to an orthorhombic lattice, which occurs in FeSe at T n ≈ 90 K and leads to a domain structure, nor magnetic ordering, as in FeTe below T AFM ≈ 75 K. Probably, the spatial inhomogeneity in FeS arises from the interplay of different types of electronic ordering, similar to organic superconductors.

Conclusions
In Figure 3, we present the experimental data on the temperature dependence of resistivity R(T) in thin FeSe samples, produced by cutting the bulk samples using FIB in the form of the microbridges shown in Figure 1. The SC transition temperature T c strongly increases as the sample thickness decreases. We explain this effect by calculating the percolation probability via the SC islands as a function of SC volume fraction φ in different directions. The anisotropy of the percolation threshold arises from the finite sample size and its flat shape. The thinnest microbridges are only a few times thicker than the SC domain size, and, in contrast to large samples, the percolation probability along the shortest sample dimension (thickness) is much higher than along its length (see Figure 4). Similar effects appear in organic superconductors [35]. Our calculations of percolation probability for the relevant sample shape and size, combined with our experimental data in Figure 3, suggest several important properties of the SC onset in FeSe. (i) The SC onset in FeSe is spatially heterogeneous and proceeds in the form of isolated SC domains, which become phase-coherent at lower temperature T c , corresponding to SC transition of the entire sample. This is similar to many other high-Tc superconductors [1][2][3][4]; (ii) The SC domain height d z is several times smaller than the thinnest microbridge thickness L z ≈ 40 nm: d z ∼ 20 nm; (iii) The SC volume fraction at T = 12 K in FeSe is rather large, φ ∼ 0.1.
The small SC domain size d z means that the in-plane SC domain size d x ∼ 7d z is smaller than the penetration depth λ of magnetic field to FeSe superconductor. Hence, the estimates of the temperature dependence of the volume fraction φ 1 (T), obtained from the diamagnetic response and shown in Figure 4d of Ref. [32] or in Figure 4 of Ref. [33], is underestimated by ∼ 20 times. As a result, the SC domain aspect ratio a z /a x in Ref. [33] is, probably, overestimated ∼ 5 times. The combined new analysis suggests that the averaged aspect ratio a z /a x ≈ 0.15, and the corresponding in-plane SC domain size d x ≈ 100 − 200 nm, which is comparable to the nematic domain width d n in FeSe. Thus, the hypothesis that SC domains are inside the nematic domains is consistent with the combined transport and diamagnetic experiments in bulk FeSe [32], as well as with our R zz (T) measurements in thin FeSe samples and corresponding percolation calculations. Notably, the proposed method of estimating the averaged SC domain aspect ratio a z /a x from the anisotropy of SC transition temperature T c in finite-size samples turned out be very precise-with an accuracy ∼ 10%-and is applicable to many other heterogeneous superconductors.
We have also generalized the analytical formulas for conductivity in heterogeneous anisotropic superconductors in the case of elongated SC domains of two perpendicular orientations with equal volume fractions, corresponding to the nematic domain structure in various Fe-based superconductors (see Equations (15)-(17)). These formulae are useful for the analysis of anisotropic excess conductivity at T > T c in order to obtain useful information about the shape and volume fraction of SC domains.
In this paper, we focused on FeSe, although our method and discussion are applicable to other materials, including various cuprate and Fe-based high-T c superconductors, organic metals, and other compounds. The obtained knowledge about the SC domains, electronic structure and SC properties of FeSe during the heterogeneous SC onset may help us to better understand the SC mechanisms and the properties of Fe-based superconductors, as well as to search for novel high-T c superconductors.

Abbreviations
The following abbreviations are used in this manuscript: