Supplemental Information : Unconventional localization of electrons inside a nematic electronic phase

The exfoliation of rectangular-shaped flakes onto pre-patterned contacts leads to samples with rather non-ideal geometries. The current flow is inevitably quite inhomogeneous and the contacts extend a long way underneath the flake, violating the assumptions of the van der Pauw approach commonly used for irregularly-shaped samples. As a consequence the longitudinal and Hall resistivities have been estimated from finite difference solutions of the transport equation for a realistic sample and contact geometry. For our two-dimensional problem the equation to be solved is given by:

Supplemental Information: Unconventional localization of electrons inside a nematic electronic phase

Geometrical corrections
The exfoliation of rectangular-shaped flakes onto pre-patterned contacts leads to samples with rather non-ideal geometries. The current flow is inevitably quite inhomogeneous and the contacts extend a long way underneath the flake, violating the assumptions of the van der Pauw approach commonly used for irregularly-shaped samples. As a consequence the longitudinal and Hall resistivities have been estimated from finite difference solutions of the transport equation for a realistic sample and contact geometry. For our two-dimensional problem the equation to be solved is given by: where ⃗ E is the 2D electric field and ⃗ J the 2D current density,ẑ is unit vector perpendicular to the sample and is the direction of the applied magnetic field, ⃗ B=µ 0 ⃗ H; σ xx and ρ xy = (B)/(n 2D e) are the sheet conductivity and Hall resistivity of the flake respectively, with n 2D the carrier concentration. Fig. S1 illustrates an example used to solve Equation S1 in terms of the potential, V , on a 144 × 66 site grid. We impose the following boundary conditions along the left and right edges of the flakes where J x = 0, the two components of the electric field satisfy E x = −(σ xx ρ xy )E y , whereas along the top and bottom edges where J y = 0, the two components of the electric field satisfy E y = (σ xx ρ xy )E x . Once the solution has converged the current passing through the device, I, is calculated by summing the following expression for J x down a vertical line through the centre of the device.
We assume that the current leads are strongly coupled to the sample and arbitrarily fix V 1 = 0V and V 4 = 1V everywhere in the flake directly above them. In contrast we assume that the voltage contacts are only weakly coupled to the flake reflecting the strong anisotropy of our layered materials (a factor larger than 4 which increases with the reduction in the flake thickness [1]), and calculate the potential at each of the voltage leads as the average of the potential in the flake directly above it. It is then straightforward to show that the scaling factors (which are typically in the range 1-3) that must be applied to obtain the longitudinal resistivity, ρ xx = 1/σ xx , and the Hall resistivity, ρ xy , from the experimentally-measured longitudinal resistance, R xx , and Hall resistance, R xy , are: Here g = σ xx ρ xy =0.01 is a parameter that is used in the numerical calculation to effectively define the magnetic field used to solve for the Hall voltage. Using the definitions in Equations S3 then the conversion factors for each resistivity component are G xx = I/(V 3 − V 2 ) and G xy = I · g/(V 2 − V 6 ). Since Equation S1 yields a Hall voltage that is linear in magnetic field, H, the scaling factor for the Hall resistivity calculated in Equation S3 does not depend on this.

Two-band model
Considering a multi-band system in which current is applied along the x in-plane axis of a sample with a magnetic field applied along the z out-of-plane axis, the total resistivity is given by ρ = i ρ −1 i −1 assuming parallel network resistor which in an applied magnetic field B leads to: where the conductivity σ i = |n i eµ i | and the Hall coefficient R i = −1/n i e contain the carrier densities n i and mobilities µ i of i number of bands.
For a two-carrier system, the different components of the resistivity tensor are given by: and FeSe in the tetragonal phase can be described as a compensated two-band system, and we assign the conductivities to correspond to a hole σ h and electron band σ e , respectively. Compensation requires that R 1 = -R 2 and n = n e = n h leading to a simplified form of the expression above and By fitting simultaneously the magnetic field dependence of both the longitudinal magnetoresistance, ρ xx , and the transverse Hall component, ρ xy , to the above equations, the carrier density n and the two mobilities µ e and µ h can be extracted.

The mobility spectrum
The mobility spectrum has been developed to eliminate the need for making a priori assumptions on the transport parameters in a multicarrier system. The mobility spectrum analysis has been used extensively in multi-band semiconductors, in multi-band topological semi-metals [2] and more limited in multi-band iron-based superconductors [3,4]. The conductivity of a multi-band system is given in terms of a mobility spectrum, based on a method described in Ref. 5 This work is part of a separate report related to the mobility spectrum of FeSe 1−x S x [6].
where s(µ) = eµn(µ) is the mobility spectrum of zero-field conductivity and the charge carrier densities are functions of mobility (with electrons defined to have negative µ and n(µ)). Thus the fundamental object of interest to identify the properties of charge carriers in a material becomes the conductivity spectrum s(µ), or equivalently the carrier-density spectrum n(µ).
The mobility spectrum constructed for t=125 nm in Figure 2(a) is expanded in different components in Fig. S9. The initial parameters proposed by the mobility spectrum, together with compensation of the charge carriers, are used to extract the discrete fitted parameters in Figs. 2(d),(e) and (f). The behaviour of the mobilities and carrier densities as a function of temperature are consistent between the two approaches. The mobility spectrum could be a powerful tool to identify the conduction in multiple band systems with different mobility which appear as distinct peaks in s(µ), as long as the condition µB < 1 is satisfied. Its fundamental description of the electrical transport does not require advance knowledge about the band structure or the scattering mechanisms. A distribution of relaxation times results in a distribution of mobilities and broadening of the corresponding peak in s(µ). Therefore, the shape of the peak in s(µ) reveals not only an average of relaxation times ( ⟨τ ⟩ 2 /(τ ) 2 ), but the actual distribution of relaxation times [5].
Approximations are often used to extract the form of the mobility spectrum for a material having a finite data set of conductivity or resistivity measurements as s(µ) is never fully constrained on the infinite set of basis functions. Discrete approximations can often add ghost peaks in the spectrum, as well as the introduction of bias with a necessary predetermined range of interest for mobility. Other techniques, used in the case of iron-based superconductors, fit analytic forms to the data which can be directly transformed [3,4]. These approaches have the disadvantages of losing some fidelity of the data and require an analytic continuation of the fitted function to infinite magnetic field values.

The Lifshitz-Kosevich equation
The amplitude of quantum oscillations can be described by the Lifshitz-Kosevich equation [7,8]: where R T , R D , and R S are damping terms. The first sum extends over all extremal Fermi surface areas A k,i perpendicular to the applied field, while the second over all p harmonics of the fundamental oscillation frequency. These oscillations are periodic in 1/B with a frequency which is determined by orbits on the Fermi Surface that enclose locally extremal momentum space area, where F i = ℏ 2πe A k,i . The first damping term, R T , accounts for the thermal broadening from the Fermi-Dirac distribution with temperature and is given by: where m * is the quasiparticle effective mass. R T depends on the ratio X ∝ k B T /ℏω c , where ω c = eB/m * is the cyclotron frequency. These two equations contain the temperature dependence of the amplitude of oscillation, and by fitting data to Equation S12 the effect mass of m * can be extracted.
FIG . S2. Estimation of the quantum scattering time from Dingle plots. Dingle plots to extract the Dingle temperature, TD, from the slope and accounting for the RT and RD field dependence for (a) a t = 58 nm device (TD ∼ 4.5(4) K) and (b) a bulk single crystal S1 (TD ∼ 1.6(2) K).
The second damping term, R D , is accounts for the impurity scattering of the electrons and is given by where τ is the scattering time, m b is the band mass, and T D = ℏ/2πk B τ is the Dingle temperature. This equations is determined by the impurity scattering rate which acts to exponentially dampen the amplitude of the quantum oscillation in magnetic field.
As the band structure calculations that provide the band mass cannot capture the Fermi surface of FeSe correctly [9], we use m * = m b to determine the Dingle temperature. The values of the bulk were taken from previous reports [10,11]. The third damping term, R S , accounts for the Zeeman splitting of Landau levels.
FIG. S3. Resistivity in zero-magnetic field of thin flakes devices. Temperature dependence of the zero-field resistivity for four thin flake devices with different thicknesses (a) t =125 nm, (b) t = 58 nm, (c) t = 28 nm and (d) t =14 nm. The insets show the low temperature superconducting transition region and highlights that at low temperature the resistivity has a rather linear dependence for most of the thin flakes samples. Additional resistivity curves for other thin flakes were previously reported in Ref. [1]. A small increase in resistivity of the thinnest flakes is detected in (c) and (d), which is often interpreted as a signature of Anderson localization due to disorder. However, insulating behaviour was detected in the ultra-thin limit of FeSe flakes below 9 nm [1,12] and the number of charge carriers remain unchanged upon reducing thickness (see Fig. 2(d) and Ref. [12]).    [14] with k00 = 0.1Å −1 , the interplane parameter of k10 = 0.01Å −1 , and the in-plane two-fold symmetric parameter (a) k02=0, (b) k02=0.01, (c) k02=0.03, (d) k02=0.05Å −1 . (bottom panels) The evolution of the Fermi surface and the scattering length of the in-plane four-fold symmetric parameter (e) k04=0, (f) k04=0.01, (g) k04=0.03, (h) k04=0.05 A −1 . The Hall conductivity is linked to the 2 × A ℓ , area swept by the scattering path length, ℓ k . In the cases (g) and (h) the Hall conductivity can change sign as compared with the other cases. The scattering time is assumed isotropic, τ = 1 ps [15].