Resolving the topological classification of bismuth with topological defects

See allHide authors and affiliations

Science Advances  01 Nov 2019:
Vol. 5, no. 11, eaax6996
DOI: 10.1126/sciadv.aax6996


The growing diversity of topological classes leads to ambiguity between classes that share similar boundary phenomenology. This is the status of bulk bismuth. Recent studies have classified it as either a strong or a higher-order topological insulator, both of which host helical modes on their boundaries. We resolve the topological classification of bismuth by spectroscopically mapping the response of its boundary modes to a screw-dislocation. We find that the one-dimensional mode, on step-edges, extends over a wide energy range and does not open a gap near the screw-dislocations. This signifies that this mode binds to the screw-dislocation, as expected for a material with nonzero weak indices. We argue that the small energy gap, at the time reversal invariant momentum L, positions bismuth within the critical region of a topological phase transition between a higher-order topological insulator and a strong topological insulator with nonzero weak indices.


Topological materials host unique boundary states that cannot be realized as stand-alone electronic systems but only on the boundaries of a topologically classified bulk. These include the chiral edge modes of the quantum Hall state, the helical modes on the surface of a topological insulator (TI), Fermi-arc bands on the surface of a Weyl semimetal, and Majorana modes in a topological superconductor (18). Bulk topological defects in the crystalline structure also exhibit a similar bulk boundary correspondence. The surface termination of a bulk topological defect necessarily hosts unique crystallographic structures. One example is the surface step edge (gray shaded in Fig. 1A) that emanates from the surface termination point of a bulk screw dislocation (SD). The only way to realize a step edge that terminates at a point on the surface, away from the surface boundaries, is by its continuation into the bulk in the form of an SD. The length of the Burgers vector, which enumerates the winding of the bulk SD, dictates also the height of its corresponding step edge (9). Crystalline topological defects, such as SDs and disclinations, also serve as effective boundaries and can potentially bind topologically protected electronic boundary modes. Soon after the discovery of the symmetry-protected topological electronic phases, a general theory of topological defects and their associated gapless modes has been developed (1016). Examples of these bound states are the Majorana modes at vortex ends in chiral superconductors (17, 18) and one-dimensional (1D) helical modes at SDs in weak TIs, as proposed by Ran et al. (10). The latter case of an SD has been reinvestigated later on by Slager et al. (14). They proposed a classification scheme that takes into account not only time-reversal symmetry but also crystalline space-group symmetries (19). Hence, it takes into account additional cases where gapless modes are formed but are not included in the formulation proposed by Ran et al. (10). In this classification, dislocations play a fundamental role as probes to distinguish different classes of topological states (14, 16). Although theoretically it is quite established that SDs should bind helical modes in certain types of TIs, direct experimental signatures for their existence are still lacking. Recently, an experimental signature of increased conductivity has been reported in Bi1−xSbx single crystals that was attributed to the introduction of SDs by plastic deformation (20). Here, we use scanning tunneling microscopy (STM) to locally detect SDs on the (111) surface of bismuth (Bi). Using careful spectroscopic mappings, we explore the local response of the various boundary electronic states in Bi(111) to the unique surface perturbations imposed by SDs. We provide first direct spectroscopic evidence for the existence of helical modes on an SD. The correspondence between the structure of the topological defect and the topology class of its associated electronic mode (10, 12, 14, 15, 19, 20) further allows us to determine the topological classification of Bi as a strong TI with nonvanishing weak indices.

Fig. 1 Surface and edge spectrum of Bi(111).

(A) Schematic illustration of an SD in a lattice. (B) Bulk Brillouin zone (BZ) of rhombohedral Bi and its (111) surface–projected 2D BZ and zigzag edge–projected 1D BZ, with all TRIM points indicated. (C) Bulk band structure of Bi obtained from ab initio calculations using the modified Becke-Johnson (MBJ) exchange potential and generalized gradient approximation shown by the solid and dashed lines, respectively. Inset shows the band structure around the energy gap (5 meV) at L calculated with MBJ. (D) Large-scale STM topography of a Bi terrace with crystallographic step edges. Inset shows atomically resolved image overlaid with Bi bilayer crystal structure. (E) dI/dV spectroscopy given in arbitrary units (au) as measured at a point on the surface far away from impurities. (F) dI/dV map measured normal to the step edge (marked with a black arrow) along the gray dotted line in (D) showing the modulated LDOS and the 1D topological edge mode (marked by curly brackets). (G) Fourier transform of (F) showing the energy dispersion of scattering wave vectors along Γ¯M¯ and the corresponding calculated spin-selective scattering probability (SSP). (H) Calculated zigzag edge–projected LDOS showing 1D edge states dispersing along ΓM (see fig. S4).

Bi, an extensively studied low-density semimetal, has a rich history of intriguing quantum phenomena like electron fractionalization (21), nematic quantum Hall phases (22, 23), and unconventional superconductivity (24). Because of its large spin-orbit coupling, it plays a fundamental role in several topological materials, such as the first realized 3D strong TI, Bi1−xSbx (3, 25, 26), which shows almost identical 3D band structure to that of pure Bi (27), and the canonical TIs, Bi2Se3 and Bi2Te3, which host a single surface Dirac band (2831). Nevertheless, the topological classification of pure Bi has remained, thus far, rather ambiguous. Many theoretical models (3, 3234) indicate that Bi has a trivial Z2 topological classification. Since bulk Bi is inversion symmetric, its strong topological index, ν0, and the weak indices, ν = (ν1, ν2, ν3), can be calculated from the parity eigenvalues of the occupied bands at time reversal invariant momentum (TRIM) points (Γ, T, and three inequivalent L and X points, shown in Fig. 1B) (3, 32). While most theoretical calculations suggest that both the strong and the weak indices are zero (3, 32, 34), some calculations argue for nonvanishing strong and weak indices (35, 36). The origin of the ambiguity in the topological classification of Bi lies in the extremely small energy gap at the L TRIM point (34), which in calculation is found to be on the order of 5 to 15 meV (Fig. 1C). Experimentally, several studies demonstrate phenomenology, which implies nontrivial topology in Bi, both in bulk crystals and in thin films, as well as in its 2D limit of a bilayer (35, 3742). Direct visualization of the electronic structure of Bi around the M¯ surface energy gap (39, 42) suggests that the 2D surface bands may have topological origin. A bilayer of Bi was among the first materials to be classified as a 2D TI (43). Later, 1D edge modes were observed in thin films of Bi (40, 4446) and on step edges of cleaved (111) surfaces of bulk Bi crystals (38). Various interpretations have been ascribed to the step-edge 1D edge modes. These include trivial spin-split bands of the 1D Rashba type (47), helical edge modes of a 2D TI (38), and, more recently, helical hinge modes of a 3D higher-order TI (41). Helical 1D edge modes are counter-propagating time-reversed partners of each other, which form a Kramer’s pair. As a consequence, it is impossible to introduce any backscattering between the two states, unless time reversal symmetry is broken. This is the origin of the topological protection of helical edge states. We show here, by the aid of SDs, the robust behavior of the 1D edge mode in Bi both at the SD and along the step that emanates from it. This phenomenology of the 1D edge mode at the SD is consistent with a topological classification of Bi as a strong TI with nonzero weak indices. The small L gap positions Bi on the verge of a topological phase transition accompanied by a gap closure between a strong TI with nonzero weak indices and a higher-order TI. We argue that approximately helical 1D edge modes, on the step edge and the SD, appear at the critical vicinity of the topological phase transition on either side of it (48, 49).


We measured Bi samples derived from crystals of GdPtBi grown in Bi flux (see Materials and Methods) (50). We have found that cleaving these crystals under ultrahigh vacuum conditions repeatedly exposes surfaces of residual crystalline Bi due to its anisotropic rhombohedral structure. The slight lattice mismatch with GdPtBi lends the opportunity to investigate the properties of bulk Bi in the vicinity of SDs. Nevertheless, apart from the nucleation of isolated SDs, the Bi atomic surface structure and band structure we measured were found to be rather unaffected by the GdPtBi substrate (see section S1). The samples were measured at 4.2 K in a commercial Unisoku STM. A representative topographic image of a freshly cleaved surface, consisting of a flat terrace terminated by step edges along high-symmetry crystallographic directions, is shown in Fig. 1D. The atomically resolved topographic image shows the underlying triangular lattice indicative of a pristine (111) cleave plane of Bi, thus allowing its spectroscopic studies.

Surface electronic structure of Bi(111)

We next identify the spectral features of the (111) surface of Bi. A typical differential conductance (dI/dV) spectrum measured on the terrace, away from any scatterers, is shown in Fig. 1E. The various peaks in the dI/dV spectrum represent the extrema of the 2D surface bands expected on Bi(111) (38). The intermediate dips around 100 and −100 meV correspond to a suppressed local density of states (LDOS) within the bulk gaps at the T and L TRIM points, respectively (calculation in Fig. 1C). To resolve some of these surface bands, we have measured the dI/dV map, shown in Fig. 1F, taken normal to the step edge (along the gray dashed line in Fig. 1D). It shows strong spatial modulations that disperse in energy. These modulations are quasiparticle interference (QPI) patterns the electrons embed in the LDOS due to scattering off the step edge. Its Fourier analysis (Fig. 1G, left) reveals clear dispersing modes that agree well with the calculated spin-selective scattering probability (shown on the right panel). This detailed comparison indicates that the dispersing QPI patterns originate from the 2D surface bands that disperse between the Γ¯ and M¯ surface TRIM points in full agreement with previous studies of the (111) surface of bulk Bi crystals (38, 51). The topological nature of these 2D surface bands depends on the topological classification of Bi: Trivial surface states in the case of Bi are higher-order TI, while topological Dirac states in the case of a strong TI. The qualitative difference between the two is restricted to the immediate vicinity of the small L gap and merely amounts to whether they all disperse there into the bulk valence band or, across the gap, between the valence and conduction bands. This phenomenological distinction in the 2D surface states is hard to resolve spectroscopically due to the small L gap. However, the response of the 1D edge mode, found on Bi step edges, to SDs allows lifting of this ambiguity and revealing the topological classification of Bi.

We therefore investigate the 1D edge mode on step edges and its topological origin. Amid all the energy-dispersing QPI modes in Fig. 1F, a nondispersing increased LDOS is bound along the step edge. This 1D bound state spans the energy interval of both the T bulk gap at 30 to 180 meV and the L bulk gap at lower energies around −100 meV (both are marked with curly brackets). This extended energy dispersion is consistent with our ab initio calculations showing (red line in Fig. 1H) that the edge mode indeed disperses from Γ to M (the edge projections of the bulk T and L TRIM points, respectively; see Fig. 1B) and extends throughout the −200- to 200-meV energy range. The topological nature of such a 1D edge mode in Bi has been discussed in several previous reports (37, 38, 41, 46, 52), where various distinct origins have been ascribed to its formation including 2D TI (38), higher-order TI (41), and nontopological Rashba metal (47). However, all of these accounted for the appearance of the helical edge mode within the T gap alone and disregarded its extended energy dispersion that we note here down to lower energies where the L energy gap resides. This is true also for the energy dispersion of the 2D surface modes in Bi(111). These energy and momentum extents suggest that the L energy gap has an important role in the topological classification of Bi.

Spectroscopic mapping of an SD

To resolve the topological origin of the extended edge mode in Bi, we examined its response to crystallographic defects. A topographic image of a single Bi bilayer high step edge, terminating at an SD on the Bi(111) surface, is shown in Fig. 2A. The height of the step edge associated with an SD can generally be of several unit cells, nu. This height reflects the SD Burger’s vector, b = nu, which classifies the topological invariant of the defect. We map the dI/dV spectrum along the last atomic row of the step edge, across the SD and onward, extending into the terrace (along the dotted line in Fig. 2A). An atomically resolved topographic profile along this line cut is shown in Fig. 2B, further revealing a missing atom along the step edge. In the corresponding dI/dV spectrum (Fig. 2C), we find an increased density of state all along the step edge at the same energies that show a suppressed density of states within the T and L gaps on the open terrace (marked with curly brackets). The added density of states, contributed by the edge mode, is nicely captured in the dI/dV profiles in Fig. 2 (D and E) from the two energy ranges (marked by dotted lines in Fig. 2C), signifying the presence of the 1D edge mode. The density of states of the edge mode remains fairly constant down to the last atom before the SD, as well as at the vicinity of the vacancy on the step edge, although both serve as strong scatterers for the 1D electronic states. The small amplitude modulations in dI/dV along the step edge in Fig. 2 (D and E) are commensurate with the atomic positions resolved in topography in Fig. 2B and accordingly do not disperse in energy. The robustness of the 1D edge mode stands in sharp contrast to the scattering of the 2D surface states we find on the open terrace both from the remote step edge and from the SD, giving rise to a faint quantized pattern in Fig. 2C (see also fig. S2).

Fig. 2 Topographic and spectroscopic mapping of an SD.

(A) Topographic image of a clean (111) surface consisting of an SD (solid black arrow) with a single Bi bilayer long Burger’s vector. (B) Topographic line cut along the black dotted line shown in (A). (C) Spectroscopic map measured along the dotted line in (A) showing QPI pattern for the 2D surface states and a nearly unhindered edge mode on the atomically corrugated step edge terminating at the SD. (D and E) dI/dV line cuts along the black dashed lines shown in (C). The shadowed regions in (B) to (D) correspond to the segment where the STM tip images the atomic vacancy states and the segment at which it descends from the top to the bottom terraces where the tunneling matrix element from tip to sample may change.

Since the spatial distribution of the wave function of the edge state can vary with energy, we also spectroscopically map the surrounding of the SD. We overlay in Fig. 3A the 2D dI/dV mappings of the energy-resolved LDOS at the SD vicinity on a topographic height map (see also fig. S3). At high energies (Fig. 3B), we find metallic surfaces on the SD, step edge, and adjacent terraces. Below 180 meV (Fig. 3C), we detect the edge mode localized on the edge of the top terrace, all across the step edge, down to the SD. In a broad range of energies, we observe an increased LDOS right on the SD (Fig. 3D), which agrees with the continuation of the mode into the SD. We observe similar behavior also at lower energies around −100 meV corresponding to the L energy gap (see fig. S3), where the mode seems to be localized slightly lower on the step edge.

Fig. 3 Imaging the localization of the edge mode in the vicinity of the SD.

(A) Topographic image of the region around an SD. (B to D) Spatially resolved differential conductance map in the vicinity of the SD. The increased LDOS on the step edge and the SD shown in (C) and (D) correspond to a segment of the 1D edge mode. The black dashed line represents the position of the line cut shown in Fig. 2C.

The absence of backscattering and hybridization of the 1D edge mode along the step edge followed by its apparent channeling into the SD provides two important observations: It presents first spectroscopic signature for the long sought prediction of perfectly conducting channels on an SD, and it reflects the phenomenology of a strong TI with weak indices rather than of a higher-order TI. For a material characterized by nonvanishing weak indices, ν = (ν1, ν2, ν3), a bulk SD with Burger’s vector b is theoretically expected to bind a helical mode if ν · b = π(mod 2π) (10, 14). When the space-group symmetry is fully taken into account, there are cases in which the gapless SD states are formed even when this condition is not satisfied (14). Yet, for Bi with a space group of R3¯M, the above condition can be simply applied for the formation of helical modes (20). Hence, the expected phenomenology on a single unit cell high step edge, terminated by an SD, is of a 1D helical mode that runs along the step edge and continues into the SD line, with no backscattering and without gapping, as guaranteed by time reversal symmetry. This is consistent with our observations. In contrast, a higher-order TI is characterized by a null strong index, ν0 = 0, and null weak indices, ν = (0,0,0). Hence, an SD in a higher-order TI is not expected to bind a helical mode, unless, as was shown recently, its Burger’s vector is of a fraction of a unit cell height (53). Therefore, the expected phenomenology in the case of a higher-order TI is of even number of helical modes on the step edge that would presumably gap out or shows some signatures of back scattering near the impurity and the SD. Even if the pair of hinge modes avoids hybridization along the step edge, the SD would necessarily scatter them one to another and lead to their gradual gapping in its vicinity. The absence of these signatures, and moreover the existence of a slight increase in the LDOS on the SD, favors their identification as helical modes of a strong TI with nonzero weak indices rather than as higher-order TI hinge states.

While the phenomenology in the vicinity of the SD strongly suggests the existence of nonvanishing weak indices, this is inconsistent with most of theoretical models of Bi (3, 32, 34, 41), including our ab initio calculations in the absence of strain. These find the same band ordering at all TRIM points in Bi, leading to trivial Z2 topological indices, ν0 = 0, ν = (0,0,0). These conclusions, however, rely on the L energy gap that in calculation is found to be as small as 5 meV (Fig. 1C). This small energy gap makes the calculation of the band ordering at the L point challenging. Recent studies suggest, for instance, that this energy gap is in fact inverted with respect to the energy gap at the T TRIM point, which classifies Bi as a strong TI, ν0 = 1, with nonzero weak indices, ν = (1,1,1) (35, 36, 42). In this case, the weak TI-associated helical edge modes are expected to bind to SDs, and the 2D surface states of Bi are identified with the Dirac surface band of a strong TI, providing a compelling unified resolution to the combined boundary phenomenology observed in Bi.

Strain analysis

We next characterize in detail the SD and its effect on its surroundings. We extract the strain field induced by the SD directly from the gradient of the measured topography, ∇z(x, y). Its phase, arg ( ∇z), shown in Fig. 4A, confirms the azimuthally oriented strain profile that winds around the SD axis. Its magnitude, ∣∇z∣, displayed in Fig. 4B, shows a radially decaying strain profile as the displacement of the lattice, b, distributes over growing circumference. Cuts of the strain magnitude along the two directions, marked with arrows in Fig. 4C, demonstrate the anisotropy of the strain distribution. Ninety percent of the strain relaxes along the direction of the step edge within 8 nm (blue arrow), while 60° away from this direction (orange arrow), it relaxes within only ~3 nm. This anisotropic strain pattern on the surface suggests that the SD axis is tilted away from the <111> direction normal to the surface. A single bilayer long Burger’s vector along the <111> direction is incompatible with the ABC stacking of the Bi bilayers (Fig. 4C, inset). A surface normal SD would thus result in excess strain and would be energetically unfavorable. The two observations are consistently resolved by asserting an SD oriented with the <001> direction, along which the periodicity is indeed of a single bilayer. Such an SD is also expected to bind a helical mode as ν · b = π.

Fig. 4 Strain analysis in the vicinity of the SD and modeling.

(A) Phase and (B) magnitude of the shear strain field around an SD calculated from the gradient of the topography in Fig. 3A, showing its azimuthal radially decaying profile. The amount of strain within the dashed line in (B) is sufficient to invert the band ordering around the L energy gap and to induce a topological phase. (C) Radial cuts along different directions from the SD of the strain magnitude. Inset shows the (001) direction of the Burger’s vector piercing the ABC stacked Bi bilayers. (D) Ab initio calculation of the band structure evolution around the L gap under uniform strain. A topological phase transition occurs at about 3.5% strain.

Our ab initio calculations for the response of Bi to a uniform strain (Fig. 4D) show that strain can also promote band inversion at the L TRIM point (36, 54, 55) that would substantiate a strong TI, ν0 = 1, phase with nonzero weak indices, ν = (111), in ab initio calculation (see fig. S5 for more details). The 3 to 4% strain needed for inducing the topological phase transition in calculation is a small fraction of the shear strain applied at the vicinity of the SD. This exemplifies the ability of a bulk topological defect to embed around it a cylindrical volume of an induced topological phase, as marked by the dashed line in Fig. 4B. Nevertheless, in our spectroscopic mappings, we do not identify any signature of gap closing and reopening away from the SD axis, which would necessarily accompany such a shear-induced transition. We thus conclude that Bi is uniformly classified as a strong TI with nonzero weak indices. We also note that recent experiments on single crystals of Bi seem to identify valley symmetry breaking due to spurious strain fields in the samples (22), which alone may promote the Z2 topological classification in Bi.


We argue that the observed phenomenology on all boundaries of Bi including surfaces, step edges, and SDs is strongly influenced by the small energy gap at the L TRIM point. This small energy gap positions Bi in the vicinity of a topological phase transition between a higher-order TI and a strong TI with nonvanishing weak indices. The surface and edge states that disperse between the corresponding projections of the bulk T and L TRIM points on the Z2 trivial side of the phase transition are the precursors of the helical boundaries modes that would be firmly established beyond the critical transition point (for instance, under application of strain or Sb doping). Doping Bi with Sb was shown to invert the bulk energy gap at L while hardly affecting the surface state dispersion around the Γ¯ point (27). Hence, those edge modes are approximately helical and protected from backscattering already on the Z2 trivial side of the topological phase transition. Their penetration depth, the critical length associated with the phase transition, is expected to diverge as the energy gap narrows (39, 56). Their exact topological nature is, therefore, ill-defined in a finite sample (57). A similar boundary mode evolution was indeed reported across the topological phase transition from the 3D TI BiSe2 to a trivial insulator with the substitution of Se with S (5860). Similar reasoning was also applied to the edge modes of a bilayer silicene (61). In both cases, approximate helical boundary modes were shown to emerge on approaching the topological quantum phase transition from the trivial phase.

In summary, we have mapped the surface and edge spectra of Bi(111). We find surface and edge modes that disperse in energy and momentum between the projections of the bulk T and L TRIM points. The edge mode withstands crystallographic irregularities without showing any signs of backscattering and gapping. In particular, the edge mode seems to bind to bulk SD as expected for TI with nonvanishing weak indices. These observations provide first spectroscopic signature for the existence of helical mode on an SD and, by this, allow us to determine the topological nature of Bi. Our findings disfavor the recent identification of the edge mode as a hinge mode of a higher-order TI class, which, in the absence of protection, would be highly susceptible to hybridization by these structural perturbations. We argue that the small scale of the bulk L energy gap, which in our ab initio calculation amounts to 5 meV, positions Bi within the critical region of a topological phase transition to a strong TI with nonvanishing weak indices. It also facilitates strain-induced transition into the topological phase in the vicinity of the SD, as well as throughout the crystal by spurious crystallographic defects. This demonstrates that SD can serve not only as a probe for the topological classification of materials but also as a tool to manipulate and alter this classification.


Sample growth

Single crystals of GdPtBi were grown by the solution growth method from Bi flux in Ta tube. The extra Bi flux was separated from the crystals by decanting the Ta tube well above the melting point of Bi. The complete procedure of crystal growth is described in (50). In this procedure, most of the crystals were found at the bottom of the crucible where they stick. The grown crystals were removed with the help of a stainless steel knife after cutting the tube for further characterizations. After a careful study, we realized that Bi inclusions are invariably present in the 10th of a millimeter of the bottom of the GdPtBi crystals.

STM measurements

The STM measurements were carried out using commercial Pt-Ir tips. The tips were characterized on a freshly prepared clean Cu(111) single crystal. This ensured a stable tip with reproducible results. The dI/dV measurements were performed using the standard lock-in technique. The map parameters for Fig. 1F are V = 200 mV, I = 900 pA, f = 413 Hz, and ac root mean square (rms) amplitude = 1 mV. The map parameters for Fig. 2C and Fig. 3 (B to D) are V = 200 mV, I = 500 pA, f = 1373 Hz, and ac rms amplitude = 3 mV.


Supplementary material for this article is available at

Section S1. Effect of substrate magnetic order on Bi

Section S2. dI/dV mapping of the SD

Section S3. Edge states of Bi nanoribbon

Section S4. Model of shear strain

Fig. S1. Effect of the magnetic order of the GdPtBi substrate on Bi.

Fig. S2. QPI of 2D surface states from SD.

Fig. S3. Direct visualization of the topological edge mode localized on the SD.

Fig. S4. The model used to simulate Bi step edge.

Fig. S5. Effect of shear strain on the electronic properties of Bi.

Reference (62)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge valuable discussions with P. W. Brouwer, L. Trifunovic, B. Andrei Bernevig, R. Ilan, and J. Cano. Funding: H.B. acknowledges support from the United States–Israel Binational Science Foundation (BSF; grant number 2016389), the Helmsley Charitable Trust (grant number 2018PG-ISL006), the German-Israeli Foundation (GIF; I-1364-303.7/2016), and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 678702). Author contributions: A.K.N. and J.R. acquired and analyzed the data. H.B. and N.A. conceived the experiments. H.B., N.A., and A.K.N. wrote the manuscript with substantial contributions from all authors. B.Y., R.Q., and H.F. modeled the system. C.S. grew the material in C.F.’s group. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article