Evidence for unidirectional nematic bond ordering in FeSe

The lifting of $d_{xz}$-$d_{yz}$ orbital degeneracy is often considered a hallmark of the nematic phase of Fe-based superconductors, including FeSe, but its origin is not yet understood. Here we report a high resolution Angle-Resolved Photoemission Spectroscopy study of single crystals of FeSe, accounting for the photon-energy dependence and making a detailed analysis of the temperature dependence. We find that the hole pocket undergoes a fourfold-symmetry-breaking distortion in the nematic phase below 90~K, but in contrast the changes to the electron pockets do not require fourfold symmetry-breaking. Instead, there is an additional separation of the existing $d_{xy}$ and $d_{xz/yz}$ bands - which themselves are not split within resolution. These observations lead us to propose a new scenario of"unidirectional nematic bond ordering"to describe the low-temperature electronic structure of FeSe, supported by a good agreement with 10-orbital tight binding model calculations.

The search for the understanding of unconventional superconductivity in the Fe-based systems has lead to a focus on the origin and nature of the ordered phases found in close proximity to the superconducting phase. Particular attention has been drawn to the "nematic" phase [1], where the four-fold symmetry of the lattice is broken at a temperature T s which is higher than the striped antiferromagnetic ordering temperature T N in some materials [2]. There has been long discussion about whether these two transitions are both magnetic in origin, or whether T s corresponds to a separate orbital instability [1]. FeSe is an exceptional case, since it undergoes a structural transition at T s = 90 K without long range magnetic order at any temperature, enabling detailed study of the symmetrybroken state [3][4][5][6][7][8][9][10], and has therefore attracted much theoretical attention [11][12][13][14][15][16][17].
Angle-Resolved Photoemission Spectroscopy (ARPES) measurements of FeSe [18][19][20][21][22][23][24][25][26][27] provide direct experimental access to the evolution of the electronic structure through T s , which can be linked to theoretical models of the underlying symmetry-breaking order. Most previous ARPES studies have paid particular attention to the electron pockets around the M point, and inferred a ∼50 meV splitting of d xz and d yz bands, similar to previous claims in NaFeAs [28] and Ba(Fe 1−x Co x ) 2 As 2 [29]. In this scenario, which could be interpreted as a ferro-orbital ordering [22], the bands of primarily d xz and d yz character are degenerate exactly at M in the high temperature phase [30], but split below T s . Then, the dispersions observed by ARPES all arise from d xz or d yz bands, taking into account the presence of twin domains in the sample, with the assumption that outer electron band with d xy character does not contribute [22]. While alternatives to the ferro-orbital ordering scenario have been recently proposed [11,25,[31][32][33], the splitting of d xz and d yz bands at the M point has been considered a hallmark of the nematic phase, until now [20][21][22]24].
In this letter, we present a high-resolution ARPES study of the evolution of the electronic structure of bulk FeSe through the nematic transition at T s . We use curvature analysis and a new procedure to fit the Energy Dispersion Curves (EDCs) at the M point to extract band positions. We also observe the k z dependence of the electron pockets, and present new Fermi surface maps covering the whole Brillouin zone above and below T s . While the hole pockets of FeSe do undergo significant symmetry-breaking distortions below T s , the changes to the electron pockets arise from an additional separation of existing d xz/yz and d xy bands, without any resolvable lifting of d xz -d yz band degeneracy at the M point. These observations exclude the ferro-orbital scenario, and place new strong constraints on the nature of the orbital ordering. Finally we use a ten-orbital tight binding model to propose that the nematic state of FeSe is characterised by rotational-symmetry breaking within Fe-Fe hoppings, which may be described as a unidirectional nematic bond ordering.
High-quality single crystal samples of FeSe were grown by the chemical vapor transport method [3,22,34]. ARPES measurements were performed at the I05 beamline at the Diamond Light Source, UK.
In Fig. 1a-d) we present temperature-dependent ARPES data for a high symmetry cut centred on the M point, also showing the EDCs at M. At low temperatures (Fig. 1d)), these EDCs display two prominent peaks with a separation ∆ M = 50 meV, previously attributed to orbital splitting in the nematic phase, i.e. ∆ F O M = E yz − E xz . However, curvature analysis [35] enhances weak features in the data and provides a different perspective, as shown in Fig. 1g-j). Above T s , in Fig. 1g) both the expected d xz and d yz dispersions and also sections with d xy character are observed. By comparing Fig. 1g-j), we observe the similarity in the dispersions above and below T s , and no extra features arise. Therefore here we propose a new scenario, where the low temperature dispersions simply consists of the d xz/yz and d xy bands as expected above T s , but with an increased separation between them. We assign the peaks in the EDCs to the two bands, i.e. ∆ M = E xz/yz − E xy , equal to 50 meV at low temperatures but also finite above T s . Within experimental resolution (11 meV FWHM for the relevant peak in Fig. 1d)) we do not detect any subsequent splitting between d xz and d yz bands, although since this is allowed by symmetry in the orthorhombic phase a small splitting may occur as a secondary effect.
This new scenario can be tested by the extraction of band Binding Energy (eV)  positions as a function of temperature from the EDCs at the M point. However while the separation ∆ M of the two features in the EDC is unambiguous at low temperatures, it becomes more difficult to define at higher temperatures, as features become broader. Here we take a new approach: we fit the EDC with two asymmetric pseudo-Voigt functions and the Fermi function at 61 K where the peaks are still distinct. Then, at higher temperatures we fit using the same peak profiles, only allowing the peak positions to vary (see Supplemental Materials, SM). We find that the fitted peak positions are separated by 20 meV even above T s , but increase substantially when the system enters the nematic phase (Fig. 1e,f). The sharp upward shift of the d xz/yz bands at T s is the strongest feature, although the d xy band position also adjusts downwards below T s , leading to a total separation of 50 meV at 10 K. In Fig. 1m-p) we present new measurements of the electron pockets using a photon energy of 56 eV, where the outer d xy electron dispersion is already clearly visible in the data above T s even without curvature analysis. The 56 eV photon energy is chosen to correspond to the A point at the top of the Brillouin zone where the warped quasi-2D electron pocket is largest, giving the best momentum-resolution of features, as can be seen in Fig. 1l). Further details of the photon-energy dependence are presented in SM, where it is also shown that ∆ M is independent of k z . In Fig. 1q,r) we extract the temperature-dependence of both k F vectors observed, from peak fitting of the Momentum Distribution Curve (MDC) at the Fermi level (SM). The k F values are well-defined at all temperatures and also demonstrate that the d yz and d xy dispersions are separate above T s and undergo additional separation below T s .
In Fig. 2 we present further support for the new interpretation including the observation of the d xy electron band. These ARPES data are obtained with the scattering plane at 45 • to the Fe square lattice to mitigate matrix elements effects [36], again using 56 eV photon energy. At 100 K, in the tetragonal phase above T s , the Fermi surface maps (Fig. 2a,e) reveal essentially the whole structure of the expected electron pockets at the A point, including the outer pocket with d xy character. Since all expected bands contribute at high temperature, we conclude that the whole structure of the low temperature Fermi surface is also observed. The comparison between Fig. 2e) and i) demonstrates that the low-temperature Fermi surface is elongated compared to the high temperature case, but retains the same basic structure and symmetry. Thus the electron pockets retain fourfold symmetry and the bands from both structural domains will overlap in experiments. The evolution of the electron pockets through T s contrasts strongly with the behaviour of the hole pockets. At the Z point ( Fig. 2d), above T s the inner hole pocket just crosses the Fermi level, making a small 3D pocket [22,24], and there is a ∼20 meV splitting due to spin-orbit coupling at the Z point between the d xz/yz bands. However below T s there is an extra splitting of ∼15 meV associated with a d xz/yz orbital polarisation, in addition to the spin-orbit splitting. This pushes the inner band completely below the Fermi level (Fig. 2h). Therefore there is a single hole band at low temperatures which displays an elliptical distortion, however due to the presence of twin domains below T s , two crossed ellipses are superposed in the ARPES data in Fig. 2b,c).
The breaking of fourfold rotational symmetry of the hole pockets below T s is therefore well-established. However, we have found symmetry-preserving changes to the electronic structure at the M point, which on first sight are difficult to reconcile with the tetragonal-to-orthorhombic phase transition at T s . Another constraint is that back-folded bands are not observed, suggesting that translational symmetry is preserved at T s . The challenge is to identify an orbital order parameter which globally breaks fourfold rotational symmetry but is consistent with the strong constraints provided by these observations. Ferro-orbital ordering could account for the symmetry-breaking at the Γ point, but requires a splitting of d xz/yz orbitals at the M point, which we have argued is not the case (see also SM). Moreover, ferro-orbital ordering is not consistent with the direction of distortion of the hole band, as revealed by ARPES measurements on detwinned crystals [25]. A Néel-type antiferro-orbital ordering would preserve the d xz/yz degeneracy at M, but it cannot explain the extra splitting observed between d xz/yz bands at the Γ/Z point [31]. The recently proposed d-wave bond nematic order predicts a splitting at the M point [31,32], as does the microscopic model of Ref. [33]. However, we suggest that a "unidirectional nematic bond ordering" is compatible with all the observations.
In Fig. 3 we present tight-binding model electronic structures with and without the proposed unidirectional nematic bond order. We use a 2D ten-orbital tight binding model including spin-orbit coupling [14,37] (SM). Within this model, we add an order parameter to the inter-site hopping block of the Hamiltonian as h = ∆ S (n yz − n xz ) cos(k x ), where n xz indicates the inter-site hopping operator distinct from on-site occupations n xz . This order parameter therefore describes a symmetry-breaking within the inter-site d xz/yz hopping terms on the Fe-Fe bonds in the x direction. This order parameter has the desired properties of giving an extra splitting in addition to the spin-orbit coupling splitting of d xz/yz bands at the Γ point, but a symmetric shift up of the d xz/yz bands at the M point without losing the degeneracy. It still globally breaks fourfold rotational symmetry yet it can also be simply shown that it does not break translational symmetry (SM). In Fig. 3b,c) we present the results of a calculation with ∆ S = 20 meV; we also include an adjustment of -10 meV to the d xy orbital which is motivated by the experimental results in Fig. 1 and may be required to maintain the charge balance of the system. Within this fairly simple model we obtain Fermi surfaces and dispersions which reproduce all qualitative features of the low-lying electronic structure. Therefore we suggest that the changes to the electronic structure of FeSe in the orthorhombic phase may be primarily described by a unidirectional nematic bond ordering.
Since the Fermi surface changes at the M point are not of a symmetry-breaking nature, and ∆ M is finite even above T s , there is an important question about how tightly the evolution of the electron pockets is linked to T s . However we have shown that there is a sharp increase in ∆ M which onsets exactly at T s (Fig. 1j). Additionally, the deviation in k F values of both the inner and outer electron pocket branches also follows a sharp order-parameter-like behaviour which onsets at T s , and this deviation of the inner pocket was shown to behave as an order parameter of the structural transition across the FeSe 1−x S x series [24]. This indicates that the changes at the M point are still fundamentally linked to the orthorhombic lattice distortion. This may be understood since in the unfolded 1-Fe Brillouin zone, the electron pockets with d xz and d yz character are located in different parts of k-space, and therefore the cos (k x ) term ensures that they shift symmetrically, such that degeneracy is not lost in the folded 2-Fe zone (SM). Therefore the apparently non-symmetry breaking band movements at M are linked to an ordering which globally does break tetragonal symmetry concomitant with the tetragonalorthorhombic lattice distortion at T s .
Within the longstanding debate about the roles of orbital and spin degrees of freedom in Fe-based superconductors [1], FeSe is often considered as an example where orbital interactions may be dominant [3-5, 9, 22]. However we have excluded any significant ferro-orbital ordering, and furthermore we have found that the primary order parameter is bondcentered and that translational symmetry is not broken, in contrast to some proposals [11,38]. Intriguingly, a bond-centered ordering is also observed in cuprates [39], although in that case it has an incommensurate modulation. On the other hand, there are mixed reports on how relevant magnetic interactions are to the structural transition in FeSe [3][4][5][6][7][8][9]12]. The magnitude of energy shifts (e.g. 20 meV for the d xz/yz bands at M) are similar to the spin-orbit coupling value of 20 meV in FeSe [22], suggesting that a spin-driven scenario is not ruled out. Thus our data excludes several proposed orbital ordering scenarios, and instead points to the existence of a unidirectional nematic bond ordered phase in FeSe, breaking rotational but not translational symmetry, which is distinct from the known striped antiferromagnetic phase in other Fe-based superconductors.
In conclusion, we have presented a high-resolution ARPES study of FeSe, and argued that the symmetry-breaking distortion of the hole pockets at the Γ point but symmetrypreserving changes of the electron pockets at the M point below T s point can be explained by a unidirectional nematic bond ordering. These measurements provide a fresh perspective on the nature of the nematic phase of Fe-based superconductors. The extraction of band positions at the M point has been an important part of previous studies of FeSe, where it has been claimed that a split peak structure arises from a single peak (i.e. the degenerate d xz/yz bands at M) either at T s [20] or above T s [21,23]. In the main text we argued that this is not the correct scenario, and that the EDC above T s must arise from two separate features, even if they are not well-distinguished, and there is a well-defined transition at T s . Other studies have typically used a second derivative analysis to determine the band positions, but this can give conflicting results [21,23]. As described in the main text, for this study we decided to take a novel approach to the problem by directly fitting the EDC to extract band positions.
At 61 K (Fig. 4a), the peaks are clearly distinct, and here we perform a fit to the data with two asymmetric pseudo-voigt functions with additional background, multiplied by the Fermi function -that is, with 14 free parameters such that the peak profiles are very well fitted. However, once this fit is converged, at higher temperatures we use exactly the same parameters (amplitude, peak width, asymmetry parameters and background), except for an overall peak width term which is fixed to increase linearly with temperature to account for the experimental increase in peak width with temperature. We therefore extract the band positions above 61 K from a final fit which has only the peak positions free -i.e. a fit with only two free parameters. Representative fits are displayed in Fig. 4. We decided that the slight loss in fit quality is compensated by the greater confidence in the peak positions when these are the only free parameters. Indeed, this procedure gives a sensible and systematic temperature-dependence of band positions as shown in Fig. 1j of the main text, and is complementary to the result of the MDC fitting analysis. Note that this fitting method in principle could converge with both peaks at the same position and summed intensity above T s which would be compatible with ferroorbital ordering and is therefore unbiased between the two interpretations, but in fact the method converges to be consistent with the revised interpretation with a 20 meV splitting at the M point at higher temperatures between d xz/yz and d xy bands. In Fig. 5 we present MDC fit profiles of the high-symmetry cut through the M point at 56 eV, at representative temperatures, in order to show how the temperature-dependence of k F was determined in Fig. 1p) of the main text. For the fit procedure, we use two pairs of Lorentzian peaks with equal widths. We presented a similar analysis using 37 eV data in Refs. [22,24]. However at that photon energy the outer branch has much weaker intensity (Fig. 1 of the main text) and cannot be tracked up to T s .

Photon Energy dependence at M
Here we present further details of the photon energydependence of the electron bands at M point, at low temperatures in a twinned sample. The photon energy dependence of the hole pockets of FeSe has been previously reported [22], revealing a significant k z dependence of the outer hole pocket. At the M point, the Fermi surface shows a smaller warping effect, as the k F of the outer pocket shrinks from 0.193Å −1 (56 eV -A point) to 0.142Å −1 (42 eV -M point). This supports the analysis of quantum oscillation frequencies [10,22] where it was determined that the electron pocket was less dispersive in k z compared to the hole pocket. It is interesting to note that the energy splitting at the M point appears to be k z -independent, as shown in Fig. 6b). This indicates that the orbital ordering has no significant k z dependence. Thus the band positions at the M point do not change significantly with k z , but the effective mass of the dispersion of the outer electron band from the M point does rise with k z , giving rise to the warping of the Fermi surface. This can be directly observed by comparing the cuts at different photon energies shown in Fig. 6c-e), where the most substantial difference between them is the d xy band dispersion. Note that one previous measurement of the k z dependence at the M point did not find any significant k z dependence at the M-A point [27], which we find surprising. Taking into account the k z dispersions of both the hole and electron pockets, we determine that the inner potential (V 0 ) appropriate for FeSe is approximately 12.2 eV.

Tight binding model
The model: In the main text, we presented Fermi surfaces and band dispersions for a 10-orbital tight binding model, in the form originally proposed in Ref. [37]. The parameters of Ref. [14] were used as a starting point, with some adjustments made in order to better match the ARPES dispersions at high temperature. In practice we adapted the parameters until the low-lying bands at Γ and M are described well; however some details such as the dispersion of the d xy hole band are less well captured.
Motivation of order parameter: We now discuss how the inclusion of a unidirectional nematic bond order term in our tight binding model is required in order to understand the ARPES data. Motivated by the experimental data, we added the following term to the 10-orbital model Hamiltonian in the inter-site blocks: Here the x direction corresponds to a Fe-Fe direction and the longer a axis of the orthorhombic unit cell. We make the distinction here between n xz the on-site orbital number operator, and n xz which is the operator for inter-site hopping between d xz orbitals. We describe it as "unidirectional nematic bond order" since there is a cos(k x ) but not cos(k y ) term.
In the unfolded 1-Fe unit cell (i.e. 5-orbital model), the two crossed ellipses that occur at the M point in the crystallographic 2-Fe unit cell are unfolded to different positions in momentum space; the pockets containing d xz and d yz character exist at (0, ±π) and (±π, 0) respectively. Experimentally, we have determined that these two ellipses must distort symmetrically. This could not occur for an order param- eter defined homogeneously across the Brillouin zone, e.g. a ferro-orbital ordering, since the two pockets would distort in opposite directions due to their different orbital characters. Therefore we considered momentum-dependent order parameters, and in particular a cos(k x ) term ensures that the orbital splitting has different signs at (π, 0) and (0, π). Due to the distinction in the orbital character at each electron pocket, this makes band shifts at each point occur symmetrically, and when folding back to the M point of the 2-Fe Brillouin zone, the Fermi surface has become more elongated as in experiment but the degeneracy of d xz − d yz bands is not broken.
A cos(k y ) momentum-dependence could also give symmetric shifts at the M point, however this would lead to the opposite distortions of the hole pocket compared to measurements on detwinned single crystals [25]. Therefore we cannot rule out the existence of an additional bond nematic ordering with cos(k y ) modulation and significantly smaller magnitude than our primary (n yz − n xz ) cos (k x ) term. However experimentally the magnitude of the band shifts at M and the extra splitting at Γ are similar and therefore a single unidirectional order parameter is sufficient. Preservation of Translational Symmetry: A simple argu-ment can be made to show that this order parameter does not break translational symmetry. We can rewrite Eq. 1 as: h = ∆ S 2 (n yz − n xz )(cos(k x ) + cos(k y ))+ ∆ S 2 (n yz − n xz )(cos(k x ) − cos(k y )) (2) The first term here corresponds to "extended s-wave bond nematic order" [32] and the second is an allowed hopping term which doesn't break any symmetries (see also Fig. 7). Therefore the symmetries which are broken at T s are those of the extended s-wave bond nematic order, i.e. breaking rotational but not translational symmetry (thus it is truly a "nematic" order parameter). However although this decomposition is mathematically valid and noteworthy, it seems to us that our description of a unidirectional nematic bond order is more intuitively related to the underlying physics than using the description given by this linear combination. In particular, since the extended s-wave bond nematic order vanishes at the M point, this description would rely on symmetry-allowed hoppings to give the order-parameter-like behaviour at the M point [24], which seems unphysical.
Comparison with other orbital orderings: In Fig. 7 we present a selection of Fermi surfaces in the presence of alternative orbital ordering schemes, sugggested by various authors for describing the nematic state. The simplest case is onsite ferro-orbital ordering [14] in b) which is not viable since it breaks d xz − d yz symmetry at the M point. The d-wave bond nematic order [31,32] also fails due to the splitting at M and absence of splitting at Γ. The extended s-wave bond nematic order was introduced in [32] as a possible nematic order, but has no signature at the M point. In e) we show that the "allowed hopping terms" introduced above which break no symmetry could account for the evolution of the M point alone. Thus a linear combination of the extended s-wave bond nematic and the extra hopping could account for the low-temperature Fermi surface and are mathematically equivalent to our unidirectional nematic bond order. Finally we reproduce the unidirectional nematic bond order parameter from the main text.
On-site vs bond order: In a 5-orbital model, one could choose an order parameter (n xz − n yz ) cos (k x ) which would reproduce the experimental Fermi surface. Since the 5-orbital model does not distinguish between on-site number and intersite orbital hopping operators, this leaves ambiguity in the interpretation. However when using the full 10-orbital model, the order parameter must be placed in the inter-site hopping terms to reproduce the experimental Fermi surface. The interpretation of on-site stripe antiferro-orbital ordering (or any other pattern) is forbidden since this would require translational symmetry-breaking and additionally has the wrong periodicity.
The d xy band: Experimentally, the sections of Fermi surface with d xy character also shift, although relatively less. We accounted for this with a 10 meV downwards shift of the d xy orbital energy in the tight-binding model. One way to understand this is by noting that the total volume of the Fermi surface must be approximately conserved in the nematic phase since the distortion of the hole pocket is approximately symmetry-preserving. Therefore a downward shift of the d xy electron band dispersion could simply be a necessary consequence of the upward shift of both the d xz and d yz bands. However alternative interpretations could be possible.
Dirac points: Finally we note that our model indicates that spin-orbit coupling gaps out the d xy − d xz band crossing slightly away from M, which has elsewhere been proposed as a Dirac point. This can also be directly observed in the data, where band discontinuities are observed at the crossing points of d xz and d xy bands, e.g. in Fig. 5. In fact the spin-orbit induced splitting is of comparable magnitude to the Fermi energy and therefore we question whether the Dirac description can be appropriate for these points.