Observation of an anisotropic Dirac cone reshaping and ferrimagnetic spin polarization in an organic conductor

The Coulomb interaction among massless Dirac fermions in graphene is unscreened around the isotropic Dirac points, causing a logarithmic velocity renormalization and a cone reshaping. In less symmetric Dirac materials possessing anisotropic cones with tilted axes, the Coulomb interaction can provide still more exotic phenomena, which have not been experimentally unveiled yet. Here, using site-selective nuclear magnetic resonance, we find a non-uniform cone reshaping accompanied by a bandwidth reduction and an emergent ferrimagnetism in tilted Dirac cones that appear on the verge of charge ordering in an organic compound. Our theoretical analyses based on the renormalization-group approach and the Hubbard model show that these observations are the direct consequences of the long-range and short-range parts of the Coulomb interaction, respectively. The cone reshaping and the bandwidth renormalization, as well as the magnetic behaviour revealed here, can be ubiquitous and vital for many Dirac materials.

D irac materials 1 are a novel class of solid-state systems, in which the low-energy electronic excitations are described by pseudo-relativistic massless Dirac fermions (DFs) with linear energy dispersion around the Fermi energy E F . Triggered by the studies in two-dimensional (2D) graphene 2 and the surface of three-dimensional topological insulators 3 , extended now to three-dimensional Weyl and Dirac semimetals with strong spin-orbit coupling [4][5][6] , many intriguing properties of DFs have been revealed and have constituted active topics in modern condensed matter physics. The role of Coulomb interaction is one of such issues of particular interest [7][8][9][10][11][12][13][14][15][16][17] . For instance, in charge-neutral 2D massless DF systems, composed of two gapless points at E F , the long-range (LR) part of the Coulomb potential V(q) (q: wave vector) is unscreened owing to the vanishing density of states (DOS) at E F . Consequently, the LR character of the interaction (V(q)p1/|q|) is preserved at low energy, which couples to the fermionic excitations and induces a logarithmic correction to the Fermi velocity v F and associated physical quantities [7][8][9] . The logarithmic velocity renormalization induces a nonlinear reshaping of the cones around each of Dirac points (DPs), as observed in graphene near the charge neutrality point 11-14 . However, graphene is a special case of 2D massless DF systems, in which isotropic Dirac cones with vertical axes have the DPs at particularly symmetric points on the Brillouin zone boundary 2 . Indeed, theoretical studies have revealed that massless DFs possessing anisotropic cones and the DPs at arbitrary k-points emerge more generally in a broad class of materials 1,18-20 . A typical example in 2D is the organic-layered compound a-(BEDT-TTF) 2 I 3 (a-I 3 ) (BEDT-TTF ¼ bis(ethylenedithio)tetrathiafulvalene) (Fig. 1a), which has a pair of Dirac cones occurring at two distinct points (±k D ) in the 2D first Brillouin zone (Fig. 1c) 20- 30 . The electronic structure of a-I 3 , described on the base of molecular orbitals as usual in this type of compounds 31 , is rather involved compared with graphene due to the presence of four sites per unit cell ( Fig. 1b) with anisotropic hopping amplitudes 32, 33 . The system has only the inversion symmetry [32][33][34] , which, in conjunction with the anisotropic hopping, brings about a tilt of Dirac cones and drives the 2D DPs away from high crystallographic symmetry positions 20,21,23,25,35,36 (Fig. 1c). A remarkable feature is that, because of the 3/4-filled nature of the electronic bands [23][24][25][26]33 , the two gapless points are anchored at E F by this band filling in a-I 3 .
Another issue of great physical interest in a-I 3 in terms of the Coulomb interaction problems is that, within the pressuretemperature (P-T) phase diagram ( Fig. 1g) 32,[37][38][39][40][41] , the 2D massless DF phase appears in the vicinity of an insulating phase with charge order, as first pointed out by transport measurements 22 . This contrasts with the case of graphene, in which no phase transitions have been observed at least in the absence of a quantizing magnetic field 10 . The charge-ordered phase in a-I 3 , which is induced by the strong short-range (SR) electron correlations in this 3/4-filled system 40,41 , is suppressed when applying a P above a critical value of P C E1.2 GPa (Fig. 1g) and turns into the 2D massless DF phase 22,39 . Once the high-P phase is reached, the Dirac cones become stable against further pressurization; in fact, the gapless point is fixed at E F on varying the hopping integrals in a finite range by virtue of the 3/4-filled nature of the electronic bands, as revealed by band-structure calculations 23,25,26,33,[42][43][44][45] . The presence of such a phase transition in this system potentially offers the possibility to test the impact of the SR electron correlations on the behaviours of 2D massless DFs. Moreover, the tilt of anisotropic Dirac cones 36 coupled with the SR and LR parts of the Coulomb interaction opens new possibilities in the physics of 2D massless DFs. For instance, it is predicted to bring about a non-uniform reshaping of titled cones 46 , novel non-Fermi liquid behaviours near the quantum critical point 16,17 , where two DPs merge 47,48 , enhanced shot noise for quantum transport 49 and anomalous charge/spin textures inside the unit cell 48, 50 . Studying the electronic structures and the role of the Coulomb interaction in pressurized a-I 3 , which remains unclear up to date, is thus of primary importance to understand the various effects of the Coulomb interaction in 2D massless DFs. In this article, we focus on the 2D massless DF phase in a-I 3 emerging under a hydrostatic pressure (P4P C ) and present experimental evidence for interaction effects of massless DFs. Employing site-selective nuclear magnetic resonance (NMR), we uncover three distinct interaction phenomena induced by the electron-electron Coulomb interaction. First, NMR-shift measurements in conjunction with renormalization-group (RG) analyses reveal a T-driven cone reshaping around each of the DPs due to the LR part of the Coulomb interaction. Because of this reshaping, tilted cones become effectively isotropic at low energies. Second, quantitative RG analyses establish that the best fit to the data inevitably requires a strong bandwidth reduction inherent to the SR electron correlation, as often discussed in strongly correlated materials. Finally, an anomalous ferrimagnetic spin polarization is observed, which is accounted for by the onsite Coulomb repulsion, as revealed by a simulation based on the Hubbard model presented here. These experimental and theoretical investigations demonstrate that a-I 3 under P is an intrinsically interacting 2D massless DF system, in which both the LR and SR parts of the Coulomb interaction strongly influence the electronic behaviours.

Results
Basic principles to probe tilted Dirac cones. Our strategy to investigate tilted Dirac cones in a-I 3 is as follows. The crystal structure of a-I 3 has a 2D unit cell with four molecular sites (dubbed sites A, A' ( ¼ A), B and C), each of which constitutes a sublattice in the crystalline ab-plane (Fig. 1a,b). The four molecular orbitals on these sites form a pair of tilted Dirac cones near E F (Fig. 1c). Around the gapless point at E F , a very unique situation is realized where the Bloch state has different weights in the amplitudes of the four molecular orbitals. The band-structure calculation 25 revealed that these weights, dubbed site-spectral weights hereafter, show anisotropic k dependence around each of 2D DPs with a clear contrast between non-equivalent sites. The corresponding site-spectral weight for the sublattice j ¼ A ( ¼ A'), B and C around the DP at k D , n z j q ð Þ (equation (11)), is shown in Fig. 1d Fig. 1c). (For details, see Methods.) Notably, the anisotropy of the site-spectral weights makes a particular distinction between the site B and C. Namely, the Bloch electrons with large v F (in the steep slope of the tilted cones) have a large weight on the B-site wavefunction, whereas the Bloch states with small v F (on the opposite side of the cones in the gentle slope) have a large weight on the C-site wavefunction (Fig. 1e,f and Supplementary Fig. 1a,c). Thus, if one probes the local electronic states on sites B and C separately by means of a site-selective measurement, it is possible to reveal the electronic nature of the Bloch states in the steep part and the gentle part of the tilted Dirac cones individually. Taking advantage of this feature, we performed a site-selective NMR in this compound to separately elucidate the electronic states in the two slopes of the Dirac cones. Specifically, the Knight shift, derived from the NMR line shift measured at a temperature T, is converted into the local electron-spin susceptibility on the site j, w j s T ð Þ, which is given by a thermal average of the site-spectral weight n z j q ð Þ around E F summed over all q for both electrons (z ¼ þ ) and holes (z ¼ À ). Hence, the site-selective NMR in a-I 3 works as an effectively q-resolved probe of 2D DFs thermally excited around each of DPs. Indeed, the electronic excitations in the steep and gentle parts of the Dirac cones can be almost independently probed by w B s and w C s , respectively, as we will demonstrate below.
NMR observation of tilted Dirac cones. To address the electron interaction issues of 2D massless DFs, we have carried out 13 C-NMR measurements in a-I 3 at P ¼ 2.3 GPa (4P C ; see Fig. 1g) on the four molecular sites in the unit cell, j ¼ A, A' ( ¼ A), B and C. The two 13 C nuclei (spin I ¼ 1/2) introduced at the centre of BEDT-TTF (ET) molecules (inset of Fig. 1a) are used for 13 C NMR, which are known as a sensitive probe of electronic states at E F in this class of compounds 51 . Figure 2a shows the typical NMR spectra observed at a magnetic field of  H ¼ 6 T, applied in the crystalline ab-plane. Eight lines are observed in the spectra, associated with the four molecular sites in the unit cell. By changing the field orientation in the crystalline ab-plane and examining the angular dependence of the line positions, the eight lines are to be assigned to two doublets from sites B and C and one quartet from the site A ( ¼ A') 37 . (The site A and A' are not distinguished hereafter). Figure 2b shows the typical angular dependence of the NMR total shift for each molecular site, defined as the centre-of-mass position of the doublet or the quartet. The j-site total shift for a given T and a field-angle c (measured from the crystalline a axis; Fig. 2c) is expressed as S j T; c ð Þ¼ A j c ð Þw j s T ð Þ þs j c ð Þ, where the first term is the conduction-electron term (Knight shift) and the second term stands for the core-electron contribution (chemical shift). Note that the so-called hyperfine coupling constant, A j c ð Þ, and s j (c) are strongly c-dependent (with little T and P dependence), while  , with alternatingly stacked BEDT-TTF (ET) and triiodide (I 3 ) layers. Balls-and-stick diagram represents the molecular structures, where sticks indicate bonds; red, blue and green balls stand for carbon atoms; grey balls are sulphur atoms; and big purple balls indicate iodine atoms. One electron per two ETs is donated to a I 3 molecule due to charge transfer, constituting a quasi-2D (hole) conducting system in (ET) 2 þ layers and non-magnetic insulting I 3 À layers. Inset of (a): The structure of a ET molecule. 13 C isotopes, substituted for 13 C-NMR, are indicated by arrows. (b) 2D unit cell with four ET sites in the ab-plane, distinguished as the site j ¼ A ( ¼ A') (blue), B (green) and C (red). Same colour correspondence as in a. The dashed rectangle indicates the unit cell, and crosses stand for the inversion centre. (c) Electronic band structure of a-I 3 at high pressures (P4P C ; see g) 25 . As the band is 3/4-filled owing to the charge transfer, the Fermi energy E F ( ¼ 0) is present between the first (z ¼ þ ) and second (z ¼ À ) bands from the top. Gapless points appear at E F and locate at ± k D in the first Brillouin zone, around which a pair of tilted Dirac cones are visible (indicated by circles). Wave vectors (k x and k y ) are given in the r.l.u.
(d-f) The calculated titled Dirac cone around E F ( ¼ 0) plotted as a function of the wave vector q ¼ k À k D ¼ (q x , q y ). The label bar stands for the size of the site-spectral weight, n z j q ð Þ, for the and C (f) 25 (equation (11) in Methods). The magnitude is normalized to the maximum value of n z j q ð Þ for j ¼ B. (g) The pressure-temperature (P-T) phase diagram. The charge ordering temperature T CO (squares) is determined from the electrical resistance (R) measurements 39 , defined as the maximum of À d(ln R)/T. Error bars follow those provided in ref. 39. At T ¼ 0, the phase boundary locates at P C E1.2 GPa as determined from a linear extrapolation (solid line), corresponding to dT CO /dPE-110 K GPa À 1 . r.l.u., reciprocal lattice unit.
the susceptibility w j s T ð Þ is isotropic in this compound 52 . Because the electronic excitations around the gapless point (at E F ) are vanishingly small at low temperatures, the total shift S j at the lowest T ( ¼ 3 K in the present experiment) is expected to provide the chemical shift term s j . Thus, subtracting s j (c) from S j T; c ð Þ and employing the value of the hyperfine coupling constant A j c ð Þ reported at ambient pressure 37 , the total shift S j T; c ð Þ is converted into the local electron-spin susceptibility w j s T ð Þ (for details, see Methods). Figure 2d shows the temperature dependence of w j s T ð Þ at the site j ¼ A, B and C, which arises from the inter-band thermal (electron-hole) excitations across the gapless point at E F and is proportional to the k B T average of the j-site DOS around E F : w exhibits T-independent features down to TE200 K, followed by a rapid decrease with a clear difference in the size of the susceptibility between non-equivalent sites, w C s 4w A s 4w B s (see also Supplementary Fig. 2), and finally becomes vanishingly small at all sites. The observation, in particular the crossover from w j s Bconst. to w j s B0 on cooling, indicates that there is an energy-independent large D j (E) at high energies (inset I of Fig. 2e) and a vanishingly small D j (E) at low energies around E F (inset II of Fig. 2e). This is consistent with the band-structure calculations 24,25,50 , where a flat DOS is predicted at high energies above the van-Hove singularity (locating at about 12 meV off E F ) and linear energy dependence is suggested around the band-crossing point at E F . Note that D j (E) around the gapless point (E ¼ E F ¼ 0) is given by a q summation of the site-spectral weight n z j q ð Þ ( Fig. 1d-f) at a given energy E for the band z (equation (14)). Then, the anisotropic q-dependence of n z j q ð Þ around the DPs of tilted Dirac cones (Supplementary Fig. 1) is expected to bring about D C (E)4D A (E)4D B (E) at low energy. This is in excellent agreement with the observed relation w C s 4w A s 4w B s below TE200 K ( Fig. 2d) and is the direct consequence of the fact that the site B and C selectively probes the steep and gentle slopes of titled cones, respectively, consistent with the prediction of the effective tight-binding (TB) model given in ref. 25. All these observations demonstrate the existence of tilted Dirac cones with E F located around the gapless point.
Fermi velocity renormalization. However, the strong nonlinear temperature dependence in w j s T ð Þ below TE150 K does not comply with the expectation of TB calculations, which leads to w j s pT at low temperatures 25,50 . To better visualize this point, we plot the first derivative of the susceptibility (dw j s /dT), as shown in Fig. 2e. With decreasing T, dw j s /dT exhibits a peak at T B flex E120 K for the site B, and at T A;C flex E50-60 K for the site A and C, and then drops continuously to zero at all sites towards low temperatures. These features are in striking contrast to the TB calculation 25 , where dw j s /dT increases on cooling but saturates at low T (Supplementary Fig. 3d-f). Indeed, w j s T ð Þ varies almost quadratic in T at all sites below the inflection point (T j flex ), which suggests a nonlinear suppression of D j (E) around E F ( ¼ 0) in an energy range of DE j % k B T j flex . As the total DOS, D(E), is proportional to the inverse square of v F in 2D massless DF systems for the non-interacting case D E ð Þ¼ E j j=p' 2 v 2 F (refs 2,20), a suppression of DOS in turn corresponds to an enhancement of v F . Thus, the observed peak structure in dw j s /dT strongly indicates that a T-driven renormalization of v F grows below TET j flex . The most probable origin of this effect is the LR part of the Coulomb interaction between electrons, which is unscreened at E F in charge-neutral massless DF systems and is known to cause a logarithmic correction to v F either driven by tuning carrier densities 1,7-13 or temperatures 9, 46 . We recall, however, the value of T j flex is twice higher for the site B (T B flex E120 K) than for the site A and C (T A;C flex E50-60 K). At first glance, this may add an extra complication to the data interpretation but in fact turns out to be a direct consequence of the anisotropy of the site-spectral weight n z j q ð Þ ( Fig. 1d-f) and the tilt of Dirac cones, as we shall see below.
Renormalization-group analyses. To further understand the nonlinear temperature dependence of w j s T ð Þ at each site, we have examined the self-energy correction effect due to the LR Coulomb interaction. For this, we employed a RG approach based on an effective Hamiltonian near the gapless point 20,25 , whose energy-momentum dispersion is given by where w 0 ¼ (w 0x , w 0y ) and v ¼ (v x , v y ) are velocities reflecting the tilt and anisotropy of the cone, respectively (for details, see Methods). At the one-loop level, the self-energy correction leads to a renormalization of v but does not affect w 0 (Supplementary Fig. 4a) 46 . The RG flow of v is expressed as where N ¼ 4 is the number of fermion species, corresponding to two DPs and two spin projections, is a momentum cutoff of the size of the inverse lattice constant 33 and is circular around the DP, Þ=g j and e is the dielectric constant. (Note that equation (2) is obtained in the leading order in 1/N assuming N4 41, which is valid both for the weak and strong Coulomb interaction.) Assuming the four velocities given by the effective TB calculation (w TB   0 and v TB ) 25 as initial velocities at q ¼ L, we have calculated the RG correction effects on w Here, we note that the temperature is used as an explicit scale parameter in the calculation of w j s T ð Þ that determines the RG flow. To get a reasonable agreement between the calculation and experiment, a phenomenological parameter u is introduced to adjust the velocities at q ¼ L such that w 0 0 ¼uw TB 0 and v 0 ¼ uv TB . The two parameters in the calculation, (u, e), are then optimized from a least-square fit to the experimental susceptibilities. Good agreements are obtained in the fit especially at the site A and C . The calculation demonstrates that the nonlinear T dependence of w j s T ð Þ below TET j flex can be properly ascribed to the logarithmic renormalization of v F . In Fig. 4a-c, the calculated shape of D j (E) is shown around the gapless point at E F . A strong suppression from the E-linear DOS is seen at low energies due to the renormalization. Figure 4d-f, depicts the corresponding energy spectrum around the DP (at k D ), where the colours indicate the magnitude of the site-spectral weight, n z j q ð Þ, in Fig. 1d-f, respectively. A nonlinear reshaping of the tilted cone induced by the renormalization is clearly visible around the gapless point. It should be stressed that a good agreement is NATURE COMMUNICATIONS | DOI: 10. accomplished only when a small value of u (E0.35) is used. The fact that we have uo1 indicates a reduction of the initial velocities or of the hopping amplitudes between the adjacent lattice sites. In conventional strongly correlated materials, the SR part of the Coulomb interaction is well known to induce this sort of hopping (or bandwidth) renormalization due to the frequency dependence of the self-energy 53 . We believe that the observed u-reduction effect in the initial velocities occurs because of this self-energy correction due to the SR Coulomb interaction, which is not considered implicitly in the original RG calculation.
Remarkably, the calculation well reproduces the observed difference in the thermal energy scale of the T-driven renormalization in w  Supplementary Fig. 1a,b). Namely, the energy cutoff is large for the large-v F DFs, dominantly probed by the site B, while it is small for the small-v F DFs, having a large weight on the site C (and A). Hence, the renormalization starts from a higher T in w B s (T) than in w C s (T) (and w A s (T)), producing the observed energy scale difference DE B 4E A,C , consistent with the previous RG calculation of Isobe et al. 46 . It is also worth mentioning that tilted cones become more isotropic at lower energies near E F because of the non-uniform velocity renormalization around each DP, as reflected in Fig. 4d-f. This is because the anisotropic term in the Hamiltonian is small 25), and we have |w 0 |o o|v| around the DP due to the RG flow ( Supplementary Fig. 4a), leading the tilting term (w 0 ) to be effectively negligible near E F .
From all these, it is concluded that our RG analyses appropriately capture many of the essential parts of the experimental results. They constitute experimental evidence for the bandwidth renormalization (the u-reduction effect) due to the SR repulsion between electrons as well as the T-driven logarithmic renormalization of v F by the LR part of the Coulomb interaction. Nevertheless, we note that, at low temperatures, the agreement is less satisfactory for the site B compared with the other sites ( Fig. 3c), suggesting the presence of another correlation effect. Indeed, we will clarify this point by a lattice-model simulation, as described below.
Emergent ferrimagnetic spin polarization. The temperature dependence of w B s (T) in the experiment is appreciably stronger and more complex than what is predicted by the RG calculation ( Fig. 3c). Indeed, the experimental w B s (T) exhibits an anomalous sign change at TE60 K and an upturn with a negative slope below TE40 K (inset of Fig. 3c), while the RG calculation shows monotonic temperature dependence. The observation of w B s o0 is in sharp contrast to w A s 40 and w C s 40 in the experiment (Fig. 5a), indicating an emergent ferrimagnetic spin polarization in which the local magnetic field points antiparallel to the applied field at the site B while it is parallel to the field at all other sites (Fig. 6). To further understand this sublattice-scale magnetism, we have investigated the Hubbard model with an onsite repulsive (Hubbard) interaction, U, at a mean-field level within the random phase approximation (RPA). For the RPA calculation, we have considered both the inter-band and intra-band contributions to the spin susceptibility with a wave vector Q ¼ 0 (for details, see Methods). Figure 5b presents the calculated temperature dependence of the RPA spin susceptibility at the site B. Using U ¼ 0.14 eV, the RPA calculation (in particular the inter-band term) clearly reproduces the observed negative spin susceptibility for site B (w B s o0) at similar temperatures. Moreover, the negative susceptibility appears only at site B in the calculation (Supplementary Fig. 3a-c) in good agreement with the experiment ( Fig. 5a and Supplementary Fig. 2). These facts show that (symbols). Thin arrows indicate the inflection point, T j flex . As the non-interacting reference in the RG theory, we assumed the effective TB model of ref. 25. The susceptibilities are calculated as discussed in the text by employing the momentum cutoff L ¼ 0.667 Å À 1 and the optimum parameters (u, e) ¼ (0.35, 1), determined from the fit to the data (Methods). The dotted lines (w j s pT) indicates the expected T dependence for gapless excitations around the DP in the absence of the RG correction 68 and B (f) in the unit cell plotted as a function of the wave vector q ¼ k À k D ¼ (q x , q y ) (in Å À 1 ). The outer grey cone stands for the linear spectrum in the absence of the LR Coulomb interaction, where the bandwidth reduction effect is considered (corresponding to the E-linear DOS in a-c Generally speaking, RPA tends to overestimate the effect of correlations, as it does not consider the self-energy correction to the energy bands 54 . In fact, the onsite interaction we use, U ¼ 0.14 eV, is chosen smaller than what would be typically used (UE0.4 eV) (refs 25,26). This discrepancy brings another support that the self-energy correction due to the SR interaction is of great importance in this system, consistent with the bandwidth reduction effect we discussed in the RG calculation. Finally, we note that RPA is unable to reproduce the observed nonlinear T dependence of w j s T ð Þ at low temperatures. This is because the self-energy correction due to the LR Coulomb interaction, which is the main origin of the nonlinear w j s T ð Þ, is not taken into account in RPA (for details, see Methods and Supplementary Fig. 3d-f).

Discussion
So far, we have demonstrated three distinct Coulomb-interaction phenomena in the 2D massless DF phase of a-I 3 , which develop systematically at different temperature scales (or energy scales of thermal excitations), as summarised in Table 1. A bandwidth renormalization (or the u-reduction of v F ) occurs due to the SR Coulomb interaction which appears to exist from room temperature down to lowest T. At temperatures TrT j flex (or in the corresponding energy range around the gapless point at E F ), a T-driven logarithmic renormalization of v F and the resultant non-uniform reshaping of the Dirac cones appear, due to the LR part of the Coulomb interaction. With further decreasing T, a ferrimagnetic spin polarization shows up in the unit cell because of the onsite Coulomb repulsion between electrons.
First, we mention that the observed non-uniform reshaping of titled Dirac cones in a-I 3 should affect other physical observables at low temperature or at low magnetic field. As the cones become effectively isotropic around each of DPs ( Fig. 4d-f), the shape of the cross-section of the cone is different at high energy and close to the gapless point, which should cause a T-dependent anisotropy of the in-plane electrical conductivity. Another experiment that would be able to see the reshaping is infrared spectroscopic measurements, known as a powerful tool to reveal the Landau level (LL) structure in graphene 55 . In a perpendicular magnetic field (H > ) normal to the 2D plane, the massless DFs in a-I 3 exhibit the LL spectrum, E z;n ¼z' l À 1 refs 56,57), where z ¼ ± distinguishes the electron and hole bands crossing at the gapless point, n ( ¼ 0, ± 1, ± 2, ?) is the LL index, l B ¼ (' /eH > ) 1/2 is the magnetic length and is a tilting parameter. As E F locates at the gapless point in a-I 3 and the n ¼ 0 LL is half-filled, one may be able to detect, for instance, the dipole transition n ¼ 0-n ¼ 1 in the absorption line at the energy logarithmic correction may appear in DE 10 at low field; v x (Ev y ) will increase by a factor of four by changing ln(L/q) from 1 to 5 ( Supplementary Fig. 4a) and (1 À O 2 ) 3/4 increases by a factor of two (from B0.5 to B1.0). One would expect to see this change in DE 10 at low temperatures where the LL broadening, mainly caused by the thermal scattering of carriers in the n ¼ 0 LL (refs 27,28), becomes sufficiently small. Theoretically, the renormalization of the coupling constant of the Coulomb interaction (the RG flow of v F ) makes the system flow to a weak coupling regime at low energies (or at low T in the case of a T-driven RG flow as in our experiment) 10 . As organic compounds are very clean and are little influenced by impurities, this consequence of the RG flow implies that low-T electron correlation effects, if any, would be induced by the SR part of the Coulomb interaction, as usually the case in conventional strongly correlated materials. In fact, the insulating phase possessing charge order 22,26,31-41 emerges next to the 2D massless DF phase on the P-T phase diagram in a-I 3 (Fig. 1g), suggesting the vital role of the SR Coulomb interactions in the massless DF phase of this 3/4-filled system, in line with recent mean-field calculations 48 . In the half-filled system graphene, strong electron correlations are predicted to stabilize Mott insulators and charge density waves 10 . However, the typical experimental conditions seem to locate away from these situations 1,10,58 , and no phase transition has been reported yet. Under a strong magnetic field, on the other hand, a gapped liquid phase is observed with a quantized Hall conductivity at fractional filling factors, stabilized by the SR part of the Coulomb interaction 10 . Similar physics may occur in thin films of a-I 3 as well, where the integer quantum Hall effects have recently been observed 59 .
To conclude, our NMR measurements combined with theoretical calculations have demonstrated three T-dependent Coulomb interaction effects of 2D massless DFs (Table 1) in pressurized a-I 3 (P4P C ), having tilted Dirac cones and gapless points fixed at E F . We found that the LR part of the Coulomb interaction, which is unscreened around the gapless DPs, causes a T-driven renormalization of the velocity and induces a nonuniform reshaping of tilted Dirac cones. Quantitative analyses of the cone reshaping based on the RG approach further necessitates a large bandwidth reduction due to the SR electron correlation. Moreover, we showed that the onsite Coulomb repulsion gives rise to a ferrimagnetic spin polarization as unveiled by the numerical calculations using the Hubbard model. These findings can be distinguished from the case of weakly interacting 2D massless DFs in graphene with vertical Dirac cones and are consistent with the emergent correlated phase on the verge of the massless DF phase in the P-T phase diagram. Continuing this study to the vicinity of the phase transition at P C (E1.2 GPa) would be of particular interest, which may connect the physics of the massless DFs and conventional strongly correlated materials.
To perform 13 C-NMR measurements, the central carbon atoms of BEDT-TTF molecules connected by a double bond were 99% enriched by carbon-13 ( 13 C) isotopes (inset of Fig. 1a) with a nuclear spin I ¼ 1/2.
Pressurization scheme. A hydrostatic pressure of P ¼ 2.3 GPa was applied to the sample using a BeCu/NiCrAl clamp-type pressure cell, with the Daphne 7373 oil as a pressure medium. At this pressure, the oil locates close to the liquid-solid phase transition point at room temperature 61, 62 . To avoid applying uniaxial strains to the sample, the cell was kept at a sufficiently high temperature (TE50°C) during the pressurization. With decreasing temperature from T ¼ 300-3 K, a pressure reduction of DPB0.1 GPa occurs inside the cell 61, 62 . The inner pressure at the lowest T we measured (3 K) is, however, substantially higher than the transition pressure to the charge-ordered phase at TE0 (P C E1.2 GPa; see Fig. 1g) 22,30,39-41,63,64 , indicating that the P-reduction will not change the physics. Hence, we neglect this effect throughout the paper. 13 C-NMR measurements. 13 C-NMR measurements were performed in a-I 3 in a magnetic field H of 6 T applied parallel to the crystalline ab-plane (Fig. 1a). To get NMR signals, the standard spin-echo techniques were used with a commercially available homodyne spectrometer. Spin-echo signals were recorded at a fixed radio frequency after the conventional spin-echo pulse sequence of t p/2 À t À t p (with t ¼ 5-25 ms and t p/2 À t p ¼ 0.6-1.2 ms) and were converted into the NMR spectra via Fourier transformation. The resonance frequency of the natural abundance 13 C nuclei in TMS (tetramethylsilane (CH 3 ) 4 Si) was used as the origin of the NMR shift. We note that the present work is targeted at lower T and higher P compared with the earlier NMR studies 65,66 , which is more suitable for exploring the nature the low-T 2D massless DF phase emerging at P4P C on the phase diagram ( Fig. 1g) 22,[27][28][29][30]39 .
Line assignments of the NMR spectra. The details of line assignments for the 13 C-NMR spectra in this compound are given elsewhere 37, 65 . Here, we shall describe only the essence. The 13 C-NMR spectra of a-I 3 show large temperature dependence ( Fig. 2a and Supplementary Fig. 9) and field-angle (c) dependence ( Supplementary Fig. 10). They have eight 13 C lines that consist of two doublets (from the B and the C molecules) and one quartet (from the A ( ¼ A') molecule) (Fig. 2a). The doublet and the quartet are caused by the nuclear dipole-dipole interaction in a ET molecule (between the two 13 C nuclear magnetic moments around the molecular centre; see the inset of Fig. 1a) 37,51, 65 . The NMR total shift S j for a particular site j is determined from the centre-of-mass position of the corresponding 13 C lines, which is expressed by a sum of the Knight shift K j (¼ A j w j s ) and the chemical shift s j terms, S j T; c ð Þ¼K j T; c ð Þþs j c ð Þ, as mentioned in the text. (The orbital van-Vleck contribution is negligible in ET-based salts because of the low lattice symmetry 51 .) Both K j T; c ð Þand s j (c) are highly anisotropic in this system 37,65-67 , causing clear c-dependence of the spectra. The anisotropy of the shift S j T; c ð Þis, however, largely different between the molecule A(A') and molecules B and C, reflecting the fishbone-like arrangement of ET molecules in the crystalline ab-plane (Fig. 1b). Thus, by rotating the magnetic field H in this plane, one can assign the 13 C lines to the different sublattices 37 , with the aid of the X-ray diffraction data under pressure 33 (Supplementary Fig. 10).
Conversion of the NMR shift to the spin susceptibility. In a-I 3 , the T-dependence curve of the j-site total shift S j T; c ð Þshows a largely distinct feature for different field orientations. For instance, the T dependence of S A A 0 ð Þ T; c ð Þis large at c ¼ 60°showing a prominent decrease with decreasing T ( Supplementary  Fig. 11a), whereas it is small at c ¼ 120 o with an increase on cooling ( Supplementary Fig. 11b). This difference can be accounted for by the anisotropy of the T-independent hyperfine coupling constant A j c ð Þ, where A j c ð Þ can be either positive or negative depending on the value of c, as we will describe below. First, we note that the chemical shift s j (c) is T-independent 51 as well as little affected by P (see Supplementary Discussion). Thus, the observed T dependence of the total shift S j can be ascribed to the T-variation of the Knight shift K j ¼ A j w j s. Here, A j is the hyperfine-coupling constant averaged for the two central 13 C nuclei in the molecule j (Fig. 1a), which is c-dependent reflecting the anisotropy of the coupling tensor 33,37,66, 67 . The principal values of the tensor are, however, weakly affected by the change of T and P ( Supplementary Figs 12 and 13). Moreover, the spin susceptibility w j s is expected to be isotropic in this compound 52 . These points clearly indicate that K j can practically be expressed as K j (T, c)¼ A j c ð Þw j s T ð Þ (for details, see Supplementary Discussion).
Notice that there are no excitations around the gapless DP in the ground state of massless DF systems. This means, at the lowest T, the Knight shift K j (T, c) is expected to become vanishingly small 68 , and the total shift S j T; c ð Þresumes to the T-independent chemical shift s j (c) (Supplementary Fig. 12). For the site B, there is a negative slope in the T dependence of the total shift S j T; c ð Þbelow TE40 K, showing a small increase of the shift (of B5 p.p.m.) towards lower T (inset of Supplementary Fig. 11b). This is to be associated to the ferrimagnetic spin polarization, which causes a local magnetic field that points opposite to the external field only at the site B (Fig. 6) 42 . The effect is, however, negligible at the lowest T since the thermally excited polarization vanishes as T-0 ( Fig. 5b and Supplementary Fig. 3c). Hence, we fitted to the angular dependence of S j T; c ð Þfor all sites at the lowest measured temperature (T ¼ 3 K) and assumed this fitted curve  Fig. 2e). T j flex stands for the peak temperature in the first derivative of the electron spin susceptibility w j s (Fig. 2e) at the non-equivalent site j ¼ A ( ¼ A'), B and C in the unit cell (Fig. 1b), o indicates the frequency and u represents the phenomenological parameter introduced in the RG fit analyses (see the text). NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12666 ARTICLE to be the chemical shift value at each angle c. The total shift S j T; c ð Þfor a particular site is then converted to K j (T, c) by subtracting this s j (c). The subtraction is done at a field orientation where the total shift becomes close to the maximum in order to minimize the ambiguity of the chemical shift, namely, c ¼ 60°for the site A(A') and c ¼ 120°for the site B and C ( Supplementary Figs 11  and 12). The value of the hyperfine coupling constant A j c ð Þ is calculated for these angles by employing the coupling tensors given at ambient pressure 37 by means of X-ray diffraction data reported at high pressure 33 Þ¼9.9 in kOe m B À 1 . In terms of these coupling constants, the Knight shift K j (T,c) is eventually converted to the local electron-spin susceptibility w j s T ð Þ (Fig. 2d).

. This yields
Effective tight-binding model for the tilted Dirac cone. In order to construct reasonable arguments for the fitting analyses to the observed w j s T ð Þ (Fig. 2d), we have introduced a four-band band-structure calculation of ref. 25. The low-energy effective model based on this calculation shall be used as a non-interacting reference to the data analyses in particular in the RG calculation. (For details, see Methods: RG calculations). The rationales behind this choice are described here.
In a-I 3 , it has been shown that the gapless DPs at ±k D (Fig. 1c) are fixed at E F for a certain parameter region in the four-band TB parameter space with and without finite site potentials 23,33,[42][43][44][45] . However, it is difficult to use bare TB parameters as adjustable variables in the fitting analyses of the low-T part of w j s T ð Þ (in Fig. 2d). This is because the TB calculations lead to linear dispersion around the DP 21,23 , which causes w  (Fig. 2d,e) and a negative w B s (T) at low T (rT B flex /2E60 K; see the inset of Figs 3c and 5a). A model calculation based on simple linear dispersion can hardly account for these features. Thus, instead of fitting the data with TB parameters, we will use them as a minimal non-interacting reference and perform more sophisticated RG calculations based on a continuum model derived from that reference.
For the minimal model in our study, we use the band-structure calculation of ref. 25, which is practically based on a non-interacting TB model with adjustable site-dependent potentials, associated to the four molecular sites in the unit cell, j ¼ A(1), A'(2), B(3) and C(4) (Fig. 1b). (This model shall be dubbed the effective TB model throughout this study.) Strictly speaking, this model takes into account the electron-electron Coulomb interaction up to the nearest-neighbour terms and, at first glance, appears not to be a non-interacting model. However, as one works within a mean-field framework, the interaction merely ends up in additional site potentials in the expression of the Hamiltonian 40, 41 . In this sense, ref. 25 can be also considered as a TB model with adjustable site-dependent potentials from a practical point of view. Importantly, it is well known that the presence of the gapless point at E F is unaffected by modulations of this kind of site potentials within a range in this compound 19,33,42,43 . The chosen values of the site potentials in ref. 25, which are given by with the kinetic terms The hopping integrals are better determined by ab initio calculations such that the resultant electronic bands become compatible with experimental observations in this system. For this, we employed the hopping integrals reported by the firstprinciple density-functional calculation at T ¼ 8 K (ref. 24), as in ref. 25, which are given for the nearest neighbours by (in the unit of eV) and for the next nearest neighbours as (For the definition of the integrals, see Supplementary Fig. 8.) The largest integrals, t LT b1 and t LT b2 , are known to vary about 15% by raising T from 8 to 300 K (ref. 24), though the variation is less than a few per cent below 100 K. As our fitting analyses primarily focus on this low-T region, it is reasonable to omit the T-dependence and keep using the hopping integrals estimated at T ¼ 8 K, {t LT p ; p ¼ a1 À a4'}, at all T as is done in ref. 25. By diagonalizing the Hamiltonian (equation (4)) in conjunction with the hopping integrals (equations (5)-(7)) and the site potentials (equation (2)), one obtains the four energy bands with tilted Dirac cones at E F , as shown in Fig. 1c.
We note that it is very important to use these hopping integrals in equations (3) and (4) 70.) From all these, we used the effective TB model in ref. 25 as our non-interacting reference to the RG analyses. It should be stressed that we do not mean to incorporate interaction effects at this level, and indeed the model we assumed is a purely non-interacting TB model with acceptable site-dependent potentials. Note that these values of site potentials are realistic because they lead to a site-dependent charge differentiation which is compatible with the observed X-ray and Raman scattering results in the conducting phase; see refs 25,32,34. Generalized Weyl Hamiltonian and site-spectral weight. Around the band-crossing DPs, where the Fermi energy E F (put as E ¼ 0 hereafter) is fixed in a-I 3 due to the stoichiometry, the low-energy continuum model is shown to be given by the generalized Weyl Hamiltonian 20,23,25,35 which describes the electronic states in the vicinity of one of the DPs. Here, w 0 ¼ (w 0x, w 0y ) and v ¼ (v x , v y ) are effective velocities describing the tilt and the anisotropy of the Dirac cone, respectively; s 0 is the 2 Â 2 unit matrix; (s x , s y ) are the Pauli matrices; and q ¼ (q x , q y ) is the 2D wave vector measured from the DP at k D . Note that the twofold valley degeneracy associated to the two DPs at ± k D (Fig. 1c) will be omitted for simplicity, and we will hereafter focus on a single valley (at k D ). (When the valley degeneracy is to be included in some of the expressions, we will specifically mention it.) The Hamiltonian (equation (8)) is defined in a space spanned by the Luttinger-Kohn bases 71 : F LK1 D and F LK2 D . These bases are the two degenerate Bloch states at k D , which are given by a (normalized) superposition of the highest occupied molecular orbital h j on each of the four different BEDT-TTFs (ETs) in the unit cell 36 where j ¼ A(1), A'(2), B(3) and C(4) represents the different sites (see Fig. 1b). Diagonalization of equation (8) yields the eigenvalue E z (q) (equation (1)) in terms of the two bands (z ¼ ± ) and the eigenstates (Goerbig, M.O., private communication.) where tan j q ¼ v y q y /v x q x . We note that the two states F LK1;LK2 D in equation (10) have an equal weight for any value of q, as in the two-band model of the graphene Dirac cone, whereas the four states h 1;2;3;4 in equssation 10 have not necessarily the same weight 25,72 . To see this, we define a (normalized) q-dependent site-spectral weight by taking a projection of c z where f Using these velocities, the phase j q in equation (11) can be approximated as j q % arctan q y =q x À Á j ð Þ, where j is the angle between q and the k x -axis. It is shown from equation (11) This causes an oscillation of equation (11) as a function of j around the DP, with a large amplitude and an opposite phase on the j ¼ B-site and the j ¼ C-site, whereas the j dependence is small on the j ¼ A(A') site ( Supplementary Fig. 1c). As the cone is tilted in the k x -direction (inset of Supplementary Fig. 1a), the Fermi velocity becomes highly anisotropic around the cone 20,25,50 ( Supplementary Fig. 1a), and can be expressed as The most striking consequence of this anisotropy, in conjunction with equation (11), is that there is a large asymmetry in the site-spectral weight for the site A(A'), B and C around the DP. Namely, the site B predominantly reflects the large-v F electrons in the steep slope of the cone (jEp); the site C mostly probes the small-v F electrons in the gentle slope (jE0); and the site A(A') probes the entire electronic states on average around the DP (see Fig. 1d-f and Supplementary  Fig. 1a,c). This asymmetry of the site-spectral weight results in a clear difference in the size of the j-site DOS, D j (E, z), for the different sites- 25)-where D j (E, z) (per valley and ET molecule) is defined as the q summation of n z j q ð Þ at a given energy in the band z, which is expressed as Here, V C is the 2D unit-cell volume in the conducting ab-plane. The contrasting features of the B and C site-spectral weights around the DP provide unique opportunities to probe the excitations of large-v F Dirac electrons (in the steep slope) and small-v F Dirac electrons (in the gentle slope) separately in terms of a site-selective local measurement such as NMR.
We note that in the original effective TB calculation by Katayama et al. 25 , the velocities are given in the unit of energy (meV), w 0 ¼ (w 0x , w 0y ) ¼ ( À 38.9, 4.8) and 43.9), because both the primitive vectors and the reciprocal lattice vectors are set to have a unit length in their notation. To recover the ordinary physical unit (length/time), one has to multiply the velocities either by a/' or b/' , using the values of the lattice constants a and b at the current pressure (2.3 GPa), where ' is the Planck constant divided by 2p. By linearly extrapolating the X-ray diffraction data obtained at 1.76 GPa (ref. 33) to 2.3 GPa, the lattice constants are estimated as a ¼ 8.567 Å and b ¼ 10.282 Å. The velocities in equation (12) are obtained in terms of these lattice constants.
Renormalization-group calculations. The observed nonlinear temperature dependence of the spin susceptibility below the inflection point T j flex in Fig. 2d,e, cannot be understood within the non-interacting Dirac-fermion picture, as we mentioned above (see Methods: effective TB model). In this temperature range, the screening effect should be weak reflecting the vanishing thermal excitations around the gapless point at E F . For this kind of situation, it is well known from the RG study of graphene Dirac cone that the LR part of the unscreened Coulomb interaction among electrons causes a logarithmic divergence of the Fermi velocity v F around the DP 7-14 . A similar argument has been recently proposed for the tilted Dirac cone in a-I 3 (ref. 46), in which a RG flow of the Fermi velocity is suggested as a function of T. Hence, the most straightforward and reasonable way to understand the low-T nonlinear feature of Fig. 2d would be to attribute it to the T-driven renormalization of v F due to the LR Coulomb interaction.
To check this hypothesis, we have performed a RG calculation based on the generalized Weyl Hamiltonian (equation (8)) and tried to fit the data. A circular momentum cutoff of L ¼ 0.667 Å À 1 around the DP (q ¼ 0) is introduced in the RG theory, which is of the size of the averaged inverse lattice constant, L ¼ 2p/L with L ¼ (a þ b)/2 ¼ 9.425 Å at 2.3 GPa. For the initial values of the velocities at the cutoff momentum q ¼ L (that is, v and w 0 in equation (8)), we employ velocities derived from the effective TB model of ref. 25, w TB 0 and v TB (equation (12)), as discussed in the previous subsection (Methods: Generalized Weyl Hamiltonian and site-spectral weight). In the one-loop order large-N expansion of the RG theory, v ¼ (v x ,v y ) are renormalized following equation (2) (given in the main text) and grows logarithmically as functions of L/q (where q is measured from the DP), whereas w 0 ¼ (w 0x ,w 0y ) are not renormalized (Supplementary Fig. 4a). We note that equation (2) takes into account the screening effect of the Coulomb interaction including the polarization bubbles in the self-energy. It is applicable to any size of the coupling g j ¼2pe 2 N=ð16 e ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi v 2 x sin 2 j þ v 2 y cos 2 j q Þ, where e is the dielectric constant and N ¼ 4 is the number of fermion species corresponding to the two DPs in the Brillouin zone and two spin projections (that means the twofold valley degeneracy is considered).
Reflecting the renormalization of v, the eigenenergy with the RG correction becomes (with the band index z ¼ ±) Correspondingly, the local DOS and the electron-spin susceptibility for the site j are respectively given by the expressions where V C ¼ 88.086 Å 2 is the 2D unit-cell volume in the conducting ab-plane estimated at 2.3 GPa, N ¼ 4 is the number of fermion species (again, the twofold valley degeneracy is included), E Z ¼ gm B H is the electron Zeeman energy and f(E) is the Fermi distribution. The integration with respect to q in equations (16) and (17) is done up to the momentum cutoff (q ¼ L). Note that T plays the role of the flow parameter in equation (17) (that is, a T-driven RG flow).
In Supplementary Figs 4-6, we present the calculated profiles of the velocities ( Supplementary Fig. 4), the j-site DOS ( Supplementary Fig. 5) and the j-site electron-spin susceptibility ( Supplementary Fig. 6) based on the RG equation (equation (2)) and equations (15)- (17). Here, to get a reasonable agreement with the experiment, we introduced a phenomenological parameter, u (r1), in the calculations which is defined by the expressions This parameter reflects a suppression of the velocity or a reduction of the hopping amplitude t ij between the lattice site i and j due to the SR part of the Coulomb interaction 53 , as mentioned in the main text. Then, the RG flow, which is determined by equation (2), is controlled by two parameters-the dielectric constant e and the phenomenological parameter u. Supplementary Fig. 4b,c, presents the parameter dependence of the flow of the velocities v x and v y . The dielectric constant e affects the power of the flow function v ¼ v(L/q) ( Supplementary Fig. 4b), while the parameter u determines both the initial values of the velocities and the power of the flow (Supplementary Fig. 4c). Supplementary Fig. 5b,c, depicts the e and u dependence of the DOS d D j E; z ð Þ for the site j ¼ A(A'). It is clearly seen that a prominent suppression of the DOS develops around E ¼ E F ( ¼ 0) by reducing the dielectric constant e (which corresponds to an increase in the coupling g j at the cutoff momentum q ¼ L). The reduction of the parameter u, on the other hand, leads to an enhancement of the DOS because the DOS is linked to the inverse square of the velocities 20 . In Supplementary Fig. 6b,c, we show the parameter dependence of the calculated electron-spin susceptibility. As the susceptibility probes the k B T-average of the DOS near E F ( ¼ 0), the suppression of the DOS for small e causes a reduction of the susceptibility at low temperatures ( Supplementary Fig. 6b). The enhancement of the susceptibility for a small value of u can be also understood in the same fashion ( Supplementary Fig. 6c).
Supplementary Fig. 7a shows the result of the variance analyses of the T-driven RG-flow fit to the experimental spin susceptibilities for the site j ¼ A(A') and C (Fig. 2d) based on equation (17). Here, the variance, plotted in the u À e À Var(w j s ) space, is defined by where T i is the experimental temperature points (i ¼ 1,2, ?, n), and w j s(Ti) and d w j s T i ð Þ are the measured and calculated susceptibilities at T ¼ T i , respectively. Note that the variance is defined only for the results on the site A(A') and C, since the results on the site B has less good agreement with the RG-calculation (see Supplementary Fig. 7b-e). This is because of the emergent ferrimagnetic spin polarization-the negative susceptibility on the site B at low temperatures (inset of Fig. 3c)-which is due to the SR part of the Coulomb interaction, not considered in this calculation. (For details, see Methods: Simulations with the Hubbard model.) The calculated susceptibilities at selected values of u ande are shown for three different sites j ¼ A(A'), B and C together with the experimental data in Supplementary Fig. 7b-e. Except for the site B, the agreement with the calculations and the experiments is pretty good for small values of u and e. The variance Var(w j s) has a minimum around (u, e) ¼ (0. 35,1) with little e dependence up to eE10 1.5 (Supplementary Fig. 7b,c) and increases rapidly as one moves away from this minimum especially when u is increased. Thus, we take these values of (u, e) ¼ (0.35, 1) as the optimal parameters hereafter. (Here, the result for e ¼ 1 is chosen because the e-dependence is small and affects the result little.) The small value of u ¼ 0. 35(o1) is in agreement with the aforementioned reduction of the hopping amplitude due to the SR repulsion 53  Lastly, it may be worthwhile to mention that the flow of the velocities v ¼ (v x , v y ) with w 0 ¼ (w 0x , w 0y ) staying constant in Supplementary Fig. 4a leads to a situation where v x , v y c|w 0x |, w 0y is realized for a large value of the momentum scale L/q (c1) (that is, in the vicinity of the DP). This means that the tilting term (w 0 ) becomes effectively negligible with respect to the anisotropy term (v) at low energy in the generalized Weyl Hamiltonian (equation (8)). Moreover, the values of the two velocities, v x and v y , remain very close to each other down to the vicinity of the DP (for instance, ln(L/q) ¼ 5 in Supplementary Fig. 4a corresponds to q ¼ 0.0045 Å À 1 in terms of L ¼ 0.667 Å À 1 ). These points together suggest that the anisotropy of the Dirac cone becomes very small near the DP and the cone is practically isotropic at low energy, as one can see in the calculated cone in Fig. 4d-f. This is can well account for the observed site dependence of the susceptibility, the root of which is linked to the tilt and the anisotropy of the cone. The site dependence becomes vanishingly small at low temperatures ( Supplementary Fig. 2), which reflect the renormalization of the velocity that makes the cone to be more and more isotropic at low energy.
Simulations with the Hubbard model. In the previous subsection, we have shown that the RG calculation based on the low-energy continuum model (equation (8)) well reproduces the observed nonlinear T-dependence of w j s T ð Þ at the site A(A') and C (Supplementary Fig. 7b-e) when the velocity suppression due to the SR Coulomb interaction is phenomenologically taken into account. However, the agreement is less good at the site B; in particular, the negative susceptibility w B s o0 below TE60 K cannot be reproduced at all, suggesting the presence of other interaction effects.
To understand the origin of the observed w B s o0, we have performed another simulation of the susceptibility based on the RPA. For this, we start from the standard Hubbard model with the onsite Hubbard interaction U, by extending earlier study 50 where a w jas is the creation operator on the site j ¼ A(1), A'(2), B(3) and C(4) in the unit cell a( ¼ 1, y, N u.c .) with the spin s( ¼ m, k), N u.c. is the total number of the unit cell, and t i a :jb is the nearest-neighbour and next nearest-neighbour hopping amplitude between the lattice site (i, a) and (j, b) (see Supplementary Fig. 8).
(It should be noticed that this Hamiltonian (equation (20)) lacks the site-dependent potential term (equations (3) and (4)) we employed above. From qualitative points of view, the inclusion of these potentials does not alter the results much, and, for the sake of simplicity, we shall omit this term in this subsection.) The amplitude t ia:jb at finite temperatures are estimated from the hopping integrals given by the first-principle calculations 24 at 8 K, {t LT p ; p ¼ a1 À a4'} (equations (6) and (7)), and at 300 K, {t HT p ; p ¼ a1 À a4'}, as given in the following (in eV) in combination with the interpolation formula given in ref. 50 Within the mean-field approximation, the diagonalization of the Hubbard Hamiltonian (equation (20)) yields where we define E ij k ð Þ ¼ P dij t ij e ikÁdij , d ij is a vector connecting the nearestneighbour lattice sites i and j, E Zs k ð Þis the eigenvalue (E 1s 4E 2s 4E 3s 4E 4s ), d iZs (k) is the corresponding eigenvector, f(E) is the Fermi distribution and m is the chemical potential. Note that the average electron number N js is determined self-consistently from the condition P js N js ¼6, reflecting the 3 4 -filling of the band. In the normal state, one has N j" ¼ N j# ; thus, the spin s is omitted hereafter. We introduce the bare site-spin susceptibility matrix,ŵ ð0Þ , whose (ij)-element is given by in terms of a form factor where id (d40) is an infinitesimally small imaginary part. Within the RPA approach, the spin fluctuations are estimated using the expression 50 (whereÎ is the 4 Â 4 unit matrix), and the total RPA j-site-spin susceptibility (for Q ¼ 0 and o ¼ 0) is given by Now, we decompose the bare susceptibility w ð0Þ ij into two parts (for the reason that will become clear below): w and the inter-band contribution to the total RPA susceptibility is expressed as w j;inter In Supplementary Fig. 3a, the calculated temperature dependence of the intra-band RPA susceptibility w j;intra RPA is shown for U ¼ 0.14 eV in comparison to the non-interacting case (U ¼ 0) for the site j ¼ A(A'), B and C. It is seen that the intra-band susceptibility w j;intra RPA for a finite value of U becomes always larger than the case for U ¼ 0 at all temperature and all site j. The inter-band contribution, on the other hand, is found to provide a site-dependent correction to the susceptibility ( Supplementary Fig. 3b). Namely, the inter-band RPA susceptibility w j;inter RPA gives a positive contribution on the site A(A') and C, whereas the contribution is negative on the site B. The negative inter-band contribution on the site B (w B;inter RPA o0) develops with increasing U and in turn leads to a negative total RPA susceptibility (w B RPA o0) above a critical U value of U C E0.12 eV. The position of the minimum shifts towards higher energies with increasing U (Supplementary Fig. 3c). By taking a value of U ¼ 0.14 eV, the minimum of the total RPA susceptibility w B RPA appears at around TE50 K, which agrees well with the experiment (inset of Fig. 3c). These results demonstrate that the SR part of the Coulomb interaction between electrons causes a ferrimagnetic spin polarization at low temperature. This leads to a situation where the site B is subjected to a negative local magnetic field, pointing opposite to the external field, while the other sites (A, A' and C) feel a positive field, as schematically illustrated in Fig. 6.
To have an overall comparison of the RPA calculations with the experiment, the first derivative of the susceptibility is analysed for the case U ¼ 0, the total RPA susceptibility w j RPA (equation (29)) for U ¼ 0.14 eV, and the observed spin susceptibility (Fig. 2e), as depicted for the non-equivalent sites in Supplementary  Fig. 3d-f. The calculations capture the experimental features relatively well on the site A(A') and C above the peak temperature (TE50 K), whereas the calculation does not agree with the experiment at all on the site B. At low temperatures, on the other hand, the agreement between the calculation and experiment becomes worse even for the site A(A') and C. That is, in the low temperature limit, the calculations show a saturation both for U ¼ 0 and the finite U, while the experiment exhibits a monotonic decrease towards zero (Supplementary Fig. 3d and f). We believe that this disagreement arises because the present RPA calculation does not incorporate the T-driven v F -renormalization effect due to the LR part of the Coulomb interaction, as discussed in the previous subsection (see Methods: RG calculations).
The v F -renormalization results in a super-linear temperature dependence in the spin susceptibility w j s ( Supplementary Fig. 6), which causes a decrease of the first derivative of w j s with decreasing temperature as reflected in the experiment. Finally, we comment on the comparison of the experiment with an orthodox RPA fitting. This is done by assuming a simplified RPA (s-RPA) expression for the spin susceptibility, which is defined by  (26) and (28), respectively, for Q ¼ 0 and o ¼ 0 with the Hubbard interaction U in equation (20) replaced by the site-dependent parameter U j .) Supplementary Fig. 14 presents the least-square fitting results to the experiment using the s-RPA expression, which yields U A(A 0 ) ¼ 0.6, U B ¼ 1.3 and U C ¼ 0.4 (in eV). It is clearly seen that the observed nonlinear temperature dependence of the susceptibility cannot be reproduced by the s-RPA fit at all sites. In particular, there is an unphysical divergence in the calculation, which is linked to the large U j values used in the calculation that are too large compared with the typical values applicable to this system 24,26,31,35,40,41,48,50 . It is thus concluded that one has to consider the full-matrix elements of equation (28) in order to obtain a reasonable agreement with the experiment.
Data availability. The data that support the findings of this study are available from M.H. upon requests.