Correlation-driven electronic nematicity in the Dirac semimetal BaNiS2

Significance Topological electronic properties, usually understood in an uncorrelated electron picture, can be guaranteed by crystal symmetries. However, a material’s electrons can spontaneously form a pattern less symmetrical than the underlying crystal due to their mutual correlations, a complex problem underpinning phenomena such as high-temperature superconductivity. In the topological semimetal BaNiS2, we use atomically resolved microscopy to discover a spontaneous rotational symmetry breaking, or “nematic,” electronic phase. We elucidate the unusual mechanism by which this state emerges through electronic correlations and how the topological band crossings of the semimetal persist within this symmetry-broken environment. This helps to elucidate the interplay between topology and electronic correlations, foundational but often discordant concepts for quantum materials.


I. FURTHER RESULTS FROM STM TOPOGRAPHY
As indicated in Fig. 1(b) of the main manuscript, cleavage can occur between adjacent BaS layers revealing (001) facets. The STM topography shown in Fig. 1(g) is repeated in Fig. S1(a) below. Its Fourier transform F[T (r)] is shown in Fig. S1(b) and the primary set of reciprocal lattice peaks represent a periodicity of 4.52Å, close to the reported value for the surface lattice vectors a = b = 4.43Å [1]. Although the Ba and apical S are uppermost at the cleaved surface, because the densityof-states at low energies is dominated by Ni orbitals, the atomic corrugation likely represents the Ni square net. Moreover, because the observed reciprocal lattice vector indicates only one Ni atom per unit cell is observed, the corrugation therefore represents specifically the uppermost Ni A square net.
An additional topography image that includes a step between two atomic terraces is shown in Fig. S1(c), and a height profile, shown in Fig. S1(d), shows an inter-layer separation of 9Å, close to the lattice parameter of c = 8.89Å [1].

II. LANDAU LEVEL SPECTROSCOPY
Under high magnetic fields, in suitably two-or quasitwo-dimensional metals, Landau quantization can be distinguished in conductance spectra measured using STM, in which Landau levels appear as conductance peaks [2][3][4]. Figure S2(a) shows dI dV (E) spectra acquired under various magnetic fields. The magnetic field was applied approximately parallel to the c-axis. (Surface tilt was measured to be less than 1 • .) Each curve represents the average over the same field of view as for Figs. 1(h) and 1(i) in the main text. The magnetic field was increased in increments of 1 T, up to 6 T, and then in increments of 0.5 T, up to the maximum field of H = 12 T.
A subtle change in the spectra appears in the region between E = 0 and 50 meV, and this is clarified upon taking the zero-field conductance curve as a baseline and subtracting it from all other curves, as shown in Fig. S2 (b). Two peaks are seen to become prominent approaching H = 12 T. The peak at lower energy shifts significantly with magnetic field, while the peak at higher energy appears to have only a small shift.
Knowing the energy of the zeroth Landau level of an ideal topological, massless Dirac fermion system should have no field-dependence (other than due to the Zeeman effect), we anticipate that the upper of the two peaks in Fig. S2 (b) can be assigned LL 0 , and the lower of the two can then be assigned as LL −1 . In the more general case in which a Dirac cone may exhibit a finite mass m * , even the energy of the zeroth Landau level can exhibit a small linear relation with magnetic field H. This is quantified within the derivation for the energy of the n th Landau level given by Fu et al. [4]. We follow from this derivation, but neglect the Zeeman term, because no Zeeman splitting is discernible in the dI dV (E) curves despite the two-fold degeneracy of the Dirac points in BaNiS 2 . We write the resulting relation as (1) where E D is the Dirac point energy, is the reduced Planck constant, ω c = |e|µ0H m * me is the cyclotron frequency, e is the electronic charge,v is the average velocity on the contour of the band from which the Landau level forms, and µ 0 is the vacuum permeability. The phase factor γ is given by γ = 1 2 − ΦB 2π where Φ B is the Berry phase accumulated around the Landau orbit. This signifies the topological character of the system and takes the value π if the orbit encloses a non-trivial Dirac node. In the case that γ = 0, the energy of the zeroth Landau level (LL 0 ) remains at E D independent of magnetic field.
Using this framework, and adopting the assumption that γ = 0, we estimate E D and m * by fitting to the data as shown in Fig S2(c). For this procedure, after a twopoint (0.5 meV) Gaussian filter is applied to smooth each curve, the Landau level peak energies are initially guessed  Fig. 1(g). Each curve is offset vertically by an amount proportional to the applied magnetic field. (b) The field-dependence as represented by subtracting the zero-field curve [dashed curve in (a)] from all others. Each curve is offset vertically as in (a). Anticipating the following result, we label the peak near to ED and the first peak below it as LL0 and LL−1, respectively. (c) Results of fitting the two labeled peaks using Eqn. 1, and using field-dependence curves from 6 to 12 T in increments of 0.5 T.
. We also include the constraint that the fitted relations for the two Landau levels must be kept consistent with each other in that they converge upon the same value of energy at H = 0 T, namely E D .
Given the above assignment of Landau level indices to the peaks, and other assumptions, the fitting estimation yields E D ≈ 26.6 meV and m * ≈ 0.6.
To summarize the key conclusions from Landau level spectroscopy results, we have determined that i) the minimum in dI dV (E) at E − E F = 27 meV, seen in Fig. 1(i) of the main text, does correspond to the Dirac point energy, and also that ii) the observation of a Landau level with a weakly linear field-dependence, apparently due to finite mass, is consistent with a Berry phase of π, signifying the topologically non-trivial nature of the Dirac cones in BaNiS 2 .

III. ENERGY-DEPENDENT INTENSITIES OF BRAGG PEAKS IN FIGURE 2
Figure 2 of the main text displays intra-unit-cell sampling of the tunneling conductance dI dV (r, E) in a small field-of-view, in which nematicity manifests as a difference between curves sampled at the two inequivalent bond-centered sublattices of the uppermost Ni square net. An alternative visualization of nematicity is suggested by the asymmetry between Bragg peaks in the Fourier transforms of the L(r) images [Figs. 2(e) and 2(f) of the main text].
In Fig. S3 we shows the energy-resolved intensities of the Bragg peaks in L q (q, E), for comparison with the real-space sampled dI dV (E) and d 2 I dV 2 (E) curves in Figs. 2(i) and 2(j) of the main text. Figure S4 shows measurements across a step-terrace morphology found at a BaNiS 2 surface. Two steps were Each set of images exhibits nematicity with the same orientation. This observation is consistent with two possible scenarios. First, the nematicity may be present in the bulk as well as at the surface, and have a 'ferro'type out-of-plane correlation. Second, the nematicity in the three layers may be uncorrelated, but triggered by a global condition such as a small macroscopic sample strain. The latter scenario allows that nematicity emerges only at the surface, but is absent in the bulk.

IV. NEMATICITY ACROSS A STEP-TERRACE MORPHOLOGY
The density-wave equation calculations presented in the latter part of the main text are based on a purely 2D system which applies to both the surface and the bulk if inter-layer hopping is neglected. The effect of three dimensionality on the nematicity is an interesting issue for future investigation.

V. FITTING OF QUASIPARTICLE INTERFERENCE SIGNALS
Here we describe the procedure for fitting to the linecut along G a taken through the L q (q, E) data, displayed in Fig. 3 of the main text. This leads to the acquisition of solid red and orange lines for q 2 and q 3 shown in Fig. 3(c) of the main text. The same procedure is also applied to the linecut data along G b , in order to make the comparison shown in Fig. 4 of the main text.
A model consisting of the sum of multiple Lorentzian lineshapes and a constant offset, c q , is fitted at each point along the q path. This is written as Here N is the number of Lorentzian lineshapes included in the model, which is automatically determined for each q point using a peak finding algorithm [5] after some smoothing. Figures S5(a) and S5(b) show the raw linecuts through the measured L q (q, E) data, along the G a and G b reciprocal lattice vectors. Figures S5(c) and S5(d) show fitting results as 2D images, for comparison against the measured data. Figures S5(e) and S5(f) show the same fitting results as points with 'error bars' that signify the Lorentzian widths Γ n . The subset of points attributed to scattering within and between the Dirac cones were chosen for subsequent linear fits. The results of these linear fits are displayed as the solid red and orange lines shown in Fig. 3(c) of the main text., and provide quantitative information about the Dirac cones, as discussed in the following section.

VI. QUANTITATIVE DESCRIPTION OF DIRAC CONES DERIVED FROM QPI RESULTS
As described in the main text, the apparent dispersion of a QPI signal has only an indirect relation to the actual dispersion of the related bands in k-space, and for this reason interpretation of QPI observations is generally not straightforward. In the present case the relation is complicated by the fact that the Dirac cones are not centered at q = 0, and that only scattering between Dirac cones is observed, and not scattering within a single one.  Nevertheless below we draw a comparison between the observed and calculated features in q-space, and also a k-space-to-q-space comparison of the calculation results, specifically for the QPI signals q 2 and q 3 shown in Fig. 3 of the main text. From this we discuss whether the Dirac point energy can be inferred from QPI observations, and also derive estimates of the Dirac cone velocities parallel to the ΓM line. First, if the observed scattering branches q 2 and q 3 result from scattering between adjacent Dirac cones as depicted in Fig. 3(d) of the main text, they intersect at the energy where each of the Dirac cones vanishes to a point, so that the shortest and longest possible vectors connecting their contours (namely q 2 and q 3 ) become identical. The energy at which they intersect (E = 101 meV) would then seem to be the Dirac point energy. However, we find that this is not true. The calculated Dirac point energy is in fact about 60 meV [see Fig. 1(d) in the main text], meaning the extrapolated point of intersection gives an over-estimate, perhaps due to a finite curvature as the band nears the Dirac point. This indicates that an ex-trapolation from the branches q 2 and q 3 in the observation (which would yield E Dirac ≈ 47 meV) would similarly over-estimate the actual Dirac point energy. For this purpose we fall back on the minimum in the dI dV (E) curve [ Fig. 1(i) of the main text], and also Landau level spectroscopy above, to determine the Dirac point energy.
As described above, the point of intersection of q 2 and q 3 represents the q-vector connecting two adjacent Dirac points, which we may call q D−D . The assumption that each Dirac point lies on a ΓM line in k-space leads to |q D−D | 2 = 2 |k D | 2 where k D , the wavevector of the Dirac point, is directed from the origin along the ΓM line. This results in |k D | ≈ 0.884 nm −1 .
An approximate value for the Dirac cones' Fermi velocity v F = dEF dk along a particular axis can be inferred from the apparent velocity v q,2,3 of the scattering branches q 2 and q 3 . We start from the assumption that v q for a generic scattering branch always relates to the velocity v k of a chosen point on the band by a proportionality constant that captures geometric details of the scattering vectors allowed within the band's structure, i.e. v q = cv k . (As an example, v q = 1 2 v k in the extremely simple case of scattering across the diameter of an isotropic Dirac cone).
Next, we compare the calculated and observed velocities v obs. Using the same method of fitting as shown in Fig. S4, but for the surface spectral function shown in Fig. 1(d) of the main text, gives velocities for the inner and outer arcs of the Dirac cone, along the ΓM line, of v calc. k,inner = 5.32 × 10 5 m s −1 , and v calc.
The band velocities along the ΓM line are higher than the average Fermi velocity as determined through Landau level spectroscopy above. This is partly because the Dirac cones are elliptical, having the highest velocity along ΓM and the lowest velocity perpendicular to this line. It may also be because the Fermi contour generally has a lower velocity as the Dirac cone has some curvature, becoming lower in velocity near E D . The estimate of velocities based on QPI dispersions relies on energies between −50 and −100 meV, where the cone has steeper and more linear dispersion.

VII. ORIGIN OF ENERGY SAVINGS DUE TO NEMATIC ORDER
Here, we discuss why an energy gain results from the present finite-energy nematic order. Although the Dirac nodes are almost unchanged by the B 2g form factor, the band dispersion along the antinodal direction is significantly modified. First, we discuss the change in the free-particle energy due to the form factor f l (k) = f ×f l (k), wheref l (k) is the normalized form factor (max l,k |f l (k)| = 1). In the hole representation, it is given by where f ( ) is the Fermi distribution function, f n,k is the n-th band dispersion with the form factor. We perform the numerical study at T ∼ 0.01, where Eq. (3) is very close to the free energy because T S is small. In the present model, the relations a < 0 and b > 0 hold, and therefore ∆E f free < 0 at a finite f . This energy gain mainly originates from the electron-like dispersion around the X points (Fig. 5). In fact, the bottom of this dispersion (E 0 ≈ 0.05 eV) is given as In fact, d ≈ (E 0 − E 0 ) −1 due to the second-order perturbation theory, and E 0 ≈ −0.5 eV is the energy of the nearest valence band at the X points. Thus, the negative a in Eq. (3) is robustly obtained in the present band-structure.

VIII. TEMPERATURE DEPENDENCE
A natural next step following on from these results is to experimentally characterize the temperature dependence of the nematic order. The energy splitting of ∼12 meV corresponds to a temperature of ∼140 K, naïvely the minimum possible transition temperature for the nematic phase. Unfortunately, tunneling spectroscopy observations suffer from an energy broadening that rises faster than k B T (by a factor of 3.5), meaning that observations similar to those presented in this work probably become impractical before reaching the temperature required to observe the potential melting of the nematic phase. This means that further investigations using other tools may be necessary to characterize its temperature dependent behavior.

IX. RELATION BETWEEN NEMATICITY IN
BaNiS2 AND MAGNETISM IN BaCoS2 The antiferromagnetic phase of BaCoS 2 can be thought of as an arrangement of collinear spin chains aligned along one of the Ni A -Ni B axes and antiferromagnetically coupled to each other, or equivalently, a spin density wave with a wavevector Q SDW = (π, π) [8]. Likewise, in BaNiS 2 , the magnetic nesting vector is Q = (π, π). In the spin-nematic scenario, the director of nematicity is always aligned with the axis of the expected spin stripe order, which would lead to a director of nematicity oriented parallel to the Ni A -Ni B axis, at 45 • to the pattern observed in Fig. 2 of the main text. In contrast, in the present spin-fluctuation interference scenario, there is no such constraint on the director of nematicity, allowing the configuration as observed in Fig. 2 of the main text, with Q = (0, π) and stripes parallel to Ni A -Ni A .