Engineering topological phases in triple HgTe/CdTe quantum wells

Quantum wells formed by layers of HgTe between Hg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{1-x}$$\end{document}1-xCd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_x$$\end{document}xTe 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8\times 8$$\end{document}8×8 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][2][3][4][5][6][7][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 semiconductorbased 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 (named after its authors Bernevig, Hughes, and Zhang).Additionally, the 2D Dirac-like band structure of 3D TIs have been observed in ARPES measurements [14][15][16][17][18][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][21][22][23][24] , strain 25,26 , and temperature 25,27 .The resulting 1D helical channels at the edges lead to quantized conductance and nonlocal edge transport, which has been observed for sufficiently short distances between measurement probes 28,29 .In all cases, the quantized conductance of HgTe-base 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 30,31 , charge puddles 32,33 , and many sources of inelastic scattering 34 .Despite the deviation from theoretical results, recent improvements in sample growth has lead to measurements closer to the predicted quantization 33,35 .
Recently, it has been proposed that the control of the layer localization of topological states in bilayer graphene 36,37 and TIs 38 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 arXiv:2112.01307v1[cond-mat.mes-hall] 2 Dec 2021 Refs. [39, 40], 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 41 that these DQWs could host second-order topological insulators with excitonic nodal phases that support flat band edge states, which could lead to superconductivity 42 .Experimentally, signatures of the conductance quantization, nonlocal transport and the Landau fan diagram of HgTe-based DQWs have been recently observed 43,44 .Additionally, magnetically doped TI layers, coupled through insulating spacers, lead to a solid state realization of the Weyl semimetal 45 , 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 46,47 , which can be realized in multiple coupled HgTe-based QWs 48,49 .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/Hg 1−x Cd x Te 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 8 × 8 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 states spread throughout the three QWs as symmetric and anti-symmetric combinations.For the E-like bands the hybridization leads to large energy splitting between these bands, while the H-like bands 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 bands 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 bands 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 conducting bands, which are shown to be in agreement with our band structure calculation.This analyzes 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.
2 } basis set.The CdTe bandstructure has a normal order, where Γ 6 is a S-type conduction band, Γ 8 and Γ 7 are the P-type valence bands corresponding to heavy and light holes (Γ 8 ) and the split-off band (Γ 7 ).In contrast, HgTe has the Γ 6 and Γ 8 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 hk k k → p p p as the momentum operator, and symmetrizes [50][51][52] the Hamiltonian.Additionally, the growth of the heterostructure typically induces strain, which is considered under the Bir-Pikus Hamiltonian 53 .The resulting 8 × 8 Kane Hamiltonian and the material parameters are shown in the Supplementary Material 54 .Here, the theoretical model set the growth direction z [001].Nevertheless, for small d 0 7.5 nm we expect the results to be nearly equivalent for growth directions [001] and [013] 55 .In the (x, y) plane, the solutions are given by plane-waves, ∝ e ik k k •r r r , where k k k = (k x , k y ) is the in-plane momentum.To numerically diagonalize the H Kane (z, p z , k k k ) we use the kwant code 56 , which provides an efficient interface to build and solve the numerical problem.

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 k k k = 0 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, 57] the authors project the total Hamiltonian into the k k k = 0 DQW eigenstates.Alternatively, in Ref. [39] 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. 1(a).For easy reading, we keep here a notation similar to Ref. [39], but we introduce the index ν = {L,C, R} 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 { H ν 1± , E ν 1± }, 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 39 into the triple well case, which we label "3×BHZ", and it reads as Here, the diagonal blocks of each layer, H ν , are given by a direct sum over Kramers partners with the Pauli matrices (σ x , σ y , and σ z ) acting on the (H 1 , E 1 ) subspace of each Kramers block.Similarly, the tunnel couplings V µ for µ = {LC,CR}, referring to the left-central and central-right QW couplings, are All coefficients above are calculated following the k k k • p p p perturbation theory up to second order.These are shown in the Supplementary Material 54 .

Theoretical Results
In this section 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 d 0 and t [see Fig. 1± and E R 1± .On the other hand, the three H-like bands (for each Kramers pair) are nearly degenerate because the heavy-hole states H ν 1± are strongly localized within each well due to their larger effective mass, and only show significant hybridization for t < 1 nm.As we increase d 0 in Figures 1(b)-(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 band flips the sign of its effective mass, but it still remains nearly degenerate with the other H-like bands at k k k = 0 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 Figs.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 bands obtained here also show the "camel back" profile 57,58 , thus also presenting an indirect gap defined by . This feature appears due to the strong hybridization of the QW subbands for large thickness d 0 7 nm (or small t < 1 nm), where the subbands are close (in energy) to each other [See Fig. 2(a)].Furthermore, for large d 0 ∼ 12 nm (not shown), the indirect gap ∆E I [see Fig. 1(b)] closes and the system becomes semi-metallic (SM).Notice that in Fig. 1(e) 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 ∆E I vanishes, while the "gapless" will refer to band structures with vanishing gap at k k k = 0.

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 59 , 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 k k k = 0. Using Q = {E, H} 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 Hµ with µ = {−1, 0, 1} have the usual BHZ form with Hµ = hµ with renormalized parameters Ãµ It is evident from Eq. ( 7) that the only hybridization happens between pairs 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 Mµ 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 d 0 and barrier thickness t.First, for fixed t = 3 nm in Fig. 2(a), 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, { E µ 1± , H µ 1± }, are highlighted (red circles) in Fig. 2(a) and also in Fig. 2(b), along the t = 3 nm line.In Fig. 2(b), the solid black lines mark the parameter values that yield crossings between subband pairs { E µ 1± , H µ 1± }.These lines delimit different regions of our phase diagram, which are labeled by the their corresponding spin Chern number, i.e., 0, I, II and III.This can be clearly seen in Figs.1(b)-(e), which contain the bandstructures corresponding to the cyan stars in Fig. 2(b).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 t 1 nm.Finally, the shaded region marks the semi-metallic phase, representing the cases where the indirect band gap ∆E I closes.The color map represents the nominal value of the gap at k k k = 0, which only becomes significant for t < 1 nm.

Edge states
To illustrate the characteristics of each phase presented above and in Fig. 2(b), 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 the confined system, we consider the effective 3×BHZ 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 |y| = L y /2, with L y ∼ 1000 nm.The spectrum is then calculated using a recursive Green's function method 60 , which allow us to calculate the local spectral function, A(E, k x ), for the bulk and edges states with an efficient exponential decimation.In Figs.lines), although they are not able to reproduce accurately the "camel back" profile around |k k k| 0.2 nm −1 .In Figs. 3 (a2), (b2), (c2) we have the spectrum corresponding to the cyan stars, I, II and III in Fig. 2(b), with gray bands representing bulk bands, and orange and green bands representing states localized at edges of our system.Interestingly, it is possible to identify here two different types of edge states.The ones indicated by the orange color are topological edge states with corresponding bands connecting conduction to valence band.Conversely, the ones indicated by the green color are trivial edge states predicted previously in Ref. [61].While the topological ones arise from the non-trivial topology of our system, the trivial ones appear when we confine bands that have a strong linear dispersion 61 .For this reason, these edge states appear due to the approximately chiral symmetry of the bands 61 .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 62 .

Experimental results
Triple quantum wells based on HgTe/Cd x Hg 1−x Te with [013] surface orientation and equal well widths of d 0 = 6.7 nm and barrier thickness t = 3 nm were prepared by molecular beam epitaxy (MBE).The sample structure is shown in Fig. 1(a).The layer thickness was determined by ellipsometry during MBE growth, with accuracy of ±0.3 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. 4(b)].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 three quantum wells, 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 5% 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 2 , was deposited on the sample surface and then covered by a TiAu gate.The density variation with the gate voltage is ∼ 0.9 × 10 11 cm −2 /V, estimated from the dielectric thickness and in 6/15 comparison with the calculated frequencies of the SdH oscillations shown in the next section.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.
Three devices from the same substrate were studied.Figure 4(a) 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 g ) 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 behaviour on the electronic side, however at the hole side the resistance shows monotonic behaviour due to strong scattering between the cone and the heavy hole branches 64 .In the CNP, the electron-like Hall resistance jumps from the negative quantized value ∼ h/e 2 to the hole-like positive value ∼ h/4e 2 .The quantum Hall effect in HgTe TQWs is beyond the scope of this work and will be reported in a forthcoming publication.
Fig. 4(b) 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. 4(c) where the evolution resembles that for single well 2D TIs 65 .The resistance decreases sharply at T > 15 K while saturating below 10 K, indicating a small mobility gap of 0.8 meV (see inset).
Fig. 4(d) and (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. 4(d) displays that the temperature dependence of the SdH oscillations is well described by the Lifshitz-Kosevich formula.Surprisingly, Fig. 4(e) 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. 4(e) is displayed in Fig. 5(c).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.

Discussion
The triple well sample used in experiment falls quite close to the phase transition line with d 0 = (6.7 ± 0.3) nm and t = (3.0 ± 0.3) nm, as indicated by the shaded rectangle in Fig. 2(b).These ±0.3 nm uncertainties in d 0 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 d 0 = 6.7 nm and t = 3 nm, and consider the fluctuations qualitatively in the following discussion.
As shown in Fig. 5, the second E-like subband E 0 1 and the H-like H 0 1 subband are nearly crossing, i.e.M0 ≈ 0, forming a Dirac-like dispersion with a small gap of M0 ∼ 0.6 meV, which is close to the experimental activation energy of 0.8 meV.The remaining two H bands form parabolic bands with positive and negative effective masses.For a slightly larger d 0 or t (within the ±0.3 nm experimental uncertainty) the E 0 1 and H 0 1 bands would cross each other to yield the phase II regime ( M0 < 0).Nevertheless, the inverted gap in phase II would still be small ( M0 ∼ −0.6 meV) and these bands 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 d 0 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 a nine terminals Landauer-Büttiker geometry, one would expect the local and non-local resistances to be where 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 n = 1 pairs of well defined edge states.However, as seen in Fig. 4(b), both R L and R NL deviate significantly from this ideal model for n = 1.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. 5(b).More specifically, this characteristic gives to opposite contributions.First, notice that the R L measurements in Fig. 4(b) 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 n = 1 pairs of edge states.We interpret this as a consequence of finite bulk conductivity leading to reduced R L and R NL .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 66-69 .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 B = 0 69,70 , the difference between electronic density of different subbands 70,71 , the electron and hole effective masses 69,70,72 and also spin-orbit couplings e.g., Rashba [72][73][74] .Most recently, SdH oscillations have also been used to probe the 2D Dirac-like character of both graphene energy dispersion 75 and 2D Dirac-like surface states of 3D topological insulators 76,77 .For low temperatures (i.e., k B T µ ≈ E F ), the SdH oscillations are well described by the Lifshits-Kosevich (LK) 68 expression, which in the absence of Zeeman and spin-orbit couplings reads 68,69,[78][79][80] ∆σ e − τ j,σ /τ 0 cos 2π where ∆σ (0) is the differential longitudinal conductivity for spin σ = {↑, ↓} and band index j.For the experiments analyzed in this work, the most relevant bands are the first two nearly spin-degenerate conduction bands shown in Fig. 5(a)-(b), i.e., one linear (L) and one parabolic (P) band, and hence we set j = {L, P}, with corresponding energies ε L = ±hv F k, and ε P = h2 k 2 /2m * .In this generic form, the temperature dependent term (λ j,σ ) and the Dingle factor (τ j,σ ) are written in terms of the density of states of each band at the Fermi level E F , g j,σ (E F ).The frequency F j,σ , determining the SdH oscillations with respect to 1/B is written in terms of the charge density n 2D, j,σ (E F ) ≈ E F 0 g j,σ (ε)dε.The characteristics of the energy dispersion of each band enters within g j,σ (E F ) and n 2D, j,σ (E F ). Accordingly, for the linear band we have , while for the parabolic band we have g P,σ (E F ) = m * /(2π h2 ), and n 2D,P,σ (E F ) = m * E F /(2π h2 ).However, the main difference between the linear (Dirac-like) and the parabolic bands lies in the phase φ j , which assumes the value φ P = 1/2 for parabolic bands and φ L = 0 for linear bands.For the case of linear bands, the absence of any phase within the argument of cos(...) in Eq. ( 11) is interpreted due to the non-zero (π) Berry's phase of 2D Dirac cones [81][82][83][84] .This introduces an extra π phase within the corresponding cos(...) argument, which effectively cancels out when combined with the previous existed π phase.Overall, the generic form of F j ∝ n 2D, j tell us that the experimentally measured frequencies can be used to obtain the different electronic densities of our bands, independently of its linear or parabolic dispersion.This allows us to directly compare the experimental data given by Fig. 5(c) to a theoretical model in terms of a common total density axis, which is used to shift the lines for each SdH curve in Figs.5(c) and (d).
In practice, the experimental measurements of the SdH oscilattions are done through the total differential resistivity ∆R/R 0 ∝ ρ xx .From the experimental data within Fig. 4(a), we see that, for V g −V CNP 4 V, we have ρ xx ρ xy .This translates ∑ j ,σ g j ,σ .
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. 5(c) we plot the Fourier transformation of the SdH oscillations of Fig. 4(e) for different values of V g −V CNP .The peaks of Fig. 5(c) correspond to the different SdH frequencies F j,σ [Eq.( 14)].Although the SdH frequencies are experimentally obtained as a function of different V g , which controls the Fermi energy, here we introduce a voltage-to-density conversion factor of ∼ 0.9 × 10 11 cm −2 /V.This is required in order to perform a comparison between our theoretical results to the experimental data.Moreover, this conversion factor shows to be in agreement with previous publications in similar devices 43,44 .Accordingly, we also plot the frequencies for different total density n 2D = ∑ j n 2D, j (left-hand side axis).
In order to compare the experimental results with our theory, in Fig. 5(d) we plot the SdH frequencies (F j,σ ) predicted by our calculations as a function of E F (or the total 2D electronic density n 2D ).They are obtained using the parameters extracted from the conduction bands dispersions within Fig. 5(a).Namely, we obtain m * = 0.0127m 0 for the parabolic H-like conduction band, and hv F = 400 meV.nm for the E-like linear conduction band.For small E F , we see that both frequencies F L,σ and F P,σ are close to each other.However, as E F is increased there is an evident separation between them, stemming from their different E F -dependencies, i.e., F L,σ ∝ n 2D,L,σ ∝ E 2 F and F P,σ ∝ n 2D,P,σ ∝ E F .Surprisingly, our theoretical frequencies F j,σ , plotted within Fig. 5(d), present a good agreement with the experimental data shown in Fig. 5(c).This shows that the SdH obtained in the experiments are consistent with the coexistence of a linear and parabolic band, similarly to the case of trilayer graphene 86 .This coexistence can also be evidenced by the different dependence of F j,σ with respect to E F , which is shown by the black dashed lines in both Figs. 5 (c) and (d) as a guide to our eyes.We emphasize that the only adjustable parameter between the experiment and theory presented in Figs.5(c) and (d) is the voltage-to-density conversion factor (∼ 0.9 × 10 11 cm −2 /V).Additionally, there are other features of the experimental data that corroborates the evidence of linear and parabolic bands contributing to the SdH oscillations.We discuss them in the next paragraphs.
First, there is an apparent crossing of the two frequency peaks near F = 6 T within Fig. 5(c).Since F j,σ ∝ n 2D, j,σ , such a crossing can only happen if there is an equivalent crossover in the densities n 2D, j,σ of each band.For a pair of parabolic bands, this would require one of the bands 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 bands nearly degenerated at k k k = 0. From the condition F P,σ = F L,σ , it follows that the crossing point occurs for E F = 2m * v 2 F , where m * is the effective mass of the parabolic band, and v F is the Fermi velocity of the linear band.
Secondly, accordingly to Eqs. ( 15) and ( 16) the different amplitudes P j,σ affect the relative spectral weight of the Fourier frequencies (F j,σ ).More specifically, for the parabolic band, the DOS, g P,σ (E F ), is constant as a function of E F , while for the linear band we have g L,σ (E F ) ∝ E F .As a consequence, as we increase E F (or n 2D ), this leads to an increase of the differential resistivity of the linear band, thus increasing its spectral weight.Surprisingly, this is also seen within the experimental data of Fig. 5(c).While the spectral weight of F P,σ remains nearly constant for different E F , the spectral weight of F L,σ increases as a function of E F .
Thirdly, the temperature dependence of SdH oscillations amplitudes show an indirect evidence of the coexistence of the linear and parabolic bands.From Eq. ( 11), we see that the temperature dependence is set by the λ j,σ terms of each band.However, one cannot experimentally distinguish the individual contribution from each band.Instead, the measurements shown in the inset of Fig. 4(d) are adjusted to be fitted by an effective LK expression with a single effective λ , defined by the replacement g j,σ (E F ) → ḡ(E F ). Here, we can expect that this effective DOS ḡ(E F ) should lie within a range set by the theoretical g P,σ and g L,σ .Indeed, for V g −V CNP = 6 V, we obtain from the experimental data ḡ = 0.07 × 10 11 meV −1 cm −2 , while the theoretical expression, at the corresponding E F , gives us g P,σ = 0.05 × 10 11 meV −1 cm −2 and g L,σ = 0.09 × 10 11 meV −1 cm −2 .Moreover, for V g −V CNP = 12 V the experimental measurements result in a slightly increased ḡ = 0.09 × 10 11 meV −1 cm −2 .Since g P,σ for the parabolic band is constant as a function of E F , the increase in ḡ comes from the linear band, which is theoretically calculated and yields g L,σ = 0.14 × 10 11 meV −1 cm −2 for V g −V CNP = 12 V.In both V g −V CNP = 6 and 12 V cases, the effective ḡ lies within the expected range g P,σ < ḡ < g L,σ , and it increases with E F due to the linear band contribution.
To finish a complete characterization of the SdH oscillations, we would need to analyze the phases φ j,σ in the oscillations of Fig. 4, and also the beating patterns that typically arise from spin-orbit couplings (SOC).However, current models for the SdH oscillations on linear Dirac bands are limited.Typically, the models neglect the Zeeman splitting (for being small), and do not include the spin-orbit coupling [82][83][84] .Interestingly, it has been discussed that the Zeeman splitting introduces a phase correction to the SdH oscillations for Dirac bands 76,85,87 , but the effect of spin-orbit coupling is unknown at the moment and it cannot be inferred from knowledge of its counterpart in parabolic bands.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. 5(c) 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 a sample that falls quite close to a phase transition line.For the theoretical investigation of the phase diagram we have projected the 3D Kane Hamiltonian into an effective 3×BHZ 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 bands 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 bands given by a set of linear and parabolic bands near the phase transition.However, future work is needed to properly characterize the SdH phase φ for linear bands in the presence of strong Zeeman and spin-orbit couplings.

8 × 8
Kane modelWe consider quantum wells (QWs) made of HgTe confined by Hg 1−x Cd x Te barriers with concentration x = 0.3, as shown in Fig.1(a).Both HgTe and CdTe crystallize in the zindblende structure with low energy bands around the Γ point (k k k = 0), which are well described by the 8 × 8 Kane Hamiltonians 50, 51 H Kane , with the corresponding {

Figure 1 .
Figure 1.(a) Conduction and valence band edges of the triple quantum well schematically shown in the bottom.The widths d 0 of the HgTe wells and the thickness t of the Hg 1−x Cd x Te barriers (x = 0.3) 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, while the H-like remain nearly degenerate due to the larger effective mass.(b,c,d) Band structures for t = 3 nm and varying d 0 .(b) For d 0 = 5 nm, in the trivial regime, all E bands are above the H bands. ∆E I indicates the indirect band gap.Phase transitions occur as each E band crosses the H bands with increasing (c) d 0 = 6 nm, (d) d 0 = 7 nm, and (e) d 0 = 8.5 nm.For large d 0 (not shown) the indirect gap ∆E I closes and the system becomes semi-metallic.
1(a)], and labeled by the corresponding spin Chern number.The edge state dispersions are illustrated for representative cases of each phase.Energy subbands Figures 1(b)-(e) illustrate the band structure for triple HgTe QWs with t = 3 nm and increasing d 0 = {5, 6, 7, 8.5} nm.In Fig. 1(b) the system is in the trivial regime i.e., all three conduction subbands have a predominantly Γ 6 electron-like (E-like) character, while the three valence subbands have a predominantly Γ 8 hole-like (H-like) character.Here, we have, for each Kramers pair, three non-degenerate conduction E-like subbands ( E − 1± , E 0 1± , and E + 1± ), which arise from the hybridization of the lowest conduction subbands of the left (L), central (C) and right (R) wells, namely, E L 1± , E C 1(c)-(d) are topologically non-trivial, but gapless.Moreover, it is only when all the E-like subbands are above [Fig.1(b)] or below [Fig.1(e)]

Figure 2 .
Figure 2. (a) Crossings of the E-like k k k = 0 band edges with the H-like bands as a function of d 0 , and t = 3 nm.(b) Phase diagram as a function of d 0 and t.The black solid lines mark the E-H band crossings, and the red circles along the t = 3 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 bands, ∆E H , which becomes clear only for t < 1 nm, and also splits the first solid black line for the E − 1 -H 1 crossings.The cyan stars along the t = 3 nm line mark the points (0, I, II, III) corresponding to the phases illustrated in Fig. 1(b)-(e).The shaded rectangle near d 0 = 6.7 nm and t = 3 nm marks the parameters of the experimental sample, with the area shaped to illustrate the experimental uncertainty of ±0.3 nm in both d 0 and t.

15 Figure 3 .
Figure 3. (a1, b1, c1) Comparison of the bulk band structures calculated with the 8 × 8 Kane model (black dashed lines) and the effective 2D 3×BHZ model (solid lines) with E-like bands in purple and H-like bands in blue.The panels correspond to the phases I, II and III indicated by the cyan stars in Fig. 2(b).(a2, b2, c2) Band structures for a nanoribbon geometry of width L y = 1000 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 3×BHZ bulk bands are shown as the solid lines as a guide to the eyes.(a2) For d 0 = 6 nm the toplogical regime with one E-band below the H bands shows one pair of topological edge state (orange) and two pairs of trivial edge states (green).(b2) As a second E-band crosses the H bands for d 0 = 7 nm one gets two topological edge states a single trivial edge state.(c2) For d = 8.5 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.

Figure 4 .
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 R ∼ exp(∆/2k B T ) with an activation energy 63 ∆ = 0.8 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 V g −V CNP 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.

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

9/15 to
85xy σ xx , which yields σ xx ≈ −ρ xx /ρ 2 xy .A direct consequence of this approximation85is that, now, ∆R/R 0 can be written as a weighted sum of the LK expressions above over our different bands with corresponding different amplitudes P j,σ , i.e.,