Introduction

The paradigm of topological phases has permeated much of contemporary condensed matter physics1,2,3,4. This fundamentally new way of classifying systems according to global topological properties rather than a local order parameter yielded a deeper understanding of a host of peculiarly robust phenomena5. At the heart of these phenomena lies the bulk-boundary correspondence, which guarantees robust states localized at the perimeter of the topological materials. These boundary states, in turn, might be used as tools for measuring fundamental constants6, as components in thermoelectrics7 or in spintronics devices8.

Time reversal symmetric band insulators can be characterized by a \({{\mathbb{Z}}}_{2}\) index and can be classified into three phases in three dimensions. A normal, or trivial insulating phase, a weak topological insulating (WTI) and strong topological insulating (STI) phase9. These phases are separated by a metallic Dirac or Weyl semimetal phase10. In order to change the \({{\mathbb{Z}}}_{2}\) it is necessary to close and reopen the band gap through one of these metallic phases5. The transition-metal pentatelluride ZrTe5 is an excellent material to investigate topological phase transitions because it lies close to the boundary between a STI and WTI11. Additionally, the material has been widely studied due to its numerous exotic properties. In the monolayer limit it is predicted to be large band gap quantum spin Hall insulator12, it has high-mobility Dirac carriers13,14, it shows the chiral magnetic effect15,16 and the 3D quantum Hall effect17, as well as multiple superconducting phases have been discovered under high pressure18. Recently, the role of spin-dependent Berry phase of the bands had been traced under mechanical strain19.

The topological nature of the bulk ZrTe5 has not been unambiguously identified, with some experiments pointing to a STI and others to a WTI or a semimetal phase, as reviewed by Monserrat et al.20. For example, angle resolved photoemission studies show evidence of a STI phase21, as well as the WTI case22,23. These experiments are supported by first-principles calculations that also show the very same pattern of contrasting results.

In our contribution, we start from ab initio calculations of ZrTe5 and calculate the complete elastic response of the crystal by obtaining the elastic tensor elements, using a new approach compared to the literature24. We validate our ab initio method by reproducing the closing of the band gap at the perimeter of bubbles formed by few layer ZrTe5 on Au(111). Measuring the surface charge density of bubbles provides an almost ideal experimental platform to validate our calculations, because bubbles of few layer van der Waals materials provide a varied deformation landscape25,26,27,28. This deformation landscape leads to a local perturbation of the charge density that can be directly mapped, using a scanning tunneling microscope (STM). As STM directly measures the surface charge density, it does not need fitting or modeling for interpretation. Our calculations reproduce the measured closing of the band gap at the bubble perimeter, without the necessity of fitting phenomenological model parameters. Thus, we provide a robust ab initio method to describe ZrTe5, which is validated by STM measurements. Furthermore, the calculated equilibrium lattice parameters are also in agreement with values measured by X-ray diffraction of our sample. This allows us to establish the electronic ground state of the system and map the topological phase diagram of ZrTe5, revealing an STI phase at equilibrium. Extended data of our calculations and an interactive way to browse them are available online at tajkov.ek-cer.hu/zrte5phasediagram/.

Results

The bulk zirconium pentatelluride (ZrTe5) crystallizes in the layered orthorhombic crystal structure with space group Cmcm (D\({}_{2h}^{17}\)). Crossed trigonal prismatic ZrTe3 chains run along the crystallographic a direction and they are linked along the c-axis via parallel zig-zag chains of Te atoms, as it can be seen on Fig. 1, panel (a). The chains form two-dimensional sheets, stacked along the y-axis, forming a layered structure in the x−z plane. The corresponding unit cell consists of 2 Zr atoms and 10 Te atoms. This is presented on Fig. 1 panel (b), where the atomic resolution STM image is matched with the computer generated image of the crystal in the x−z plane. Each ZrTe5 layer is nominally charge neutral, and the coupling between the layers is of van der Waals type12.

Fig. 1: Geometry of ZrTe5.
figure 1

a The considered unit cell of the crystal. The wavelike Te–Te chains run parallel to the crystallographic c direction and to the z-axis in the corresponding Cartesian coordinate system. b Atomic resolution STM image of ZrTe5 surface. The crystal structure, from the (110) point of view, in the x−z plane is superimposed in a fading manner. The STM image calculated by DFT is overlaid on the crystal structure (right side). It highlights the Te atoms on the sample surface, which dominate the atomic resolution STM topography. STM image measured at 800 pA tunneling current, 400 mV sample bias, temperature 9 K.

Full size image

The initial geometry was fully relaxed by ab initio calculations (see section Methods for more details). The starting point of the relaxation was the lattice constants provided by our X-ray diffraction measurements of the bulk sample. The experimental and relaxed theoretical geometry parameters can be found in Table 1.

Table 1 Unit cell dimensions and positional parameters in fractional coordinates for ZrTe5 as derived from the analysis of X-ray diffraction data at 300 K (calculated standard deviation in parentheses) and from relaxed DFT calculations.
Full size table

As ZrTe5 decomposes under ambient conditions, we exfoliate the bulk crystals onto a Au(111) surface inside an inert glovebox environment, using the scotch tape method. The samples are transferred by a vacuum suitcase to the ultra-high vacuum (UHV) chamber of our STM, allowing the sample surface to remain pristine, as shown by the atomic resolution image in Fig. 1b. Bubbles form during exfoliation of van der Waals materials onto substrates and are predicted to contain inert hydrocarbon contamination29. The (a) and (b) panels of Fig. 2 show the STM topography image of two such bubbles. They have an ellipsoid form and their geometry can be parameterized by the two semi-axes (R1 and R2), as well as the height (h) (see Table 2). Next to the topography images in panel (c) and (d), we present the corresponding measurements of the dI/dV signal within the band gap, which is proportional to the local density of states (LDOS). The middle of the gap is identified as the minimum of the dI/dV signal, as shown in Fig. 2f (for more details see Supplementary Notes 3). Focusing on Bubble1 in panel (c), a halo can be observed at the perimeter of the bubble, which is more intense perpendicular to the x direction (point \({{{{\mathcal{A}}}}}_{1}\)) and absent across the z direction (\({{{{\mathcal{B}}}}}_{1}\)). This indicates an area with increased density of states within the gap, relative to the unstrained areas outside the bubble. The anisotropy of the ring is a consequence of the highly anisotropic nature of the material, as we show below. The increase in the density of states is only observed at energies within the gap. This is illustrated by the fact that no halo is observed when mapping the dI/dV signal at –0.5 V, away from the band gap well within the valence band (see Fig. 2e). For Bubble2 hosted by a thicker flake, shown in panels (b) and (d), the dI/dV map displays no gap closing (point \({{{{\mathcal{A}}}}}_{2}\) in Fig. 2g), because the deformation values are much smaller in the thicker flake.

Fig. 2: Strain induced topological phase transition in ZrTe5 bubbles.
figure 2

a STM topography image of Bubble1 and b of Bubble2. c Map of the dI/dV signal in the same area as in a), measured within the gap. Bright dots on the dI/dV image are point defects. d Map of the dI/dV signal in the same area as in b of Bubble2. e Map of the dI/dV signal in the same area as in panel c but at –500 mV sample bias, deep in the valence band. f dI/dV spectrum measured on the sample surface far away from Bubble1. Colored arrows highlight the sample bias used to measure the dI/dV image in panel c and e. The minimum of the dI/dV signal corresponds to the gap and has a value of 78 pS. g Band gap of ZrTe5 resulting from ab initio calculations, along the path depicted in c and d. The gap, and as a consequence the LDOS, is modulated by the locally varying strain within the bubbles (see Fig. 3). STM measurement parameters: 300 K, 500 pA, sample bias for a–d –100 mV and –500 mV for panel e.

Full size image
Table 2 Geometric parameters of Bubble1,2.
Full size table

It is worth noting that the reduced LDOS within the bubble, stems from the lack of direct contact with the substrate. The density of states of Au is orders of magnitude larger than that of the semimetallic ZrTe5 and increases the LDOS measured in areas where the two materials are in close contact30,31. This increase in LDOS is also present when measuring within the valence band (Fig. 2e) but for thicker flakes it becomes much reduced (Fig. 2d).

We calculated the strain pattern of the bubbles surfaces, using a unique technique that combines finite element calculations and density functional theory. After we obtain the deformation of the bubbles, we can calculate the local electronic structure and the size of the gap along the path that is depicted in Fig. 2 panel (c) and (d) of the distorted crystal, using density functional theory. Throughout the rest of this article, unless otherwise stated, gap is taken to mean global gap, especially since this is what is reflected in the STM measurements. As it can be seen in panel Fig. 2g, DFT predicts a gap closure at the edge of Bubble1, but only at the \({{{{\mathcal{A}}}}}_{1}\) point and not at \({{{{\mathcal{B}}}}}_{1}\). The absence of the band gap means that the density of states must be higher in the area around the point, in agreement with the measurements. Tracing the same path along the surface of Bubble2, no gap closure can be seen, in good agreement with the measurement.

In order to avoid computationally intensive simulations for the bubbles containing a large number of atoms, we combined density functional theory and finite element calculations. In this method, the mechanical stress field of the bubbles was determined using the finite element method, with the bubbles being considered as continuous elastic materials. The obtained stress values were then used to perform DFT electronic structure calculations in bulk samples, with mechanical distortions corresponding to the local stress values. To apply this methodology, first, we must establish a proper description of the elasticity tensor of the material. For this we calculated the stiffness tensor by fitting the free energy of the distorted crystal via DFT. For more details see Methods section. The corresponding tensor elements are depicted in Table 3 in Voigt notation32. This stiffness tensor can then be used to describe the strain pattern of the bubbles, using finite element calculations. We simulated the bubbles as a continuous anisotropic material, which can be characterized by the calculated stiffness tensor. We matched the thickness of the sample to the experiments and applied hydrostatic pressure25 to the bottom of the system until we reached the measured height and reproduce the experimentally measured height profile. For more information about the procedure see Supplementary Notes 2. The calculated strain patterns can be seen in Fig. 3. In panel (a), we show a schematic representation of a bubble from a side view and denote the most important geometrical and strain parameters. The yellow, diagonally dashed region denotes the substrate. The solid black line in the middle of the bubble shows the neutral plain. The different trapezoid quadrilaterals indicate the different strain patterns. In Fig. 3 panel (b), we show the strain pattern for the two bubbles, plotting all relevant strain tensor elements on the surface of the bubble. We have only omitted the presentation of the εいぷしろんxz component as it is practically zero everywhere in the sample. We show the exact numerical results in the supplementary material. The first row corresponds to Bubble1. It can be observed that the strain patterns in the two bubbles show the same qualitative tendency, but the magnitudes are almost 5 times larger in the first bubble, due to the larger aspect ratio25.

Table 3 Stiffness tensor obtained by DFT in Voigt notation and in the units of GPa.
Full size table
Fig. 3: Numerical results.
figure 3

a A simplified view of the strain patterns in a bubble. The yellow, black dashed region indicates the gold substrate. The middle black solid line depicts the neutral plane, which shares the same geometrical parameters as the unstrained crystal. The different trapezoid quadrilaterals indicate the different strain patterns. We have also indicated the geometric parameters. b Strain pattern on the top surface of the bubbles obtained by finite element calculations (COMSOL). This is the surface measured by STM. The top row corresponds to Bubble1, the bottom row to Bubble2. The different panels show different strain tensor components. Negative values correspond to compression. c The influence of different strain components on the electronic properties. The magnitude of the gap as the function of the strain components, calculated by DFT.

Full size image

As a final step using the strain field we obtained from the finite element method, we can calculate the band structure of the system under the influence of the strain patterns by using DFT. In panel (c), we present the effect of four different strain components on the band gap of the crystal. In the left side of subfigure c, we show the elements that are responsible for the in-plane shear. It is clear that both decrease the band gap, but the εいぷしろんxy component has a smaller contribution compared to εいぷしろんyz. Furthermore, it can be clearly observed that the magnitude required to close the gap is as high as 7.5%, which is too large to explain the effect observed on Bubble1. We can find the explanation to the LDOS halo observed in the meauseremnts by looking at the curves in the adjacent panel. The solid blue line denotes the effect of the εいぷしろんzz component on the band gap. Positive strain values close the gap at around 1.5%, but the negative values first widen the gap and after 2% start closing it, eventually closing at around 4%. We found that εいぷしろんxx and εいぷしろんzz behave identically for positive values, only showing a drastic difference for compression, namely that εいぷしろんxx closes the gap for much reduced strain values.

This allows us to explain the increase in LDOS at the Bubble1 perimeter as a closing of the gap. We consider the calculated gap values along the path in Fig. 2 panel (g). Starting from the equilibrium Ox point along the x-axis and moving towards the \({{{{\mathcal{A}}}}}_{1}\) point, the εいぷしろんxx strain element becomes a large enough negative number to close the gap. This magnitude then goes to zero, elevates again, becomes a relatively smaller positive value, but never reaches the magnitude that closes the gap again. After we reach the center of the bubble (point C) and turn upwards along the z-axis the εいぷしろんzz becomes smaller, reaches zero, and then becomes a large negative value. However, as εいぷしろんzz has a different role in the manipulation of the electronic structure of the crystal it slightly widens the gap in the \({{{{\mathcal{B}}}}}_{1}\) points compared to the Oz points, where it reaches equilibrium again. As for the second bubble, the explanation is much simpler, the magnitude of the strain components never reaches a sufficiently high value to have an observable influence on the gap. In the supplementary material we present the exact values for every element of the strain tensor along the whole path (see Supplementary Notes 3).

Since our ab inito calculation reproduces the effect of the complex strain pattern on the ZrTe5 band gap, we can take one step further and map the phase diagram of the system. The first panel in Fig. 4 shows the contour map of the band gap size (Eg) under different mechanical strains. The horizontal axis indicates the in-plane (the lattice vectors a and b were distorted isotropically) strain from –1.2% to 3.1% in 20 steps, while the vertical axis corresponds to the van der Waals direction (y-axis) from –5% to 10% in 20 steps. At every point a sign has been assigned to the gap as the \({{\mathbb{Z}}}_{2}\) invariants were calculated, positive (negative) sign indicates WTI (STI). The phase diagram shows three main domains. Around the equilibrium the system is a STI (1) and it can be tuned to the WTI (3) phase through a conducting phase (2). The black line in the conducting phase shows where the Dirac cones in the Γがんま point touch each other.

Fig. 4: Topological phase diagram of ZrTe5.
figure 4

a The phase diagram of the electronic structure of the crystal under mechanical strain. At every point the size of the gap was calculated and a sign has been assigned to it according to the topological flavor of the gap. The negative gap corresponds to the strong topological insulating phase, while the positive gap to the weak topological phase. The black solid tentative line indicates the boundary where the Dirac cones in the Γがんま point touch each other. The green dots assigned with a number denote the corresponding band structure in the b subfigure. The inset shows the corresponding Brillouin zone indicating the high-symmetry points. b The calculated band structure along the path of the high-symmetry points. The opaque band shows the size of the gap and color indicates the topological flavor.

Full size image

The conclusion that the equilibrium state is a STI, is supported by our dI/dV measurements. For a STI, the LDOS within the gap should not go down to zero because of the presence of the topological surface state33. Taking into account the noise level in our instrument, a value of zero LDOS would correspond to a tunneling conductance value lower than 1 pS. Compared to this, the conductivity at 300 K inside the gap is 78 pS (see Fig. 2f) and at 9 K is 74 pS (see Supplementary Notes 3.).

Figure 4b shows three typical band structures at points 1, 2 and 3 marked in panel (a), and total density of states in 1/eV units. The first band structure shows a strong topological insulating phase, the red opaque band highlights the 22 meV band gap. As we go towards point (2) we reach the gray area where the band gap is closed but there is no touching of the bands. At point (2) the bands corresponding to the massive Dirac fermions touch each other around Γがんま. As we go further towards point (3) the Dirac cones open up and the band between the high-symmetry points A1 and T lift up from the Fermi level. This opens a gap in the weak topological insulating phase.

The phase diagram, as defined above and presented in Fig. 4a, allows us to explain the different, sometimes contradictory, numerical results on the electrical properties of the material20. It may happen that calculations that respect the symmetry of the material give different results for the ground state geometry and therefore, for example, predict a weak topological insulating state34. We repeated our calculations using another commonly used DFT code, VASP, that resulted in a similar phase diagram to the one presented in Fig. 4. Besides the similarities, there are also differences that are instructive and bear information. A detailed discussion of these calculations and their implications can be found in section Supplementary Notes 4.

Discussion

ZrTe5 is a unique material, with electronic properties that are extremely sensitive to its lattice parameters. This over-sensitivity has led to controversies in the crystal’s literature in recent years. On the other hand, the sensitivity makes it a perfect candidate for strain engineering the material across a topological phase transition. In this contribution we have combined density functional theory and finite element calculations with direct measurements of the charge density, allowing us to identify the mechanical deformation needed to induce topological phase transitions in the material. We showed that shear in the plane perpendicular to the van der Waals direction could also lead to a phase transition. Our results pave the way for a robust understanding of the band structure and topological properties of ZrTe5, enabling future investigations into the physics of topological phase transitions as well as applying these in electronic devices.

Methods

Density functional theory

The optimized geometry and electronic properties of the crystal were obtained by the SIESTA implementation of density functional theory (DFT)35,36,37,38. SIESTA employs norm-conserving pseudopotentials to account for the core electrons and linear combination of atomic orbitals to construct the valence states. The generalized gradient approximation of the exchange and the correlation functional was used with Perdew-Burke-Ernzerhof parametrisation39 and the pseudopotentials optimized by Rivero et al.40 with a double-ζぜーた polarized basis set and a realspace grid defined with an equivalent energy cutoff of 350 Ry for the relaxation phase and 900 Ry for the single-point calculations. The Brillouin zone integration was sampled by a 30 × 30 × 18 Monkhorst-Pack k-grid for both the relaxation and the single-point calculations41. The geometry optimizations were performed until the forces were smaller than 0.1 eV nm−1. The choice of pseudopotentials optimized by Rivero et al. ensures that both the obtained geometrical structures and the electronic band properties are reliable. After the successful self consistent cycles the necessary information was obtained by the sisl tool42. The spin orbit coupling was taken into account in the single-point calculations. In every case we have calculated the \({{\mathbb{Z}}}_{2}\) invariant using a home made tool based on the numerical method developed by Fukui and Hatsugai43,44.

The simulated STM image was obtained for a three layers thick sample, that we cut from the original bulk crystal in the (110) orientation. We used Monkhorst-Pack resolution 20 × 12 × 1 in the geometry relaxation this case, where the samples were separated with a 40 Å thick vacuum in the perpendicular direction. The simulated STM image was obtained by the tools developed by the SIESTA developers, where the position of the tip was 2.5 Å away from the position of the topmost Te atom.

Stiffness tensor

We obtained the stiffness tensor elements by fitting the free energy change of the crystal under mechanical deformations. The change in the free energy of the crystal is a quadratic function of the strain tensor45. The procedure presented here is the most recent way to precisely determine the elements of the stiffness tensor using DFT46.

The general form of a deformed crystal is the following45:

$$F=\frac{1}{2}{\lambda }_{klmn}{u}_{kl}{u}_{mn},$$
(1)

where λらむだklmn is the elastic modulus tensor, uij is the strain tensor. The general expression for the free energy in the orthorhombic system is

$$\begin{array}{l}F=\frac{1}{2}{\lambda }_{xxxx}{u}_{xx}^{2}+\frac{1}{2}{\lambda }_{yyyy}{u}_{yy}^{2}+\frac{1}{2}{\lambda }_{zzzz}{u}_{zz}^{2}+{\lambda }_{xxyy}{u}_{xx}{u}_{yy}+\\\qquad +{\lambda }_{xxzz}{u}_{xx}{u}_{zz}+{\lambda }_{yyzz}{u}_{yy}{u}_{zz}+2{\lambda }_{xyxy}{u}_{xy}^{2}+2{\lambda }_{xzxz}{u}_{xz}^{2}+2{\lambda }_{zyzy}{u}_{zy}^{2}.\end{array}$$
(2)

It contains nine moduli of elasticity45. By applying mechanical strain to the relaxed geometry in the ab initio calculations using the free energy obtained by SIESTA software we can fit Eq. (2) to get the different moduli.

Finite element method

The corresponding strain patterns were calculated by numerically solving the three-dimensional equation of motion by the finite element method (FEM). The FEM calculation was performed by the MEMS Module of the COMSOL Multiphysicsⓒ47 5.6 software package. The bubbles were simulated by a 20,000 × 20,000 Å2 block. The thickness of the block was matched to the sample size in the STM measurements. The sides and the bottom of the block were fixed and we applied hydrostatic pressure on the bottom of the block on an ellipse shaped part of the bottom to match the shape of the bubbles in the STM measurements. The pressure was chosen to match the height of the bubbles.

Sample preparation and STM measurements

The ZrTe5 crystals were purchased from hqgraphene.com and exfoliated onto gold substrates48, inside an inert glovebox environment. The samples were transferred to the chamber of the UHV STM, via a vacuum shuttle. Atomic resolution images show that the ZrTe5 crystal surface remains pristine after transfer into the STM chamber (see Fig. 1).

STM measurements were performed using an RHK PanScan Freedom microscope at 300 K and 9 K temperatures in UHV, at a base pressure of 5 × 10−11 Torr. STM tips were prepared by mechanically cutting Pt/Ir (90%/10%) wire. We have used a large working distance optical microscope to place the STM in the vicinity of selected ZrTe5 crystals. dI/dV spectra were measured using a Lock-In amplifier, with a reference frequency of 1.372 kHzきろへるつ and a bias modulation of 5 mV at 9 K and 30 mV at 300 K.

X-ray diffraction

X-ray diffractometry measurements were performed in the parallel geometry using a Bruker AXS D8 Discover diffractometer equipped with Göbel-mirror and a scintillation detector with Cu Kαあるふぁ radiation. The X-ray beam dimensions were 1 × 5 mm, the 2 step size was 0.02∘, scan speed 6 s/step. We used the Diffrac.EVA program and the ICDD PDF database for phase identification.