Engineering topological phases in triple HgTe/CdTe quantum wells.

Quantum wells formed by layers of HgTe between Hg[Formula: see text]Cd[Formula: see text]Te barriers lead to two-dimensional (2D) topological insulators, as predicted by the BHZ model. Here, we theoretically and experimentally investigate the characteristics of triple HgTe quantum wells. We describe such heterostructure with a three dimensional [Formula: see text] Kane model, and use its eigenstates to derive an effective 2D Hamiltonian for the system. From these we obtain a phase diagram as a function of the well and barrier widths and we identify the different topological phases composed by zero, one, two, and three sets of edge states hybridized along the quantum wells. The phase transitions are characterized by a change of the spin Chern numbers and their corresponding band inversions. Complementary, transport measurements are experimentally investigated on a sample close to the transition line between the phases with one and two sets of edges states. Accordingly, for this sample we predict a gapless spectrum with low energy bulk conduction subbands given by one parabolic and one Dirac subband, and with edge states immersed in the bulk valence subbands. Consequently, we show that under these conditions, local and non-local transport measurements are inconclusive to characterize a sole edge state conductivity due to bulk conductivity. On the other hand, Shubnikov-de Haas (SdH) oscillations show an excellent agreement with our theory. Particularly, we show that the measured SdH oscillation frequencies agrees with our model and show clear signatures of the coexistence of a parabolic and Dirac subbands.

Introduction

The discovery of two and three dimensional (2D and 3D) topological insulators (TIs), also known as quantum spin hall (QSH) insulators, strongly impacted the field of quantum materials due to their interesting fundamental properties and technological applications[1-8]. They constitute a peculiar class of materials, characterized by an insulating bulk dispersion and gapless topological helical surface or edge states that are shown to be protected against back-scattering by the time reversal symmetry. The first theoretical prediction for a TI was proposed by Haldane[9], and it is built upon a “2D graphite” spinless toy model. Although this proposal lacked physical justifications at that time, it was later theoretically realized in spinful graphene[2] despite its spin-orbit gap being too small[10] to make them experimentally realizable. In fact, whereas a broad variety of semiconductor-based materials can host topological helical states[11], the first experimental indication of edge channels conductivity was reported in HgTe/CdTe composite quantum wells (QWs)[12,13], following the theoretical prediction of the BHZ model[3,4] with inverted subbands (named after its authors Bernevig, Hughes, and Zhang). Additionally, the 2D Dirac-like band structure of surface states in 3D TIs have been observed in ARPES measurements[14-19].

Particularly in 2D TIs, the energy ordering of electron-like and hole-like QW eigenstates leads to the topological phase transition from a trivial insulator to a TI. This is shown to be controlled by the QW width[3,4], electric field[20-25], strain[26,27], temperature[26,28], superlattice potential[29] and concentration of Bi doping[30]. The resulting 1D helical channels at the edges lead to quantized conductance and nonlocal edge transport[12,13,31], which has been observed for sufficiently short distances between measurement probes[32]. In all cases, the quantized conductance of HgTe-based QWs shows significant fluctuations, with experimental values different from the ones predicted theoretically through Landauer-Buttiker formalism. The deviation from the theoretical prediction has been attributed to many different effects, including disorder[33,34], charge puddles[35,36], and many sources of inelastic scattering[37]. Despite the deviation from theoretical results, recent improvements in sample growth has lead to measurements closer to the predicted quantization[36,38].

Recently, it has been proposed that the control of the layer localization of topological states in bilayer graphene[39,40] and TIs[41] could be used to design “layertronic” devices, where this additional degree of freedom could be used together with the spin, valley, and charge to build novel devices. For HgTe-based 2D TIs, the multilayer character arises from multiple coupled QWs. In the bilayer case, i.e. double HgTe QWs (DQWs), an intuitive 2D model for coupled QWs has been introduced by Refs.[42,43], with Ref.[23] showing that a variety of topological phases can be obtained depending upon the QW geometrical parameters. Additionally, it has been recently proposed[44] that these DQWs could host second-order topological insulators with excitonic nodal phases that support flat band edge states, which could lead to superconductivity[45]. Experimentally, signatures of the conductance quantization, nonlocal transport and the Landau fan diagram of HgTe-based DQWs have been recently observed[46,47]. Additionally, magnetically doped TI layers, coupled through insulating spacers, lead to a solid state realization of the Weyl semimetal[48], providing a platform to examine interesting topological features such as Fermi arc surface states and the chiral anomaly effect[5]. The Weyl semimetal can also be achieved in time-reversal invariant systems with broken inversion symmetry[49,50], which can be realized in multiple coupled HgTe-based QWs[51,52]. Further advance on the physics of the multilayer Dirac fermions can be achieved in trilayer systems, i.e. triple quantum wells (TQW). In contrast to the DQW case, the additional layer of the TQW allows for an interplay between the hybridization of the inner and outer layers, which can be controlled by its geometric parameters (wells and barrier widths), external electric fields and gates.

In this paper, we investigate the band structure and transport of triple HgTe/HgCdTe quantum wells, comparing theoretical transport predictions with experimental measurements of local and nonlocal resistivities, and Shubnikov-de Haas (SdH) oscillations. First, we obtain the band structures calculated with the Kane Hamiltonian, and determine the corresponding topological phase diagram as a function of the geometric TQW parameters (well and barrier widths). We find that electron-like (E) and hole-like (H) eigenstates are hybridized throughout the three QWs as symmetric and anti-symmetric combinations. For the E-like subbands the hybridization leads to large energy splitting between these subbands, while the H-like subbands remain nearly degenerate due to its larger effective mass. Additionally, we show that by increasing the QW widths, only crossings between E-like and H-like subbands with equivalent envelope wavefunction profiles along the wells lead to topological phase transitions. Using their spin Chern number, we label four different phases of our system, namely, 0 (trivial), and I, II, III, referring to the number of pairs of edge states in each phase. Particularly, for phases I and II we show that bulk band structure is gapless and the edge states always coexist with the valence bands. Only for the phase III a bulk gap opens and the three edge states become isolated from the bulk within the gap region.

Using the theoretical prediction of the topological phase diagram corresponding to our HgTe-TQW, we are able to predict that the experimental sample analyzed in this work lies in phase I, close to the transition towards phase II. Within this regime, we show that resistivity measurements do present signatures of nonlocal transport. However, these are strongly affected by bulk conductivity and scattering between different states, which leads to an inconclusive characterization of the total number of edge states conducting the current. On the other hand, the measurements of SdH oscillations shows a good agreement with the theoretically predicted SdH frequencies and also with the temperature dependence of the SdH oscillations. More specifically, we explained the SdH oscillation data as arising from an interplay of one linear and one parabolic conduction subbands, which are shown to be in agreement with our band structure calculation. This analysis is shown to be also in agreement with dependence of the peak nominal values and their spectral weight as a function of the Fermi energy.

Theoretical models

In this section we introduce the theoretical models used to obtain the bulk energy bands, edge state energies, and all the corresponding wave-functions of our system.

Kane model

We consider quantum wells (QWs) made of HgTe confined by HgCdTe barriers with concentration , as shown in Fig. 1a. Both HgTe and CdTe crystallize in the zincblende structure with low energy bands around the point (), which are well described by the Kane Hamiltonians[53,54] , with the corresponding {, , , } basis set. The CdTe bandstructure has a normal order, where is a S-type conduction band, and are the P-type valence bands corresponding to heavy and light holes () and the split-off band (). In contrast, HgTe has the and bands inverted due to relativistic fine structure corrections (Darwin, spin-orbit, and mass-velocity terms), which ultimately allows for the QSH topological phase of single HgTe QWs[3,4]. For heterostructures, one considers the Kane parameters to be position dependent, restores as the momentum operator, and symmetrizes[53-55] the Hamiltonian. Additionally, the growth of the heterostructure typically induces strain, which is considered under the Bir-Pikus Hamiltonian[56]. The resulting Kane Hamiltonian and the material parameters are shown in the Supplementary Material. Here, the theoretical model is set with the growth direction . Nevertheless, for small  nm we expect the results to be nearly equivalent for growth directions [001] and [013][57]. In the (x, y) plane, the solutions are given by plane-waves, , where is the in-plane momentum. To numerically diagonalize the we use the kwant code[58], which provides an efficient interface to build and solve the numerical problem.

Figure 1

(a) Conduction and valence band edges of the triple quantum well schematically shown in the bottom. The widths of the HgTe wells and the thickness t of the HgCdTe barriers () are indicated. Single well states located on each of left (L), central (C) and right (R) wells are illustrated in purple (electron-like) and blue (hole-like). As t decreases their hybridization leads to a large energy splitting of the E-like states at , while the H-like remain nearly degenerate due to the larger effective mass. (b,c,d) Band structures for  nm and varying . (b) For  nm, in the trivial regime, all E subbands are above the H subbands. indicates the indirect band gap. Phase transitions occur as each E subband crosses the H subbands with increasing (c)  nm, (d)  nm, and (e)  nm. For large the indirect gap closes and the system becomes semi-metallic (shown in the Supplementary information which contains the Kane Hamiltonian[53–56], the table of material parameters, the expressions for the matrix elements of the 2D model[42,43], and extra figures complementing the ones in the main text). All of the plotted subbands above present degeneracy in spin.

Effective 2D Hamiltonian for triple wells

To investigate the confinement of the subbands of our triple HgTe/CdTe quantum wells, their corresponding topological character and the characteristics of their edge states, we consider an effective 2D Hamiltonian for our system. For single HgTe-quantum wells, this is achieved by the projection of the Hamiltonian into its eigenstates, which leads to the well known BHZ model[3,4]. In contrast, for double HgTe-quantum wells there are two interesting approaches. First, similarly to the derivation of the BHZ Hamiltonian, in Refs.[23,59] the authors project the total Hamiltonian into the DQW eigenstates. Alternatively, in Ref.[42] the authors project the total DQW Hamiltonian into the subbands of the single wells (left and right), and introduce tunneling parameters for the coupling between neighboring QWs. Here, for the case of triple QWs, we follow the approach from the latter, as it provides an intuitive perturbative picture of the coupling between quantum wells. This is illustrated schematically by Fig. 1a). For easy reading, we keep here a notation similar to Ref.[42], but we introduce the index to label the quantities of the individual left (L), central (C) and right (R) QWs. Accordingly, we define the subbands of the isolated QW as , with ± labels referring to time-reversal partners. Assuming that tunnel coupling occurs only between neighboring QWs, it is immediate to extend the double QWs model[42] into the triple well case, which we label “3BHZ”, and it reads as

Here, the diagonal blocks of each layer, , are given by a direct sum over Kramers partners , each composed by BHZ-like Hamiltonianswith the Pauli matrices (, , and ) acting on the subspace of each Kramers block. Similarly, the tunnel couplings for , referring to the left-central and central-right QW couplings, are , with

All coefficients above are calculated following the perturbation theory up to second order. These are shown in the Supplementary information.

Theoretical results

In this section, we discuss the energy dispersions of the TQWs. The different topological phases of the TQW are presented in terms of a phase diagram as a function of the geometric parameters and t (see Fig. 1a), and labeled by the corresponding spin Chern number. The edge state dispersions are illustrated for representative cases of each phase.

Energy subbands

Figure 1b–e illustrate the band structure for triple HgTe QWs with  nm and increasing  nm. In Fig. 1b the system is in the trivial regime i.e., all three conduction subbands have a predominantly electron-like (E-like) character, while the three valence subbands have a predominantly hole-like (H-like) character. Here, we have, for each Kramers pair, three non-degenerate conduction E-like subbands (, , and ), which arise from the hybridization of the lowest conduction subbands of the left (L), central (C) and right (R) wells, namely, , and . On the other hand, the three H-like subbands (for each Kramers pair) are nearly degenerate because the heavy-hole states are strongly localized within each well due to their larger effective mass, and only show significant hybridization for  nm. As we increase in Fig. 1b–e, the E-like subbands cross the H-like subbands, one by one. Each time a E-like subband crosses down the three H-like subbands (with same Kramer pair index), one H-like subband flips the sign of its effective mass, but it still remains nearly degenerate with the other H-like subbands at due to the small overlap of their wavefunctions. Similarly to the case of only one well, the subband inversions produce topological phase transitions. Here, however, only some of these crossings characterize the phase transitions, and this will be discussed with more detail in the next section. For now, we should only note that the band structures in Fig. 1c,d are topologically non-trivial, but gapless. Moreover, it is only when all the E-like subbands are above (Fig. 1b) or below (Fig. 1e) all the H subbands that the system shows a well-defined gap, and as a consequence can be claimed to be either a trivial or topological insulator.

It is also important to stress that the valence subbands obtained here also show the “camel back” profile[59,60], thus also presenting an indirect gap defined by (see Fig. 1b–e). This feature appears due to the strong hybridization of the QW subbands for large thickness  nm (or small  nm), where the subbands are close (in energy) to each other (see Fig. 2a). Furthermore, for large  nm (see Supplementary information), the indirect gap (see Fig. 1b) closes and the system becomes semi-metallic (SM). Notice that in Fig. 1e the indirect gap is already smaller than the direct gap. We emphasize that throughout this work, the notation “semi-metalic (SM)” will refer to systems in which the indirect gap vanishes, while the “gapless” will refer to band structures with vanishing gap at .

Figure 2

(a) Crossings between the E-like subband edges and the H-like subbands as a function of , and  nm (All the subbands are spin-degenerate). (b) Phase diagram as a function of and t. The black solid lines mark the E–H subband crossings, and the red circles along the  nm line correspond to those in (a). Shaded region marks the SM phase, and the colors within phases I and II refer to hybridization gap between H subbands, , which becomes clear only for  nm, and also splits the first solid black line for the – crossings. The cyan stars along the  nm line mark the points (0, I, II, III) corresponding to the phases illustrated in Fig. 1b–e. The shaded rectangle near  nm and  nm marks the parameters of the experimental sample, with the area shaped to illustrate the experimental uncertainty of  nm in both and t.

Topological phase transition and Chern number

In principle, within our system we have three conduction subbands crossing three different valence subbands, yielding a total of nine different inversions. However, only three out of those nine inversions give rise to a topological phase transition, and therefore, a precise characterization of the inversions becomes important. For instance, a counter intuitive scenario occurs in InAs/GaSb type-II QWs[61], where the topological phase transition takes place as E-like states localized at the InAs layer crosses H-like states from the GaSb layer.

It is important to stress that a band inversion is a necessary, but not sufficient ingredient to have a topological phase transition. Accordingly, bands must not only invert, but also hybridize to open a gap between them after the inversion. It turns out that only three out of the nine inversions in the TQW satisfy these conditions. To identify which are the relevant crossings, we diagonalize the effective Hamiltonian from Eq. (1) at . Using to label the E-like and H-like subbands, the diagonal subbands for the case of three identical QWs read as

Projecting the Hamiltonian from Eq. (1) into this basis yields

Here, the block diagonal terms with have the usual BHZ form with , andwith renormalized parameters , and . It is evident from Eq. (7) that the only hybridization happens between pairs , and . Accordingly, it follows from the BHZ model[3,4] that the topological phase transition will only take place when the energies of these individual pairs invert, i.e. as each changes from positive to negative, the corresponding spin Chern number goes from 0 (trivial) to 1 (topological).

To understand and characterize how the QW hybridizations lead to phase transitions, in Fig. 2 we draw the corresponding phase diagram as a function of the well width and barrier thickness t. First, for fixed  nm in Fig. 2a, we see that as the well width increases, the E-like (H-like) subbands move down (up) in energy, a feature well known for the single HgTe/CdTe quantum wells[3,4]. The crossings between the E-like and H-like subbands pairs, , are highlighted (red circles) in Fig. 2a and also in Fig. 2b, along the  nm line. In Fig. 2b, the solid black lines mark the parameter values that yield crossings between subband pairs . These lines delimit different regions of our phase diagram, which are labeled by their corresponding spin Chern number, i.e., 0, I, II and III. This can be clearly seen in Fig. 1b–e, which contain the bandstructures corresponding to the cyan stars in Fig. 2b. As a consequence, the region 0 corresponds to the trivial insulator regime, while regions I, II, and III correspond to the topological insulator regimes with one, two and three pairs of topological helical edge states, respectively, despite the gapless character of the full spectrum of both phases I and II for  nm. Finally, the shaded region marks the semi-metallic phase, representing the cases where the indirect band gap closes. The color map represents the nominal value of the gap at , which only becomes significant for  nm.

Edge states

To illustrate the characteristics of each phase presented above and in Fig. 2b, we now plot and analyze the energy spectrum of representative cases in the presence of an extra confinement along the y direction. To calculate the spectrum for the confined system, we consider the effective 3BHZ 2D model presented above, which describes the effective Hamiltonian for the three lowest (highest) conduction (valence) subbands. Additionally, here we consider a hard-wall confinement at , with  nm. The spectrum is then calculated using a recursive Green’s function method[62], which allow us to calculate the local spectral function, , for the bulk and edges states with an efficient exponential decimation. In Fig. 3a1,b1,c1 we compare the band dispersions calculated with the Kane model (dashed lines) and with the effective 3BHZ 2D model (solid lines) for the same parameters of Fig. 1c–e. We see that, for  nm, these bands are in good agreement, although the effective 3BHZ 2D model is not able to reproduce accurately the “camel back” profile around  nm. In Fig. 3a2,b2,c2 we have the local spectral function corresponding to the cyan stars, I, II and III in Fig. 2b, with gray colors representing bulk subbands, and orange and green representing states localized at edges of our system. Additionally, Fig. 3a3,b3 show a zoom of the data from the corresponding Fig. 3a2,b2 for clarity. Interestingly, it is possible to identify here two different types of edge states in Fig. 3a2,b2 for phases I and II. The ones indicated by the orange color are topological edge states connecting conduction to valence subbands. Conversely, the ones indicated by the green color are trivial edge states predicted previously in Ref.[63]. While the topological ones arise from the non-trivial topology of our system, the trivial ones appear when we confine subbands that have a strong linear dispersion[63]. For this reason, these edge states appear due to the approximately chiral symmetry of the subbands[63]. Even though the trivial edge states are not protected against backscattering, it was shown that it is possible to make these edge states protected when the ribbon is reduced to a quantum dot[30].

Figure 3

(a1,b1,c1) Comparison of the bulk band structures calculated with the Kane model (black dashed lines) and the effective 2D 3BHZ model (solid lines) with E-like subbands in purple and H-like subbands in blue. The panels correspond to the phases I, II and III indicated by the cyan stars in Fig. 2b. (a2,b2,c2) Band structures for a nanoribbon geometry of width  nm shown by the local density of states of the bulk (gray tones) and edges (orange and green tones). The E-k range match the shaded regions from the panels above. The 3BHZ bulk subbands are shown as the solid lines as a guide to the eyes. (a2) For  nm the topological regime with one E-subband below the H subbands shows one pair of topological edge state (orange) and two pairs of trivial edge states (green). (b2) As a second E-subband crosses the H subbands for  nm one gets two topological edge states a single trivial edge state. (c2) For  nm the system reaches the full topological regime with three topological edge states and a well defined gap. (a3,b3) Zoom into the rectangles of (a2,b2) emphasizing the edge states. All of the subbands above are spin-degenerate.

Experimental results

Triple quantum wells based on HgTe/CdHgTe with [013] surface orientation and equal well widths of  nm and barrier thickness  nm were prepared by molecular beam epitaxy (MBE). The sample structure is shown in Fig. 1a. The layer thickness was determined by ellipsometry during MBE growth, with accuracy of  nm. The devices are multiterminal bars containing three 3.2 μm wide consecutive segments of different length (2, 8, and 32 μm) and nine contacts (see inset in Fig. 4b). The contacts were formed by the burning of indium to the surface of the lithographically defined contact pads. The growth temperature was near 180 °C, therefore, the temperature during contacts fabrication process was relatively low. On each contact pad, the indium diffuses vertically down providing an ohmic contact to all the three quantum wells together, with contact resistance in the range of 10–50 k. During AC measurements we continuously checked that the reactive component of the impedance never exceeds of the impedance, which demonstrates the good ohmicity of the contacts. The I-V characteristics are also ohmic for low voltages. A dielectric layer, 200 nm of SiO, was deposited on the sample surface and then covered by a TiAu gate. The density variation with the gate voltage was estimated to be  cm/V from the dielectric thickness, as previously reported in studies using similar devices[46,47], and in comparison with the calculated frequencies of the SdH oscillations shown in the next section. Experimentally, the density variation was measured from the classical Hall effect in the low magnetic field range in agreement with the value above. The transport measurements were performed in the range of temperatures (T) from 1.4 to 80 K by using a standard four-point circuit with 1–13 Hz AC current of 1-10 nA through the sample, which is sufficiently low to avoid overheating effects.

Figure 4

(a) Longitudinal (green) and Hall (orange) resistances at B = 3 T and T = 4.2 K as a function of gate bias. (b) Local resistance as a function of the gate voltage measured along segments with different lengths (A to C) and nonlocal result (D). Insets show device scheme and measurement configurations. The dashed lines are the expected resistances calculated using Landauer-Buttiker formalism for a device with 9 terminals including the contribution of several pairs of edge states. (c) Resistance as a function of the gate voltage for different temperatures. The inset shows the resistance at CNP as a function of 1/T with solid line corresponding to with an activation energy[64]  meV. Shubnikov-de Haas oscillations as function of (d) temperature and (e) gate voltage. Inset in (d) shows the fitting of the amplitude dependence at a fixed magnetic field of 1.2 T for two different values of where the solid line is a fitting using Lifshitz-Kosevich formula. Curves in (e) were vertically shifted for clarity and the dashed line is a guide to the eye for the change in the oscillations phase.

Three devices from the same substrate were studied. Figure 4a shows the measured transport under strong magnetic field that identifies the charge neutrality point (CNP), with near zero carrier density, in the energy spectrum. Sweeping the gate voltage (V) from positive to negative values depopulates the electron states and populates the hole states, while the Fermi level passes through the CNP. The longitudinal resistance exhibits oscillating behavior on the electronic side, however at the hole side the resistance shows monotonic behavior due to strong scattering between the cone and the heavy hole branches[65]. In the CNP, the electron-like Hall resistance jumps from the negative quantized value h/e to the hole-like positive value h/4e. The quantum Hall effect in HgTe TQWs is beyond the scope of this work and will be reported in a forthcoming publication.

Figure 4b shows local and nonlocal resistance as function of gate voltage in a representative device. In the local case, the current flows between contacts 6-1 and the voltage was measured in the different device segments: 3-2 (curve A), 4-3 (curve B), 5-4 (curve C). For the nonlocal case, the current flows between 7-5 and the voltage measured in 8-4. The resistance maximum for all curves occurs at the CNP, as identified in (a). Curve C, situation closer to the ballistic transport, presents a maximum in agreement with the Landauer-Büttiker calculation for two pairs of edge states in a device with nine terminals. Nevertheless, in the next section, we will discuss that this direct association can be misleading if the bulk contributions are not considered.

In order to identify the nature of the transport in the triple QW sample, we have measured the temperature dependence of the resistance near the CNP. The variation of the resistance with the gate voltage and temperature is shown in Fig.4c where the evolution resembles that for single well 2D TIs[66]. The resistance decreases sharply at T > 15 K while saturating below 10 K, indicating a small mobility gap of 0.8 meV (see inset).

Figure4d,e shows Shubnikov-de Haas (SdH) oscillations measured in the region of electron conductivity as function of temperature and gate bias, respectively. The inset in Fig.4d displays that the temperature dependence of the SdH oscillations is well described by the Lifshitz-Kosevich formula. Surprisingly, Fig.4e presents a strong change in the phase of the SdH oscillations, as indicated by a dashed line that follows a constant phase. Sudden changes in the phase at 8 and 10 V could be related to variations in the Berry phase across system transitions.

Fourier analysis of the magnetoresistance in Fig.4e is displayed in Fig.5c. Two peaks corresponding to branches of conduction band carriers were obtained. The oscillations frequency increased with bias voltage and we observed the splitting of the upper frequency peak for sufficiently high gate voltage. The experimental results are compared with the theoretical model in the following section.

Figure 5

(a) Bulk and (b) edge energy dispersions for the nominal experimental parameters of the triple quantum well, i.e.,  nm and  nm, which falls close to the phase transition line, shaded rectangle in Fig. 2b. The lines and colors follow the same definitions as in Fig. 3. The system is in phase I, with the and subbands hybridized with nearly vanishing mass , yielding a linear pair of Dirac subbands (all the subbands are spin-degenerate). (c) The Fourier transform of the measured Shubnikov-de Haas oscillations show two main peaks splitting with increasing and a further Rashba-like splitting at high . (d) The theoretical SdH frequencies F for the linear and parabolic subbands from panel (a) qualitatively matches the experimental measurements. In both (c) and (d), the dashed lines mark for each subband with from panel (a) and the peaks are built from gaussian broadenings for easy comparison with panel (c).

Discussion

The triple well sample used in the experiment falls quite close to the phase transition line with  nm and  nm, as indicated by the shaded rectangle in Fig. 2b. These  nm uncertainties in and t are sufficient to locally shift the system between phases I and II along the sample. Nevertheless, here we focus on the nominal geometrical parameters  nm and  nm, and consider the fluctuations qualitatively in the following discussion.

As shown in Fig. 5, the second E-like subband and the H-like subband are nearly crossing, i.e. , forming a Dirac-like dispersion with a small gap of  meV, which is close to the experimental activation energy of 0.8 meV. The remaining two H subbands form parabolic dispersions with positive and negative effective masses. For a slightly larger or t (within the  nm experimental uncertainty) the and subbands would cross each other to yield the phase II regime (). Nevertheless, the inverted gap in phase II would still be small ( meV) and these subbands would still have a nearly linear Dirac dispersion. Moreover, with a small negative gap, the corresponding edge states would not be well defined, since its localization length is inversely proportional to this gap. Therefore, within the experimental uncertainty for and t, we would expect to have only one pair of well defined edge states near the phase transition line from phase I to II.

Local and non-local resistivities

To analyze the transport properties of the system in the experimental setup, first, let us consider the ballistic regime and ignore the bulk contributions to the local and non-local transport measurements. From nine terminals Landauer-Büttiker geometry, one would expect the local and non-local resistances to bewhere n is the number of edge states in each edge of the sample, and the labels i :  and v :  indicate the contacts used to apply the currents and measure the voltages, respectively, which were presented in the previous section. As discussed above, here we expect to have only pair of well defined edge states. However, as seen in Fig. 4b, both and deviate significantly from this ideal model for . Indeed, these deviations are justified by the bulk contributions for transport, since the edge states are immersed in the nearly gapless bulk, as seen in Fig. 5b. More specifically, this characteristic gives two opposite contributions. First, notice that the measurements in Fig. 4b decreases as the distance between contacts decreases, which is expected to asymptotically approaches the ballistic regime for short distances. Second, at the shortest distance, line C falls below the expected line for pairs of edge states. We interpret this as a consequence of finite bulk conductivity leading to reduced and . In summary, we expect scattering effects to lead to increased resistivities, while the bulk conductivity contributes to reduce the resistivities. These contrasting contributions show that transport measurements are ineffective to characterize if the system is in phase I or II.

Shubnikov-de Haas oscillations

Shubnikov-de Haas (SdH) oscillations are oscillations of the 2D magnetoresistivity of a material as a function of the external magnetic field B[67-70]. They were discovered in early 1930 and currently are one of the most important tools to access and extract values for both the 2D electronic and hole semiconductor densities at [70,71], the difference between electronic density of different subbands[71,72], the electron and hole effective masses[70,71,73] and also spin-orbit couplings e.g., Rashba[73-75]. Most recently, SdH oscillations have also been used to probe the 2D Dirac-like character of both graphene energy dispersion[76] and 2D Dirac-like surface states of 3D topological insulators[77,78]. For low temperatures (i.e., ), the SdH oscillations are well described by the Lifshits-Kosevich (LK)[69] expression, which in the absence of Zeeman and spin-orbit couplings reads[69,70,79-81]where is the differential longitudinal conductivity for spin and subband index j. For the experiments analyzed in this work, the most relevant subbands are the first two nearly spin-degenerate conduction subbands shown in Fig. 5a,b, i.e., one linear (L) and one parabolic (P), and hence we set , with corresponding energies , and . Throughout this work we always consider spin-degenerate subbands, despite the fact that at high electronic densities, a Rashba-like spin split is seen in the experimental data. We comment on this at the end of this section. In the generic expressions above, the temperature dependent term () and the Dingle factor () are written in terms of the density of states of each subband at the Fermi level , . The frequency , determining the SdH oscillations with respect to 1/B is written in terms of the charge density . The characteristics of the energy dispersion of each subband enters within and . Accordingly, for the linear subband we have and , while for the parabolic subband we have , and . The last parameter in above is the phase , which we discuss at the end of this section. Overall, the generic form of tell us that the experimentally measured frequencies can be used to obtain the different electronic densities of our subbands, independently of its linear or parabolic dispersion. This allows us to directly compare the experimental data given by Fig. 5c to a theoretical model in terms of a common total density axis, which is used to shift the lines for each SdH curve in Fig. 5c,d.

In practice, the experimental measurements of the SdH oscillations are done through the total differential resistivity . From the experimental data within Fig. 4a, we see that, for  V, we have . This translates to , which yields . A direct consequence of this approximation[82] is that, now, can be written as a weighted sum of the LK expressions above over our different subbands with corresponding different amplitudes , i.e.,

Here we write the weight factor as a reasonable, but rough estimate, that should be considered qualitative only. Once we have now obtained the theoretical equations that are going to be compared to the experimental quantities, we move to the experimental data.

In Fig. 5c we plot the Fourier transformation of the SdH oscillations of Fig. 4e for different values of . The peaks of Fig. 5c correspond to the different SdH frequencies [Eq. (14)]. Although the SdH frequencies are experimentally obtained as a function of different , which controls the Fermi energy, here we use the voltage-to-density conversion factor to perform a comparison between our theoretical results to the experimental data. Accordingly, we also plot the frequencies for different total density (left-hand side axis).

In order to compare the experimental results with our theory, in Fig. 5d we plot the SdH frequencies () predicted by our calculations as a function of (or the total 2D electronic density ). They are obtained using the parameters extracted from the conduction subband dispersions within Fig. 5a. Namely, we obtain for the parabolic H-like conduction subband, and  meV nm for the E-like linear conduction subband. At small frequencies F the theory in Fig. 5d predicts an extra peak arising from the second parabolic conduction subband [with subband edge at  meV in Fig. 5a, and ]. This peak is missing in the experimental data from Fig. 5c, which could be suppressed due to disorder, smaller Dingle factor, or by smaller weight in Eq. (15). Therefore, hereafter we focus on the two higher frequency branches in Fig. 5c,d. For small , we see that both frequencies and are close to each other. However, as is increased there is an evident separation between them, stemming from their different -dependencies, i.e., and . Surprisingly, our theoretical frequencies , plotted within Fig. 5d, present a good agreement with the experimental data shown in Fig. 5c. This shows that the SdH obtained in the experiments are consistent with the coexistence of a linear and parabolic subbands, similarly to the case of trilayer graphene[83]. This coexistence can also be evidenced by the different dependence of with respect to , which is shown by the black dashed lines in both Fig. 5c,d as a guide to our eyes. We emphasize that the only adjustable parameter between the experiment and theory presented in Fig. 5c,d is the voltage-to-density conversion factor (cm/V). Additionally, there are other features of the experimental data that corroborates the evidence of linear and parabolic subbands contributing to the SdH oscillations. We discuss them in the next paragraphs.

First, there is an apparent crossing of the two frequency peaks near  T within Fig. 5c. Since , such a crossing can only happen if there is an equivalent crossover in the densities of each subband. For a pair of parabolic subbands, this would require one of the subbands to be at a higher energy, and with a heavier mass (larger DOS). This is unlikely for our sample. In contrast, this crossing of frequencies is expected from the coexistence of linear and parabolic subbands nearly degenerated at . From the condition , it follows that the crossing point occurs for , where is the effective mass of the parabolic subband, and is the Fermi velocity of the linear subband.

Secondly, accordingly to Eqs. (15) and (16) the different amplitudes affect the relative spectral weight of the Fourier frequencies (). More specifically, for the parabolic subband, the DOS, , is constant as a function of , while for the linear subband we have . As a consequence, as we increase (or ), this leads to an increase of the differential resistivity of the linear subband, thus increasing its spectral weight. Surprisingly, this is also seen within the experimental data of Fig. 5c. While the spectral weight of remains nearly constant for different , the spectral weight of increases as a function of .

Thirdly, the temperature dependence of SdH oscillations amplitudes show an indirect evidence of the coexistence of the linear and parabolic subbands. From Eq. (11), we see that the temperature dependence is set by the terms of each subband. However, one cannot experimentally distinguish the individual contribution from each subband. Instead, the measurements shown in the inset of Fig. 4d are adjusted to be fitted by an effective LK expression with a single effective , defined by the replacement . Here, we can expect that this effective DOS should lie within a range set by the theoretical and . Indeed, for  V, we obtain from the experimental data  meVcm, while the theoretical expression, at the corresponding , gives us  meVcm and  meVcm. Moreover, for  V the experimental measurements result in a slightly increased  meVcm. Since for the parabolic subband is constant as a function of , the increase in comes from the linear subband, which is theoretically calculated and yields  meVcm for  V. In both and 12 V cases, the effective lies within the expected range , and it increases with due to the linear subband contribution.

To finish a complete characterization of the SdH oscillations, we would need to analyze the phases in the oscillations of Fig. 4, and also the beating patterns that typically arise from spin-orbit couplings (SOC). Indeed, it is well known[84-87] that for parabolic subbands the phase is , while for Dirac-like subbands and extra phase arises due to Berry’s curvature, yielding for the linear subband. However, current models for the SdH oscillations on linear Dirac subbands are limited. Typically, the models neglect the Zeeman splitting (for being small), and do not include the spin-orbit coupling[85-87]. Interestingly, it has been discussed that the Zeeman splitting introduces a phase correction to the SdH oscillations for Dirac subbands[77,82,88], but the effect of spin-orbit coupling is unknown at the moment and it cannot be inferred from knowledge of its counterpart in parabolic subbands. Since these developments are beyond the scope of this paper, we will leave a discussion of the phase and the splitting of the peaks in Fig. 5c at high densities for future works.

Conclusions

In summary, we have investigated the phase diagram of triple HgTe quantum wells as a function of its geometric parameters and compared its prediction with experimental measurements. For the theoretical investigation of the phase diagram we have projected the 3D Kane Hamiltonian into an effective 3BHZ 2D model that allowed us to investigate its edge state characteristics in each topological phase. We found that phases I and II are gapless due to small hybridization of H-like states from different quantum wells, but still present, respectively, one and two pairs of edge states immersed in the bulk. It is only in phase III that all E and H subbands are inverted and the three sets of edge states form within a bulk gap.

The experimental data, for a sample with geometric parameters that fall quite close to the transition from phase I to II, allowed us to analyze the predictions of the theoretical model and the consequences of having the edge states immersed in bulk. We have seen that non-local resisitivity measurements show a reduced signal due to bulk conductivity, while the local resistivity deviates from perfect quantization due to both bulk conductivity and non-ballistic transport. Consequently, models for edge states within bulk would have to account for these features to achieve reliable comparison with experiments. More interestingly, we have seen that SdH measurements show signatures of the predicted bulk subbands given by a set of linear and parabolic subbands near the phase transition. However, future work is needed to properly characterize the SdH phase for linear subbands in the presence of strong Zeeman and spin-orbit couplings.

Supplementary Information.