Bulk Rotational Symmetry Breaking in Kondo Insulator SmB6

Kondo insulator samarium hexaboride (SmB6) has been intensely studied in recent years as a potential candidate of a strongly correlated topological insulator. One of the most exciting phenomena observed in SmB6 is the clear quantum oscillations appearing in magnetic torque at a low temperature despite the insulating behavior in resistance. These quantum oscillations show multiple frequencies and varied effective masses. The origin of quantum oscillation is, however, still under debate with evidence of both two-dimensional Fermi surfaces and three-dimensional Fermi surfaces. Here, we carry out angle-resolved torque magnetometry measurements in a magnetic field up to 45 T and a temperature range down to 40 mK. With the magnetic field rotated in the (010) plane, the quantum oscillation frequency of the strongest oscillation branch shows a four-fold rotational symmetry. However, in the angular dependence of the amplitude of the same branch, this four-fold symmetry is broken and, instead, a twofold symmetry shows up, which is consistent with the prediction of a two-dimensional Lifshitz-Kosevich model. No deviation of Lifshitz-Kosevich behavior is observed down to 40 mK. Our results suggest the existence of multiple light-mass surface states in SmB6, with their mobility significantly depending on the surface disorder level.

In Kondo insulators, the physics is controlled by the strong many-body interactions [1]. The hybridization between the localized f electrons and conduction d electrons causes the formation of Kondo singlets, which leads to a quench of local-magnetic-moment characteristics. Also, a narrow hybridization gap is developed at low temperature, resulting in a crossover from metallic to insulating behavior. In recent years, topological nontriviality is suggested to be hosted by Kondo insulators [2,3]. The opposite parity in the f band (odd) and d band (even) protects a band inversion similar to that in normal Z 2 topological insulators. In particular, the very large spinorbit coupling in the renormalized f electrons can give a system ground state with "nontrivial" topological order, i.e., a different topological invariant from that in vacuum. As a result, a gapless two-dimensional (2D) Dirac electron state, known as the topological surface state, has to exist at certain high-symmetry points in the surface Brillouin zone. Such predictions point the Kondo insulators out as promising candidates of interaction-driven topological insulators, subsequently make this family a focus of attention in condensed matter physics.
The cubic structured SmB 6 , the very first confirmed member of Kondo insulators [4], has been elaborately studied as the most feasible example of the electroncorrelated three-dimensional (3D) strong topological insulator [5][6][7]. A large amount of experimental observations on this material have been published [8], with some giving hints of the topological surface state [9][10][11], though the decisive evidence is yet to be found. The most Corresponding authors: L. Li (luli@umich.edu) striking discovery in SmB 6 is the complicated de Haasvan Alphen (dHvA) oscillations detected at low temperature where the resistivity shows insulating behavior followed by a plateau below 3.5 K. While our group has reported quantum oscillations corresponding to 2D Fermi surface (FS) and light carriers that are consistent with the expectation on a typical topological surface state in aluminum-flux-grown samples [12], another work based on the floating-zone(FZ)-grown sample claimed the oscillations have 3D characters thus bulk origin, and an abnormally enhanced quantum oscillation amplitude suggesting a deviation from the Lifshitz-Kosevich (LK) theory below 3 He temperature [13]. To make it more confusing, no quantum oscillations have ever been observed in transport measurements [12,14,15]. Several theories have been proposed to reconcile the puzzling experimental results [16][17][18] as well as to explain the enriched exotic low-temperature behaviors in SmB 6 [19][20][21][22][23][24]. The key problem to the confusion lies in the lack of a controlled study of the quantum oscillation amplitude.
In this work, we resolve the problem by mapping the oscillation amplitudes of DIFFERENT surfaces of the SAME crystal. We carried down magnetic torque measurements of flux-grown SmB 6 single crystals [25] down to 40 mK in a rotating magnetic field up to 45 T. Our result shows a broken rotating symmetry in the amplitude of the main dHvA oscillation branch, indicating a 2D nature of the electronic state. In addition, neither very high frequency oscillations with 3D behavior, nor abruptly enhanced dHvA amplitude suggesting failed a LK description is observed in any of our samples. These observations point to multiple 2D metallic states with small effective mass existing in Kondo insulator SmB 6 . FFT peaks are indexed with the same labels as in Ref. [12]. Solid symbols denote the data points at φ < 45 • , while hollow symbols are data taken at φ > 45 • . Harmonics of branches α, β and γ are presented by diamonds, circles and triangles, respectively. Dashed lines are fittings based on 2D FS model: F = F0/cos(φ − φ0). For branch β, φ0 = ±45 • and F0 = 285 T. The two "split" branches γ1 (blue) and γ2 (navy) both have symmetric axis along [100], i.e., φ0 = 0 • and 90 • , while F0 is 374 T for γ1 and 414 T for γ2.
By using the capacitive magnetic torque magnetometer shown in Fig. 1(a), we observe clear dHvA oscillations with the coexistence of different periods in our flux-grown SmB 6 samples (for details of sample preparation and experimental methods, see Appendix A). The angle-resolved field dependencies of the magnetic torques τ are shown in Fig. 1(b). We convert the measured capacitance C(B, θ) to torque by the relation τ ∝ 1/C. The absolute value of the torque signal is calibrated by a zerofield rotation in which the weight of sample m generates a change in the capacitance of ∆(1/C) ∝ mgl cos φ. Here l is the length of the cantilever beam and φ is the sample tilt angle between the magnetic field and the crystalline [001] direction as described in Fig. 1(a).
There are several interesting features in the angledependent behavior of τ (B). First, it is apparent that the non-oscillatory background of τ (B) changes its sign abruptly at φ = 0 • and 90 • , while at φ = 45 • there is another sign change but it is more smoother. We argue that this is most likely due to a bulk magnetic susceptibility anisotropy between the cubic [100] and [101] directions. Second, the dHvA oscillations on the torque curves also change sign at φ = 0 • , 45 • and 90 • . This "flipped" dHvA pattern across certain magnetic field directions (see τ (B) curves at φ = 40.5 • and 48.4 • in Fig.  1(b) as a typical reference) hardly reflects a sudden jump on the phase of quantum oscillation, while a more natural explanation is a direction change of the oscillatory torque vector τ = M × B that happens at these angles. Both a 3D electronic system with susceptibility anisotropy along [100] and [101], and a 2D diamagnetic system, can exhibit such an oscillatory torque flip at φ = ±(N/4)π (N can be any integer). Third, the amplitude of dHvA oscillation is considerably large. At above 35 T, the oscillatory torque ∆τ is roughly (0.5-1)×10 −6 N·m, corresponding to an effective magnetic moment ∆M ef f = ∆τ /µ 0 H of approximately (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) Such a large dHvA oscillation amplitude basically rules out the possibility of "false quantum oscillation" coming from the aluminum (Al) flux incorporated inside the sample (Appendix B), but also gives some difficulties to the 2D surface state interpretation. In a standard analysis on 2D electron systems, the magnetization oscillation has an amplitude upper limit of e n 2D /πm * for a unit area. This dHvA amplitude is usually much larger than that in real materials since there are several damping factors needed to be taken into account [26]. Here n 2D is the 2D density of carriers that contribute to the dHvA oscillations. Using the electronic parameters calculated in our earlier work [12], a 2D magnetic moment of ∼ 1×10 −9 A·m 2 is estimated for a surface area of 1 cm 2 . Giving the millimeter size of our sample and the effective magnetic moment ∼ 1×10 −8 A·m 2 , the discrepancy turns out to be roughly 2 order of magnitude. Moreover, the M ef f discussed above is only the component perpendicular to the magnetic field, M ef f = M ⊥ , therefore smaller than the total magnetic moment. The puzzling large amplitude of the dHvA oscillations has been, however, reported in the confirmed 3D topological insulator Bi 1−x Sb x [27], though the reason why the conventional estimation based on the surface carrier density failed there is unknown [28]. Further works looking into the peculiar magnetizing properties of topological surface state are needed to solve this question.
The fast Fourier transformation (FFT) of M ef f show results that are consistent with our previous work [12]. Three main branches F α , F β , F γ and their higher order harmonics are resolved. The FFT peak positions are plotted in Fig. 1(c) as a function of the angle between the applied field and the crystalline equivalent [100] directions in the cubic structure of SmB 6 . For F β and F γ , a fitting of F = F 0 / cos(φ − φ 0 ) can follow the behavior of F (φ) quite well, indicating a 2D nature of the related FSs. The value of φ 0 hints that pockets β and γ have the symmetric axes along the equivalent [101] and [100] directions, respectively. A "split" of peak γ is observed in a wide angle range which may indicates a subtle magnetic breakdown (for details, see Appendix C). These results are confirmed by repeated measurements in several samples and no sample dependence on the dHvA frequencies have ever been observed. Also, a recent tunneling spectroscopy study on SmB 6 single crystals reveals two Dirac-like surface bands on (100) surface and an additional one on (101) surface [29], which is in agreement with our observation of three bands. The smallest orbit α, however, has indeterminate dimension and geometry (Appendix C).
We do not resolve any dHvA frequencies higher than 2 kT (Appendix D). The high-frequency components observed in FZ-grown SmB 6 , which indicate large Fermi pockets with the size comparable to the area of Brillouin zone [13], are confirmed to be absent in flux-grown samples. It should be pointed out that there are still some similarities between the quantum oscillation spectra in our flux-grown samples and those in FZ-grown crystals (see Appendix E for details). In FZ-grown SmB 6 , low frequency oscillations were also resolved and assigned to small orbits ρ and ρ ′ [13]. ρ ′ shares the same angle range with our branch α and ρ is close to β and γ in our FFT spectra. As mentioned in Ref. [13], alternative possibilities are cylinder-like "neck" sections in a 3D electronic structure, or elongated ellipsoidal FSs. Similar elaborated comparison of the angle dependence of β/γ and ρ has been made in Ref. [17], which shows inconclusive result in distinguishing the effectiveness of 2D and 3D model. Therefore, the angular dependence of the oscillation frequencies cannot determine the origin of quantum oscillations in SmB 6 . As demonstrated later, the angular dependence of the oscillation amplitudes indicates that the oscillations most likely aries from the surface state. The angular dependence of the FFT amplitude of dHvA oscillation branch F β is plotted in Fig. 2(a). Here we use a tilt angle θ with different definition: the angle between H and the [101] direction of the crystal, as depicted in the inset of Fig. 2(a). As expected for an effective magnetization extracted from magnetic torque data, this amplitude will drop to zero if the total magnetization vector is parallel/antiparallel to the applied magnetic field which result in no "effective" component to be detected. In Fig. 2 Fig. 2(b) we summarize the angle-resolved dHvA amplitudes in another SmB 6 single crystal, S1, measured with the same experimental set-up. The twofold feature in Fig. 2(b) is much weaker, suggesting the symmetry breaking is more related to the sample instead of the cantilever magnetometry setup.
A reasonable interpretation of this symmetry breaking is the surface origin of F β . In our magnetic torque measurement, the signal from two parallel surfaces (e.g. (101) and (101)) will be picked up simultaneously, and in the field rotation in (001) plane, we can obtain the magnetic response from two individual sets of surfaces, which are perpendicular to each other. These two sets of surfaces are prone to have different plane impurity densities and subsequently differed carrier scattering rates that can apparently affect the amplitude of the quantum oscillation (Appendix F). Here, we analyze the angular dependence of dHvA amplitude of F β using a 2D LK model [26]: where M ⊥ is the effective magnetic moment picked up in torque measurement, µ(0) is the carrier mobility at θ = 0: µ(0) = eτ s (0)/m * (0), and ξ = πλB/µ(0). This model takes the 1/cos θ anisotropy of the cyclotron mass m * (θ) as the main contribution to the angular dependence of oscillation amplitude. The comprehensive simplification process of this model, in which we consider the anisotropy of each term in the LK formula in 2D case, is presented in Appendix F. We also applied Eq. 1 to the dHvA data of topological surface state as well as bulk state in Bi 1−x Sb x reported in Ref. [27], the fittings shown in Appendix Fig. 7(a) and (b) provide strong evidence that Eq. 1 is a valid model in describing the two dimensionality of electronic states and can effectively track the difference in the angle-dependent quantum oscillation amplitude between 2D and 3D system.
As shown by the fittings in Fig. 2 and Appendix Fig.  7(c), the two-fold symmetry in FFT amplitude can be well described by a difference of carrier mobility µ(0) on the two perpendicular sets of surfaces. By assuming a reduction on µ(0) of ≃ 11.5% ( Fig. 2(a)), the nearly 50% amplitude suppression is reproduced even for a 25% larger surface area (estimated from the sample geometry in the inset of Fig. 1(a)). This result substantially supports the 2D nature of oscillation branch β, as even an elongated ellipsoidal FS shows an evident deviation at high tilt angle (Appendix F). Furthermore, the effective fitting by using zero or small value of ξ reveals an ignorable Zeeman attenuation of the scattering rate (Appendix F).
The mobility difference is much larger between different samples, as in sample S2 µ(0) (Appendix Fig. 7(c)) is more than three times as large as in sample S5 ( Fig.  2(a)). This sample-dependent behavior is more likely due to the varied surface impurity level within samples. As we know, the scattering on the surface of SmB 6 is highly related to the surface disorder [30], and the dephasing length is differed by several hundred percent in different samples [10]. Also, the carrier mobility µ(0) we obtained from the fitting is only 50%-70% of that calculated from the Dingle plot (see Appendix G). This discrepancy can be addressed to the complicated electron scattering mechanism in SmB 6 which remains an enigma to be solved by further studies. We point out that the mobility attained by our dHvA amplitude fitting and Dingle analysis yield the same order of magnitude and are both much higher that those obtained from transport experiments [14,31].
Finally, the fitting based on the LK formula in Fig.  2 demonstrated again the validity of the Fermi Liquid theory in SmB 6 . The mysterious sudden enhancement of quantum oscillation amplitude in SmB 6 reported by Tan et al. [13] has attracted much attention as a rare 3D example of the deviation from LK formula, which has been suggested as a reflection of unconventional quantum oscillation [19-21, 32, 33]. In our dHvA studies down to 40 mK, however, such behavior is not repeated, even though we applied up to 45 T magnetic field, stronger than that used in Ref. [13]. In Fig. 3 we summarized the temperature dependence of dHvA oscillations at φ = 32.6 • : the capacitance C(B) (Fig. 3(a)), the oscillatory magnetic torque τ osc (inset of Fig. 3(a)) and the FFT curves of τ osc (Fig. 3(b)) all show almost no discernible difference at varied temperatures between 41 mK and 656 mK. None of the dHvA branches we resolved have any abnor- mal enhancement in this temperature range. Fig. 3(c) summarizes the evolution of the FFT amplitude. The light effective masses for both F β and F γ together with the low T give a temperature damping factor R T very close to 1. As a result, the oscillation amplitude is generally a constant with the relative change almost ignorable. We also tracked the evolution of the FFT amplitude of F β from 350 mK up to 30 K, at a tilt angle of φ = 41 • , as shown in Fig. 3((d). Assuming that the 2D LK formula the damping factor related to temperature is the same as that for the 3D case [34], the overall behavior can be fitted by LK formula with an effective mass of 0.138 m e , consistent with our former report [12]. The weak temperature dependence of dHvA amplitude below 300 mK is confirmed by measurements in different samples and at various tilt angles, which is shown in Appendix Fig.  9.
The distinct physical phenomena observed in fluxgrown and FZ-grown SmB 6 are quite intriguing and confusing. The FZ growth are reported to induce a small portion of Sm vacancies in SmB 6 single crystal [35,36,42]. Subsequently, a slight difference on the Sm valence on the surface can be established, which in turn modifies the Kondo interaction near the surface [16,17]. Evidences of incomplete Kondo coupling in SmB 6 , especially at the vicinity of the surface, have been discovered [14,[36][37][38]. Actually, in a topological Kondo insulator, both light and heavy surface states can be supposed to appear by varying the detailed band parameters [39]. Besides this valence variation scenario induced by non-stoichiometry in the framework of Kondo physics, another interpretation of the inconsistencies between our results and those reported in FZ-grown crystals by Tan et al. [13] is attributed to different disorder levels. Recent theory points out a topological Kondo insulator can be restored from an exotic Skyrme insulator, in which dHvA oscillation is contributed by scalar charge neutral particles, by introducing a certain degree of disorder [24]. The absence of high-frequency oscillations in our sample can also be attributed to the higher scattering rate induced by the impurity/defects. Furthermore the deviation of LK behavior in FZ-grown samples was taken as evidence for exotic nature of oscillations [13], we note that a similar double-step feature was observed in topological nodal semimetal ZrSiS [40], in which the Fermi surface nesting leads to two LK behavior with two different effective masses. This is a more realistic origin than the various unconventional quantum oscillation models [19-21, 23, 24, 32, 33]. Overall, the topological surface state explanation is still the most natural one for all of our observations of the magnetic quantum oscillations in flux-grown SmB 6 .
In summary, magnetic torque data of Kondo insulator SmB 6 measured in intense magnetic field at dilution refrigerator temperature has been investigated comprehensively. The amplitude of main dHvA oscillation branch β, which shows 1/cos φ angular dependence in frequency, displays a broken four-fold symmetry with field rotating in crystalline (010) plane. The angle-dependent oscillation amplitude can be fitted by a standard 2D LK model with respected to each sets of (101) planes. The carrier scattering rate obtained from the fittings is differed significantly between samples as well as surfaces on the same sample. The dHvA oscillations are also fully saturated at low temperature, implying small carrier masses. Our results indicate that multiple 2D light electron states are existed on the surfaces of SmB 6 . SmB 6 single crystals were grown by the aluminum (Al) flux method [25]. The chunks of Sm (99.95%), the powder of Boron (99.99%) and Al (99.99%) were mixed together with a mass ratio of 1:6:400, and then loaded into an alumina crucible. The entire mixture was heated to 1550 • C and then stayed at this temperature for 2 days before cooled down to 600 • C at 5 • C/h. During all the preparing and heating progress, the mixture was kept in the Argon gas. After cooled to room temperature, the samples with Al flux were soaked in the dense NaOH solution to remove the Al flux, and then washed by dilute HNO 3 solution. The samples were characterized by Xray diffraction(XRD) to determine the orientation. Upon cooling from room temperature to 3 He temperature, the resistivity of our SmB 6 samples is enhanced by 4-5 orders of magnitude, and a resistive plateau shows up below 3.5 K [15]. Data discussed in this work were mainly taken from SmB 6 sample S5, which has a size of 2.1×1.6×1.2 mm 3 and hosts large (100) and (101) surfaces. A smaller sample labeled as S6 was also measured but the results shown no significant differences and the signal quality is lower due to the smaller quantum oscillation amplitudes. All the samples we used are as-grown single crystals without cleaving or polishing.
The high magnetic field torque magnetometry measurements were carried out using the capacitance method in National High Magnetic Field Laboratory (NHMFL), Tallahassee.
The SmB 6 samples were glued to a beryllium-copper cantilever with the thickness of 0.025 mm. Variation in the capacitance between the cantilever and a fixed gold film reflects the bending of cantilever, from which the magnetic torque can be obtained. The set-up with sample S5 attached is shown in Fig. 1(a). Such devices were put into a rotator and then loaded into a dilution fridge in a hybrid magnet which can apply magnetic field up to 45 T. The cantilevers were rotated with magnetic field in the crystalline (010) plane.
We measured the capacitance change with sweeping magnetic field via two methods. The frequently used Andeen-Hagerling AH2700A digital capacitance bridge usually has a noise level of 10 −4 pF in our experimental environment, and the automatic balancing is slow, which means it may not be suitable to pick up the weak high-frequency dHvA oscillations in SmB 6 . As an alternative, we chose the General Radio analog capacitance bridge combined with the Stanford Research SR124 analog lock-in amplifier. By balancing the starting capacitance C 0 manually and reading the voltage change during field sweeping, we can achieve a better resolution with the noise level reduced by one order of magnitude. Also this allows for a continuous reading of the cantilever response. Extrinsic quantum oscillations introduced by Al flux trapped inside the sample has been reported in CaB 6 [41]. The existence of epitaxially oriented Al flakes in the flux-grown SmB 6 single crystals has also been confirmed, with a percentage of 2-4 wt% [42]. However, the dHvA amplitude of single-crystalline aluminum is known to be ∆τ ≃ 3.5 × 10 −7 N · m at 4.2 K under B = 2 T, for the strongest oscillation branch γ 5 in a sample with the mass of 45 mg [43]. Considering of the reported effective mass m * /m e = 0.18 and Dingle temperature T D =0.8 K for this band in aluminum, there is an amplifying factor of approximately 60x in torque for the condition of B = 40 T and T = 45 mK. It means the weight of incorporated Al, if it contributes to all the dHvA signals shown in Fig. 1(b), should be ∼ 2 mg. Since the weight of SmB 6 sample S5 is 12.9 mg, the amount of co-crystallized Al could be as large as 15 wt% (58.4 mol%) which is much larger than that revealed by X-ray diffraction study [42]. Consequently, the dHvA patterns in Fig. 1(b) is more likely to be intrinsic.
There are extra evidences against the flux-induced extrinsic quantum oscillations. First, there is no two-peak feature in any of the γ branches in the dHvA oscillations of aluminum [43,44], where as such a peak split has been observed in our flux-grown SmB 6 samples ( Fig. 4(a)). Second, there are evident discrepancies between the angular dependence of oscillation frequencies in SmB 6 and Al as reported before [12]. Also, the resemblance of dHvA frequencies in flux grown samples and floatingzone grown samples (which are grown free of Al, see Appendix E) supports that the oscillations are intrinsic.
Appendix C: The splitting on γ and the low frequency oscillation α There are two adjacent FFT peaks at the location of branch γ, as shown in Fig. 4(a), with a frequency interval of 43±5 T for most of the angles. We label the peak with lower frequency γ 1 and the high-frequency one γ 2 . Such phenomenon is also presented in our old data, though under lower magnetic field (up to 18 T) the splitting is less clear [12]. The inset of Fig. 4(a) gives a comparison of the FFTs with different end points, and the splitting appears to be unambiguous only for a cut-off field higher than 30 T. The origin of this observation is unclear. If the FS γ is a quasi-two-dimensional one and has a small periodic warping along k z direction, two extrema of the FS cross-sectional area can result in two quantum oscillation frequencies that are close to each other [45,46]. In this scenario, however, the two extreme areas have different curvature in the angle dependence and can cross together at certain angles. These expectations are absent in our data (see the triangle symbols in Fig. 1(c)). We suggest a more plausible explanation that the two Two peaks on the high-frequency side of the major peak β are labeled as γ1 and γ2, respectively. One unknown small feature on the φ = 48.4 • at 435 T is marked by a star. Inset: The FFT spectrum at φ = 28.6 • , between 11.4 T and a varied cutoff field. The split of peak γ is more evident as the cut-off field becomes larger. (b)The angle dependence of lowest frequency oscillation branch α and its second harmonic, plotted together with the frequency interval between split γ1 and γ2 and data of frequency ρ ′ extracted from Ref. [13]. The red curves are fittings using a 2D Fermi surface (FS) model, whereas the black curves are 3D fittings based on an ellipsoid FS. For the 3D model, we follow the previous report in which the geometry of the FSs contributing to ρ ′ are figured out to be small ellipsoidal pockets with their long axes along the cubic [101] directions [13] peaks γ 1 and γ 2 respectively come from two sets of electron orbits on (100) surfaces. Since the pockets α and γ have identical symmetric axes along [100], one possibility is that the split is originated in a magnetic breakdown between them, i.e., γ 2 = γ 1 +α.
We notice that the interval between the two split γ branches show insignificant change with varying field directions ( Fig. 1(c)). As exhibited in Fig. 4(b), the small angle dependence of γ 2 -γ 1 almost coincides with that of F α , which raises the possibility of magnetic breakdown between orbits γ 1 and α. Taking into account that the splitting of γ 2 is only apparent under high magnetic field (inset of Fig. 4(a)), magnetic breakdown is a rather promising explanation. With the orbit area differed by one order of magnitude, however, magnetic breakdown is unlikely to happen if γ and α are centered at the same momentum position in the Brillouin zone. It is probable that one of the two orbits is located away from the highsymmetry points in momentum space, but in this case additional frequencies like γ 1 +nα or α+nγ 1 are supposed to be observed according to the crystal symmetry. Actually, they are absent in our FFT spectra. At this stage, the reason of the splitting and the locations of orbits γ and α remain unclear.
The smallest orbit α resolved from the FFT of the torque curves has an indeterminate geometry. In Fig.  4 Fig. 4), which are identified as the 2nd harmonics here, only show up in the vicinity of [101] direction where the oscillation amplitude of α is the strongest. This behavior suggests those peaks in the range of 80-110 T are harmonics in nature instead of the diverging F α from one of the same set of FSs with the symmetry axis 90 • away (as β ′ and γ ′ plotted in Fig.  1(c) as well as in Ref. [12]). We also added the data re- FIG. 6. A comparison of the FFT peaks resolved from the magnetic torque data in the floating-zone grown SmB6 single crystals, extracted from Ref. [13], and the Al-flux-grown samples studied in Ref. [12]. In panel (a) and (b) the FFT peaks ρ and ρ ′ in Ref. [13] are plotted, respectively, together with the dHvA oscillation features in the flux-grown samples that within the same frequency range. In both data sets, the magnetic field is rotated from crystal [100] axis towards [011] axis. Dashed lines are fittings by the 2D cylinder FS model [12], and dot lines are fittings of the harmonics. ported in Ref. [13] to Fig. 4(b), though considering the different sample rotation direction in the two studies we only select the angles close to 0 • ([100]) and 45 • ([101]). According to Tan et. al., the smallest orbit ρ ′ is assigned to a small ellipsoid inside the "neck" connecting the large FSs [13]. However, judging from the consistency in Fig.  4(b) it is arguable that ρ ′ branch is corresponding to our α and its second harmonic.

Appendix D: Absence of high-frequency dHvA oscillations
Oscillation branches with frequencies higher than 2 kT, while have been detected in the floating-zone furnace grown SmB 6 single crystals [13], are completely missing in our measurement. Actually the frequencies between 1 kT and 2 kT are already rather weak in our FFT spectra. In Fig. 5 we show the analysis for two field orientations, i.e., φ = 40.5 • at which the low-frequency branch α is strong but β and γ are relatively weak, and φ = 84.9 • where branch β has a large spectral weight. At both angles the capacitance was measured by analog capacitance bridge. The oscillatory part of M ⊥ , the perpendicular component of magnetization, is dominated by the "slow" dHvA oscillations (Fig. 5(a) and (c)). After subtracting those "slow" components (with F < 1200 T), the residual magnetization term is barely noise with the amplitude approximately 10 −11 A·m 2 , that is, 0.1-0.2% of the total oscillatory M ⊥ (see the inset of Fig. 5(a) and (c)). No periodic small wiggles can be isolated that can indicate the existence of fast quantum oscillations. On the FFT of M ⊥ , everything with frequency higher than 2 kT sinks into the background noise that is roughly 1/1000 of the main peak height and no features can be resolved, even when we take the FFT in a high field range (Fig. 5(b) and (d)).
Appendix E: Comparison between the dHvA frequencies in flux-grown and floating-zone-grown samples In Fig. 6 we provide a comparison of the results reported by two groups on the dHvA frequencies resolved from SmB 6 single crystal samples grown via different approaches, i.e., Al-flux crystals studied by Li et. al. [12] and in floating-zone crystals studied by Tan et. al. [13]. Measurements are taken with the same technique and experimental conditions [12,13]. As is confirmed by Fig.  5, the high-frequency components with F > 2 kT are not detected in the flux-grown samples. On the other hand, the low-frequency FFT peaks in these two data sets are generally comparable, with data points falling into the same frequency range and the overall trend of angular dependence also show some similarities. It highly suggests that (i) those dHvA frequencies are intrinsic in SmB 6 (ii) the two works are looking into the quantum oscillations from the same electronic states.
Appendix F: Angle dependence of dHvA oscillation amplitude in 2D electron system For a 2D electron system, given the condition of F/B ≫ 1, the amplitude of longitudinal magnetization quantum oscillations can be approximately described by a 2D LK expression [34,[47][48][49]: where R T = X p / sinh(X p ), X p = 2π 2 pk B m * T /e B, R D = exp(−2π 2 k B m * T D /e B), T D = /2πk B τ Q the Dingle temperature, A a parameter proportional to quantum oscillation frequency F , and R S = cos(pπgm * /2m e ) the spin-splitting factor. Here m e is the free electron mass and τ Q the quantum oscillation relaxation time. In this work we study the magnetic torque, and the effective magnetization extracted from the torque signal is M ef f = M ⊥ . According to Ref. [26], which serves as the starting point of our dHvA oscillation amplitude analysis.
There are three independent parameters in Eq. F2 that are functions of the magnetic field tilt angle θ: F = F (θ), m * = m * (θ) and τ Q = τ Q (θ). The dHvA frequency F appears in the universal coefficient of all the harmonics, (1/F )(∂F/∂θ)A, in which A ∝ F . According to our fitting in Fig. 2(b), F β (θ) shows the typical characterization of a 2D FS, i.e., F β (θ) ∝ 1/cos θ. Therefore we have: The effective mass m * is included in all three amplitude factors of R T , R D and R S . For a parabolic-band system, the definition of effective mass is m * = (∂ 2 E/∂k 2 ) −1 , whereas, for a linear dispersive 2D electron system such as the surface state of topological insulator, we can use the expression as follows [28,50]: where S(E) is the cross sectional area of the 2D FS perpendicular to the field vector. Note that S(E) ∝ F varies as 1/cos θ but the energy dispersion and carrier density will not change with the sample rotation in magnetic field, assuming the band dispersion relation is invariant along all directions and there is no magnetic-field induced modification in any of the bands. It means for 2D electron systems we have m * ∝ 1/cos θ which is equally effective for the conventional parabolic and the Dirac-like band dispersion [27,51]. We checked the anisotropy of effective mass in two SmB 6 single crystals (the FFT plots of these two samples can be found in Fig S3-S5 in Ref. [12]), and the results for pocket β are m * /m e = 0.124 (0.122) at θ = 14.6 • and m * /m e = 0.140 (0.147) at θ = 34.8 • for sample S1(2). The relative offset of m * (θ) regards to the expected 1/cos θ behavior is therefore 1.9% and 4.2% for sample S1 and S2, respectively. Considering the sample misalignment and the LK fitting error, this result is quite reasonable, and we can also give an estimation of the effective mass at θ = 0 • : m 0 = m * (0) ≃ 0.120m e . To simplify the model, we take the Dingle factor R D as the only amplitude factor that is effectively influenced by the angle-dependent m * (θ). This approximation is actually sensible owing to the following reasons. Giving the small value of m 0 for pocket β and the low environment temperature, the angular dependence of R T is almost negligible. For m 0 /m e = 0.12 and T = 40 mK, the value of R T is very close to 1 with an offset smaller than 0.1% except for the angle range |θ − 90 • | < 1 • . In our measurement, the oscillation signal from β branch on one surface cannot be detected with the field θ > 75 • away from the normal direction of the relating surface, as shown in Fig. 2. Hence the factor R T can be safely treated as an constant in our fitting. We also left out the spin-splitting factor R S since we do not have a reliable estimation of the Landé factor g for the light Fermi pocket β. If this pocket is a topological surface state, this factor will hardly play any role in affecting the quantum oscillation amplitude, because R S comes from the superposition of oscillations from the split Landau levels (i.e., spin-up and spin-down) [26]. In a topological surface state there is no spin degeneracy at k = 0, accordingly the Zeeman effect shift the position of the Landau levels instead of causing the splitting [27,28,52], consequently results in no reduction on the amplitude of oscillation.
The angular dependence of relaxation time τ Q (θ) is more complicated. We separate the magnetic field into two components, the in-plane field H = H sin θ and the out-of-plane field H ⊥ = H cos θ. The first term is known to have no significant transport response from topological surface state in the absence of the hybridization between top of bottom surfaces [53], which can be completely neglected in our bulk single crystals. Theoretically, the in-plane magnetic field will only shift the position of surface Dirac point in momentum space [54,55] and cause a net in-plane spin polarization [56]. A deformation of the Fermi pocket corresponding to the spin density redistribution is also suggested [57]. In all, the spin momentumlocking and the prohibition of backscattering in topological surface states is not destroyed in an in-plane magnetic field. The case is totally different for the second term, the out-of-plane component, which can break the timereversal symmetry and lift the protection of the topological nontriviality. In this case, the back-scattering is reintroduced, and the electron-impurity scattering is enhanced by the Zeeman-energy-related spin-canting [58]. The transport scattering rate takes the form as: λ is a system parameter related to g-factor and Fermi energy.
Taking into account all the angle-dependent parameters discussed above, we can give a fitting model of the quantum oscillation amplitude of M ⊥ for the fundamental harmonic: Ref. [27]. Tilt angle θ is the angle between the magnetic field, which is rotated in the binary plane, and the crystal axis C3. The carrier mobility in each panel is calculated from the parameters obtained in the same work. (c) Fitting of the angle-dependent amplitude of FFT peak β in SmB6 sample S2 by Eq. F6, with parameter ξ = 0. Definitions of µ0 and ξ are the same as in Fig. 2. Data are extracted from Fig. S4 in Ref. [12].
where µ(0) is the carrier mobility at θ = 0: µ(0) = eτ s (0)/m * (0), and ξ = πλB/µ(0). This expression is the same as Eq. 1. One needs to mention that the transport relaxation time τ tr in Eq. F5 is different from the quantum oscillation relaxation time τ Q we used in the Dingle factor, as the former one is more sensitive to backscattering [59]. Also, the field dependence of τ tr in Eq. F5 is basically a weak-field approximation. As a simplification, in the fitting model Eq. F6 we assume that τ tr and τ Q share the same field dependence in the magnetic field range in our measurement. The detailed field effect on the scattering in topological surface states still needs further investigations.
To examine the validity of our model, we applied it to the quantum oscillation amplitudes in a well known topological insulator Bi 1−x Sb x , reported by Taskin and Ando, in which both 2D and 3D FSs can be resolved in dHvA measurement [27]. Since the data in Ref. [27] was taken by SQUID magnetometer, the oscillations were detected on longitudinal magnetization, consequently the fitting model Eq. F6 should be modified to: Here we drop the field-dependent mobility term due to lack of information. Relying on the electronic parameters presented in Ref. [27], the carrier mobilities of the 2D surface state (F 1 in Fig. 7(a)) and 3D bulk state (F 2 in Fig. 7(b)) are 5.5 and 1.2 m 2 V −1 s −1 , respectively. Taking these mobilities as fitting parameters in Eq. F7, the curve in Fig. 7(a) can roughly track the fast decrease of dHvA amplitude when field is rotated towards the bisectrix plane down to θ ≃ 20 • . However, in Fig.  7(b), the attenuation of dHvA amplitude is obviously much slower than that expected in our 2D model at θ < 50 • . It should be mentioned that the frequency F 2 comes from a highly anisotropic ellipsoidal FS in Bi 1−x Sb x : the length of its longest semiaxis (along crystal axis C 3 ) is 8.5 and 17.7 times of the other two semiaxes, respectively [27]. The fittings in Fig. 7 indicate that even for such an extremely elongated FS, our model can effectively distinguish it from a real 2D cylinder FS. The 2D LK formula we used in data fitting is theoretically a low magnetic field approximation. In a 2D system, if the number of electrons is kept constant, the chemical potential is prone to be pinned in the highest occupied Landau level and also oscillates with increasing field [26]. This chemical potential oscillation will cause strong deviation from LK theory at high field, in which the position of chemical potential is assumed to be fixed [47,60]. Such deviation has been found in quasi-2D organic compounds [34,47,61] as well as in cuprate high-temperature superconductor YBa 2 Cu 3 O 6+x [62]. In topological insulators, however, this behavior has not ever been reported or discussed. An analytical analogue of LK formula, with well defined R T and R D , has been established for 2D Dirac-like electron system [63], and is effectively used in graphene [50]. Conventional LK analysis is widely applied and accepted in the study of quantum oscillations in topological insulators [27,28,64,65]. As for our data, the upper limit of magnetic field B (45 T) is much lower than the oscillation frequencies F β and F γ , therefore the low-field condition ∆E n ≪ E F is fulfilled for these two bands (here ∆E n is the energy interval between Landau levels and E F the Fermi energy). Besides, the sharp saw-tooth-like oscillation patterns expect for clean 2D systems are missing in our measurements, and the LK description works well in fitting the temperature dependence of oscillation amplitude [12]. Given the reasonable modeling of the angular dependence of the oscillation amplitude in the surface state of Bi 1−x Sb x , we conclude that the LK model in our analysis is a correct model.
For the value of B in the fittings, we still use the averaged inverse field [66] as the effective value: B ave = 17.48 T in Fig. 2(a) and B ave = 6.54 T in Fig. 2(b) and Fig. 7(c). For all three samples, the fittings are reasonably good, though not perfect, with the parameter ξ = 0. We also made curves with finite values of ξ and other fitting parameters unchanged in Fig. 2(a). It appears that an acceptable value is ξ 0.1, corresponding to λ 5.89 × 10 −5 and the enhancement on 1/τ is less than ≃ 12% at 45 T. The small Zeeman effect in scattering rate is qualitatively consistent with the calculation for Bi-based 3D topological insulators [58]. We note that with an appreciable Zeeman effect there will be visible splitting of the peaks/valleys in the dHvA oscillation patterns. Early work also suggested that a large Zeeman effect makes the Landau Level indexing plot non-linear [67]. None of these effects are observed in our results of the quantum oscillation patterns in SmB 6 [12]. Nonetheless, the small mismatch in the fittings in Fig. 2 can be assigned to the subtle effects that are ignored in our model such as the Zeeman term, as well as the sample misalignment in the measurement.
The effective fittings by a 2D LK model (Eq. F6) is an essential evidence against the bulk origin of F β . As mentioned above, our model describes a fast amplitude damping with field rotating away from the symmetric axis of FS. The elongated 3D FS in Bi 1−x Sb x with FS cross-sectional areas 8.5 times different between two perpendicular directions shows apparent deviation from the fitting curves ( Fig. 7(b)). The supposed 3D orbit ρ in FZ SmB 6 samples [13], which shares the same frequency range with our β branch, has a cross-sectional area difference of ∼ 3.3 between [101] and [101] directions. Such a moderate anisotropy cannot give the fitting results shown in Fig. 2 and Fig. 7(c).
Appendix G: Dingle plot and scattering rate difference between samples Figure 8 gives a brief summary of the Dingle temperature T D in the three samples S5, S1 and S2, derived from the slope of ln(∆M ⊥ /R T ) versus 1/B. The linearity of Dingle plot suggests the field modification on τ Q is almost ignorable. While the quantum oscillation mobility µ = e /2πk B m * T D is 50-100% larger than the fitting parameter µ(0) in Fig. 2, the relative magnitude of mobility within the three samples is the same for the two approaches. The discrepancies between mobilities attained by different experimental methods is a famous conundrum in SmB 6 . In transport measurement, much lower mobilities have been reported, which vary from several tens of cm 2 V −1 s −1 [68] to 120-140 cm 2 V −1 s −1 [14,31], and no quantum oscillation has ever been observed. These confusing phenomena may suggest complicated scattering mechanism in this material. Nonetheless, our magnetic quantum oscillation experiment has clearly proved that light carriers with relative high mobilities reside in SmB 6 , most likely on the surfaces.
The difference of dHvA amplitudes among samples is appreciably large. While sample S5 shows ∆M ef f with a magnitude of 10 −8 A·m 2 ( Fig. 1(b)), signals from other samples can be one to two orders of magnitude smaller with the surface area within the same order of magnitude [12]. It is also apparent that the signal strength is not proportional to the related surface area on one sample ( Fig. 2(a)). Apart from the most likely reason of surface impurity effect on the carrier mobility, we are also aware of the complex surface reconstruction in SmB 6 [69][70][71]. The multiple surface phases are possible to give different contribution to quantum oscillation. This information is not included in the fittings in Fig. 2  most temperature-independent behavior between the base temperature of the dilution fridge (40-45 mK) and 300 mK. No considerable low-temperature dHvA amplitude increase [13] has been observed. Figure 9 shows the result taken at φ = 88.2 • , i.e., field close to [100] direction, in Sample S6. Similar to the observation shown in Fig. 3, the oscillatory magnetic torque curves at all temperatures overlap with each other. The dominating frequency at this tilt angle is F β = 397 T as shown in Fig. 9(b). The amplitude attenuation of this peak is within 0.2% from 45 mK to 296 mK, in contrast with the previous observation of > 80% reported by Tan et. al. [13] in floating-zone-grown crystals. At this stage we confirm that such steep increase in dHvA amplitude does not exist in our flux-grown samples.