# Interaction between interface and massive states in multivalley topological heterostructures

G. Krizman, B. Assaf, M. Orlita, G. Bauer, G. Springholz, Robson Ferreira,

L. A. de Vaulchier, Y. Guldner

► To cite this version:

G. Krizman, B. Assaf, M. Orlita, G. Bauer, G. Springholz, et al., Interaction between interface and massive states in multivalley topological heterostructures. Physical Review Research, 2022, 4 (1), pp.013179. 10.1103/PhysRevResearch.4.013179. hal-03602971

HAL Id: hal-03602971

https://hal.science/hal-03602971v1

Submitted on 26 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Interaction between interface and massive states in multivalley topological heterostructures

G. Krizman,¹ B. A. Assaf,² M. Orlita,³,⁴ G. Bauer,¹ G. Springholz,¹ R. Ferreira,⁵ L. A. de Vaulchier,⁵ and Y. Guldner⁵

¹Institut für Halbleiter und Festkörperphysik, Johannes Kepler Universität, Altenbergerstraße 69, 4040 Linz, Austria

²Department of Physics, University of Notre Dame, Notre Dame, Indiana 46556, USA

³Laboratoire National des Champs Magnétiques Intenses (LNCMI), Centre national de la recherche scientifique (CNRS)-Université Grenoble Alpes (UGA)-Université Toulouse III-Paul Sabatier (UPS)-Institut National des Science Appliquées (INSA)-European Magnetic Field Laboratory (EMFL), 25 Avenue des Martyrs, F-38042 Grenoble, France

⁴Institute of Physics, Charles University, CZ-12116 Prague, Czech Republic

⁵Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, 24 rue Lhomond, 75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 10 November 2021; revised 17 January 2022; accepted 31 January 2022; published 4 March 2022)

Topological interface states (TISs) in multivalley systems are studied to unravel their valley sensitivity. For this purpose, multivalley IV-VI topological crystalline insulator (TCI) heterostructures are explored using magnetooptical Landau level spectroscopy up to 34 teslas. We characterize the TISs emerging from the distinct L valleys in Pb₁₋ₓSnₓSe multiquantum wells grown along the [111] direction. It is shown that the shape of the two-dimensional (2D) Fermi surfaces of TISs residing at the TCI-trivial insulator interfaces are strongly affected by the valley anisotropy of topologically trivial Pb₁₋ᵧEuᵧSe barriers. This phenomenon is shown to be due to the deep penetration of the TISs into the barriers. For the valleys tilted with respect to the confinement direction, a significant interaction between topological states and the conventional massive quantum well states is observed, evidenced by the resulting large anticrossings between Landau levels. These are theoretically well described by a **k**·**p** model that considers tilt and anisotropy of the valleys in 2D. Therefore, in this paper, we provide a precise characterization of the TIS valley splitting as well as an accurate determination of the anisotropy of their Dirac cone dispersion.

DOI: 10.1103/PhysRevResearch.4.013179

# I. INTRODUCTION

Multivalley semiconductors offer additional degrees of freedom for tuning of the electronic properties for valleytronics device applications because these properties significantly depend on the structure of the valleys as well as their interaction [1–5]. It is, however, a major challenge to control and determine the valley splitting and valley anisotropy in semiconductors. For instance, qubit manipulation in multivalley systems like Si, Si/Si₁₋ₓGeₓ, or Ge/Ge₁₋ₓSiₓ heterostructures, which is based on spin manipulation, requires a fine-tuned valley splitting to ensure Pauli blockade and avoid decoherence [6–8]. Optimized valley splitting has recently allowed us to obtain qubits at an enhanced operating temperature [9,10]. The additional valley degree of freedom can also be used to control the electric current and constitutes the basis of valleytronics device applications such as valley-based classical bits [11,12] and/or qubits [13,14]. From a fundamental point of view, the valley anisotropy is also at the origin of spontaneous crystal symmetry breakings, leading to nematic phases, or the quantum Hall ferroelectric state [15–18]. The lead

salt semiconductors are a particularly interesting template for valleytronics not only because of their multivalley band structure [19,20] but also because, through alloying with tin, they can be converted into topological crystalline insulators (TCIs) featuring massless and spin-polarized Dirac electrons with the spin locked to the momentum [21–27]. Therefore, they can be envisioned for topological valleytronics in which information is stored, e.g., by the different chirality of electrons in different valleys.

In lead salt compounds, the band extrema are located at the four L points of their Brillouin zone (BZ), as illustrated in Fig. 1(a) [19,28], where in each valley, the band structure consists of two mirror bands $L_6^+$ and $L_6^-$ with opposite symmetry [see Fig. 1(b)] [19,29,30]. In pseudobinary Pb₁₋ₓSnₓSe and Pb₁₋ₓSnₓTe alloys, this symmetry can be inverted when the Sn content exceeds a certain critical value [20,22–24,26,27,31,32]. This leads to the emergence of the TCI phase, whose hallmark is the existence of Dirac (interface or surface) states protected by crystalline symmetry at the (001) and (111) crystal faces [33,34]. Due to the multivalley band structure, for the (111)-oriented surface, three Dirac cones appear at the $\bar{M}$ points of the two-dimensional (2D) BZ, and one is located at the center $\bar{\Gamma}$, as depicted in Fig. 1(a). The $\bar{M}$ and $\bar{\Gamma}$ points can be regarded as the projection of the four L points in the three-dimensional (3D) BZ to the 2D BZ of the (111) surface, corresponding to the oblique and longitudinal valleys, respectively. Extensive studies have been performed to quantify the valley anisotropy in bulk lead salt materials

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Three-dimensional (3D) Brillouin zone of cubic lead salt crystals and its two-dimensional (2D) projection onto the (111) surface. Band minima in 3D are located at the four different *L* points: three in red and one in black. The surfaces of constant energy of these band minima are plotted for the two types of valleys with respect to [111]: the longitudinal valley depicted in black and the three equivalent oblique valleys in red, tilted by $\theta = 70.5^\circ$ with respect to the surface normal [111]. They both have an anisotropic factor *K*, defined as the ratio between the major and minor axis of the ellipsoids. The projection onto the (111) surface corresponds to the location of the topological interface state (TIS) Dirac cones. Three at the $\bar{M}$ points in red stem from the band inversion in the oblique valleys and one at the $\bar{\Gamma}$ point in black originates from the longitudinal valley. Note that these cones can be gapped with a hole part TIS' and an electron part TIS. (b) Schematic band alignment of $L_6^\pm$ bands between Pb$_{0.75}$Sn$_{0.25}$Se (well) and Pb$_{1-y}$Eu$_y$Se (barrier) along Z parallel to the [111] growth direction. The confined state energy levels at $k_\parallel = 0$ are represented by dashed lines in blue for TIS and TIS' and in orange for massive states, labeled H2, H3, E2, E3, ... (c) Square of the $F_+$ component of the envelope function spinor for TIS in blue and the two first massive confined states E2 and E3 in orange calculated using the parameters for the sample multiquantum well (MQW)-39, see Table I. The curves are vertically shifted for clarity.</figcaption></figure>

[32,35–38] and at free surfaces of Pb$_{1-x}$Sn$_x$Se TCIs using surface-sensitive techniques [24,27,39]. Here, we focus on the topological states at buried interfaces of topological multilayer structures probed by magnetooptical spectroscopy to accurately characterize both interface and bulk states, thereby revealing the significant impact of the bulk on the valley splitting and anisotropy in heterostructure systems.

The electronic properties of lead salt compounds are well-described by a 4-band **k**·**p** Hamiltonian, which is like a massive anisotropic 3D Dirac Hamiltonian [32,40]. The constant energy surfaces are ellipsoids with an anisotropy factor *K* defined by the ratio between the main and minor axes of the ellipsoids [see Fig. 1(a)] [38,41]. Although *K* is identical for the four valleys, their different locations within the BZ
lead to different orientations of their major axes. One ellipsoid has its long main axis aligned along the [111] direction, and the three others are tilted by an angle $\theta = 70.5^{\circ}$ with respect to [111]. We will refer to these as longitudinal and oblique valleys, respectively, as drawn in Fig. 1(a).

Here, we focus on the multivalley features of multiquantum well (MQW) heterostructures composed of topologically nontrivial Pb$_{0.75}$Sn$_{0.25}$Se quantum wells (QWs) separated by topologically trivial barriers of Pb$_{1-y}$Eu$_y$Se, as shown schematically in Fig. 1(b) [22,24,32,42]. In addition to the usual quantum confined states localized within the QW [orange states in Figs. 1(b) and 1(c)], the switching of the topological character at each interface gives rise to topological interface states (TISs) localized at the heterointerfaces, shown in blue in Figs. 1(b) and 1(c) and already observed in our previous work [43]. The relatively small thicknesses of the wells yield hybridization between top and bottom TISs of each QW that opens a gap $2\delta$ in the Dirac cone, between the electronlike and the holelike TISs, here denoted as TIS and TIS', respectively [43–45]. Surprisingly, using magnetooptical Landau level spectroscopy, we show that TIS and TIS' are very sensitive to the valley anisotropy of the bulk barrier material due to their remarkable property of being localized mainly at the interface, contrary to trivial QW states found in ordinary semiconductor heterostructures. Even though the QW material Pb$_{0.75}$Sn$_{0.25}$Se has rather isotropic valleys, the valley anisotropy of the barrier material Pb$_{1-y}$Eu$_y$Se is sufficient to strongly impact the dispersion of the TIS. As a result, the interaction between the TIS Landau levels and those of trivial QW states is allowed and generates an observable signature in Landau level spectroscopy. Lastly, we show that this impact is unique to the $\bar{M}$ points of the 2D BZ and does not occur at the $\bar{\Gamma}$ point.

The paper is divided as follows. In Sec. II, the experimental results are introduced. We develop in Sec. III a complete theory for the confinement effect on a 3D anisotropic Dirac model, considering multivalley band structure, which applies for lead salts. We show that the anisotropy and valley tilt of the parent bulk materials allow for a finite interaction between TISs and massive confined QW states [indicated in blue and orange states in Figs. 1(b) and 1(c)], lifting valley degeneracy at some magnetic fields or in-plane momenta. The analysis of the experimental results using the theory is presented in Sec. IV. We experimentally demonstrate that this anisotropy induces interactions between these states, observable as a set of repeated anticrossings that occur between the Landau levels of topological and trivial states. Our spectroscopy also allows us to quantify the strength of the interaction potential between different levels and to accurately deduce the in-plane anisotropy of the TIS Dirac cone dispersion [20,46].

## II. PRESENTATION OF THE EXPERIMENTS

The investigated samples were grown by molecular beam epitaxy on freshly cleaved (111) BaF$_2$ substrates following the procedures detailed in our previous works [43,47]. Two MQWs were prepared with ~15 periods of Pb$_{0.75}$Sn$_{0.25}$Se (well material) and Pb$_{1-y}$Eu$_y$Se (barrier material) to which we refer as MQW-39 and MQW-25 according to their well thicknesses (see Table I). Note that the barrier material composition

<table><caption>TABLE I. Structure parameters of the Pb$_{0.75}$Sn$_{0.25}$Se bulk reference sample and the two Pb$_{0.75}$Sn$_{0.25}$Se TCI MQWs with Pb$_{1-y}$Eu$_y$Se barriers used in this paper. NA = not applicable.</caption><thead><tr><th>Parameter</th><th>Bulk</th><th>MQW-39</th><th>MQW-25</th></tr></thead><tbody><tr><td>d<sub>PbSnSe</sub> (nm)</td><td>2000</td><td>39 ± 2</td><td>25 ± 1</td></tr><tr><td>x<sub>Sn</sub></td><td>0.25 ± 0.01</td><td>0.25 ± 0.01</td><td>0.25 ± 0.01</td></tr><tr><td>d<sub>PbEuSe</sub> (nm)</td><td>NA</td><td>200</td><td>120</td></tr><tr><td>y<sub>Eu</sub></td><td>NA</td><td>0.10 ± 0.01</td><td>0.05 ± 0.01</td></tr><tr><td>Period number</td><td>NA</td><td>15</td><td>16</td></tr></tbody></table>

is $y = 0.05$ and $0.1$ for MQW-25 and MQW-39, respectively. As a reference sample, a 2-µm-thick film of Pb$_{0.75}$Sn$_{0.25}$Se was grown. The Sn content $x = 0.25 \pm 0.01$ in the three samples ensures a topological phase for temperature $\sim 100$ K [21,42,43,48,49]. The ratio between the barrier and the well thicknesses is fixed as 5 : 1 to suppress any coupling of confined states between adjacent QWs. To obtain a low carrier concentration in the three samples, Bi doping was employed during growth to compensate the native p-type character of Pb$_{0.75}$Sn$_{0.25}$Se due to Se vacancies. Carrier densities as low as a few $10^{17}$ cm$^{-3}$ and mobilities $> 10^4$ cm$^2/(Vs)$ are obtained in these samples [43,50].

Infrared magnetospectroscopy is used to characterize the electronic properties of Pb$_{1-x}$Sn$_x$Se MQWs and the TIS band structure. Experiments in magnetic fields up to 34 T are performed at the Laboratoire National des Champs Magnétiques Intenses in Grenoble. A setup as previously described in Ref. [51] is used to measure the infrared transmission of the samples with a magnetic field B// [111] (Faraday geometry). The transmitted radiation through the samples is collected in the same bath by a composite Si bolometer operating at 1.8 K, which is thus also the sample temperature.

Figure 2 compares the magnetooptical fan charts obtained on the bulk sample [Fig. 2(a)] with that of MQW-39 [Fig. 2(b)]. The dots are obtained by pinpointing the energy of the experimentally determined transmission minima and plotting these vs magnetic field, and the solid lines in Fig. 2(a) represent the Landau level transitions calculated by solving the 4-band $k \cdot p$ Hamiltonian for bulk material, as described below. In the Faraday geometry, the interband transitions obey the selection rules $n \rightarrow n \pm 1$ with a Landau level spin flip [52]. The data are in perfect agreement to our model using the parameters listed in Table II. For the given composition, we obtain a negative bulk bandgap $2\Delta = -52.5$ meV, in good agreement with previous studies [25,32,48,53]. Only one series of Landau levels with $v_\parallel = 4.50 \times 10^5$ m/s is required to explain our results, meaning that the longitudinal and oblique valleys are indistinguishable in the bulk material, i.e., the constant energy surfaces are approximately spheres ($K \sim 1$).

Comparing Figs. 2(a) and 2(b), a striking difference between the fan charts of the bulk and the MQW-39 samples is observed. Whereas the bulk sample displays a high regularity of the transmission minima vs magnetic field that evolves as $\sim \sqrt{B}$ as expected for bulk Dirac materials, in the MQW, a pronounced nonmonotonic dependence of the interband transition energies as a function of B emerges. In the following, we will show that this behavior arises from an interaction
<figure><img src="image_4.png" /><figcaption>FIG. 2. Magnetooptical transition fan charts for (a) the bulklike Pb<sub>0.75</sub>Sn<sub>0.25</sub>Se film and (b) multiquantum well (MQW)-39. Dots represent the energy of transmission minima plotted vs applied magnetic field. For MQW-39, two series of transitions are highlighted to emphasize the difference between (a) and (b). In (a), the solid lines represent the calculated interband Landau level transitions for the thick film sample that perfectly fit to the experimental data, whereas for the MQW, the dependence is highly nonmonotonic. The interband transitions are labeled by the indices *n* of the two involved Landau levels. The shaded area denotes the experimentally nonreachable range of photon energies <60 meV due to energy cutoff of the ZnSe windows used in the experimental setup.</figcaption></figure>

between TIS and massive confined states in the MQW. This coupling is an inherent consequence of the topological and valley degenerate character of Pb<sub>1-x</sub>Sn<sub>x</sub>Se.

### III. THEORY

The electronic properties of lead salts including Pb<sub>1-x</sub>Sn<sub>x</sub>Se or Pb<sub>1-y</sub>Eu<sub>y</sub>Se can be well described by a 4-band **k**·**p** Hamiltonian, written near the *L* points on the basis of the *L*<sub>6</sub><sup>+</sup> and *L*<sub>6</sub><sup>-</sup> conduction and valence bands

<table><caption>TABLE II. Band parameters determined by the fit of the magnetooptical data in the three investigated samples. NA = not applicable.</caption><thead><tr><th>Parameter</th><th>Bulk</th><th>MQW-39</th><th>MQW-25</th></tr></thead><tbody><tr><td>2Δ (meV)</td><td>-52.5</td><td>-45</td><td>-55</td></tr><tr><td>V (meV)</td><td>NA</td><td>250</td><td>175</td></tr><tr><td>m&#773; (m<sub>0</sub>)</td><td>0.35</td><td>0.35</td><td>0.35</td></tr><tr><td>v<sub>||</sub><sup>l</sup> (10<sup>5</sup> m/s)</td><td>4.50</td><td>4.60</td><td>4.70</td></tr><tr><td>v<sub>z</sub><sup>l</sup> (10<sup>5</sup> m/s)</td><td>4.50</td><td>4.40</td><td>4.40</td></tr><tr><td>K<sub>F</sub></td><td>NA</td><td>1.00</td><td>1.00</td></tr><tr><td>v<sub>||</sub><sup>o</sup> (10<sup>5</sup> m/s)</td><td>4.50</td><td>4.90</td><td>4.95</td></tr><tr><td>v<sub>z</sub><sup>o</sup> (10<sup>5</sup> m/s)</td><td>4.50</td><td>3.90</td><td>4.15</td></tr><tr><td>K<sub>M</sub></td><td>NA</td><td>1.23</td><td>1.17</td></tr></tbody></table>

[19,32,41]:

$$H = \begin{bmatrix} -\left(\Delta + \frac{p_{\parallel}^2}{2\tilde{m}} + \frac{p_z^2}{2\tilde{\mu}}\right)\mathbf{1} & v_{\parallel}\mathbf{p}_{\parallel} \cdot \boldsymbol{\sigma}_{\parallel} + v_z p_z \boldsymbol{\sigma}_z \\ v_{\parallel}\mathbf{p}_{\parallel} \cdot \boldsymbol{\sigma}_{\parallel} + v_z p_z \boldsymbol{\sigma}_z & \left(\Delta + \frac{p_{\parallel}^2}{2\tilde{m}} + \frac{p_z^2}{2\tilde{\mu}}\right)\mathbf{1} \end{bmatrix}. \quad (1)$$

In this equation, **1** is the 2×2 identity matrix and **σ** the Pauli matrices. Here, 2Δ denotes the bulk bandgap of the system, which is positive for Pb<sub>1-y</sub>Eu<sub>y</sub>Se and negative for Pb<sub>0.75</sub>Sn<sub>0.25</sub>Se at *T* = 1.8 K. Here, m&#773; represents the contribution of all the remote bands to the in-plane effective mass. It exerts a rather small influence, as these are at ~1 eV above and below *L*<sub>6</sub><sup>+</sup> and *L*<sub>6</sub><sup>-</sup> [30]. Note that the far-band correction μ&#773; to the longitudinal effective mass can be safely neglected, as it has been demonstrated to be very small in lead salt compounds [52]. In Eq. (1), the anisotropy is encoded by the in-plane and longitudinal Dirac velocities *v*<sub>||</sub> and *v*<sub>z</sub> respectively, and thus, the anisotropy factor of bulk states is given by *K* = *v*<sub>||</sub>/*v*<sub>z</sub>. The Hamiltonian is written for the *z* axis parallel to the major axis of the ellipsoid.

### A. Bulk Landau levels

To interpret the magnetooptical data with the magnetic field parallel to [111], we calculate in the following the
Landau levels. For the longitudinal valley, the z axis is aligned with the [111] direction and therefore with **B**. Here, $p_x$ and $p_y$

are expressed in terms of the ladder operators following the Peierls substitution. The Hamiltonian at $k_z = 0$ is given by

$$H = \begin{bmatrix} -\Delta - \hbar\tilde{\omega}(n-1) & 0 \\ 0 & -\Delta - \hbar\tilde{\omega}(n+1) \\ 0 & v_{\parallel}\sqrt{2e\hbar Bn} \\ v_{\parallel}\sqrt{2e\hbar Bn} & 0 \end{bmatrix} \begin{array}{c} \frac{v_{\parallel}\sqrt{2e\hbar Bn}}{\Delta + \hbar\tilde{\omega}(n-1)} \\ 0 \\ 0 \\ 0 \end{array}, \quad (2)$$

where $n$ is the Landau level index, and $\tilde{\omega} = eB/\tilde{m}$. As for the effective mass corrections, the remote bands induce a small correction to the $g$ factor that accounts for $\hbar\tilde{\omega}/2$ in this system [50,54]. The Hamiltonian in Eq. (2) leads to Dirac-like Landau levels, whose energies are given by

$$E_n = \pm \hbar \tilde{\omega} \pm \sqrt{(\Delta + \hbar \tilde{\omega} n)^2 + 2e\hbar \left(v_{\parallel}^2 + \frac{\Delta}{\tilde{m}}\right) Bn}. \quad (3)$$

These energies were computed at $k_z = 0$, where the joint density of state of the interband magnetooptical transitions is maximal.

For the oblique valleys, the main ellipsoidal axes are not aligned with the magnetic field. It is therefore convenient to introduce a new coordinate system (X, Y, Z) with the Z axis parallel to **B**, thus parallel to [111]. Here, (x, y, z) and (X, Y, Z) are connected by a rotation combined with a scale change according to relative sizes of $v_\parallel$ and $v_z$. Using the Hamiltonian transformation as described in Refs. [40,41], the Landau levels of the oblique valleys are given by Eq. (3) in which $v_\parallel$ is replaced by [54,55]

$$v_{\parallel} \left( \cos^2 \theta + \frac{1}{K^2} \sin^2 \theta \right)^{1/4}.$$

The Landau levels are thus described by Eq. (3) for both types of valleys using two different electron velocities, which account for the anisotropy effect in the bulk sample.

## B. QW Landau levels

Confinement effects in the QWs are calculated by adding a square potential $V(Z)$ into Eq. (1), where $Z$// [111] is the heterostructure growth axis. The confinement potential $V(Z)$, whose value depends on the Eu content $y$, is electron-hole symmetric for PbSnSe/PbEuSe MQWs [35,43], as shown schematically in Fig. 1(b), where the variation of the band edges is drawn vs the growth direction. The potential $V(Z)$ makes [111] a preferential direction that will clearly differentiate the properties of the two types of valleys: longitudinal and oblique.

Let us first compute the energy and wave function of the $j$th confined state at $k_\parallel = 0$. To this end, we follow the procedure carried out in Ref. [43], which consists of solving the following equation and applying the probability current continuity conditions at each interface [56–58]:

$$\begin{bmatrix} -\Delta + V(Z) - E^{(j)} & i\hbar\xi v_z \frac{d}{dZ} \\ -i\hbar\xi v_z \frac{d}{dZ} & \Delta + V(Z) - E^{(j)} \end{bmatrix} \begin{bmatrix} F_+^{(j)} \\ \xi F_-^{(j)} \end{bmatrix} = 0,$$

where $\xi = \pm 1$ represents the two spin components, and $F_\pm^{(j)}$ are the components on $L_6^\pm$ of the envelope function. Solving this equation leads to two confined states TIS and TIS' [plotted in blue in Figs. 1(b) and 1(c)], directly emerging from the topological character of the well material, as well as the usual confined states (plotted in orange) which emerge from the quantification of $L_6^\pm$. While TIS and TIS' wave functions are peaked at the interfaces, as shown in Fig. 1(c), the higher energy confined states (E2, H2, ...) remain mainly localized in the heart of the QW, thus, representing massive QW states. If we consider the parity of the $F_+$ component of the $j$th confined state envelope wave function, TIS and E3 are odd, while E2 is even, and so on for the following confined states. As $F_-$ is proportional to $dF_+/dZ$, the parity of $F_-$ follows the opposite interplay.

For the calculation of the in-plane motions of the longitudinal valley, no Hamiltonian transformation is needed. We treat the $k_\parallel$ terms (or $B$ terms if we consider an applied magnetic field) in perturbation. It requires us to calculate the matrix elements $\langle F_\pm^{(i)} | \delta W | F_\pm^{(j)} \rangle$ where $1 \le i, j \le N$ with $N$ the number of confined states obtained by solving the previous equation for $k_\parallel = 0$, and $\delta W$ is given by

$$\delta W = \begin{pmatrix} -\frac{p_{\parallel}^2}{2\tilde{m}} \mathbf{1} & v_{\parallel} \mathbf{p}_{\parallel} \cdot \mathbf{\Sigma}_{\parallel} \\ v_{\parallel} \mathbf{p}_{\parallel} \cdot \mathbf{\Sigma}_{\parallel} & \frac{p_{\parallel}^2}{2\tilde{m}} \mathbf{1} \end{pmatrix}.$$

For an applied magnetic field **B**// [111], the numerical solutions yield the Landau levels of the longitudinal valley. Those originating from the conduction band states are shown in Fig. 3(a) for MQW-39. Symmetric (i.e., related to energies of opposite sign) Landau levels are obtained for the confined hole states, as $V(Z)$ is electron-hole symmetric. Note that no interactions are obtained between Landau levels.

For the oblique valleys, an axis rotation is needed to account for their tilt with respect to the [111] orientation of the samples [55], as described in detail in Appendix A. As a result, the following Hamiltonian is obtained:

$$\begin{bmatrix} -[\Delta - V(Z)]\mathbf{1} & v_z P_Z \Sigma_Z \\ v_z P_Z \Sigma_Z & [\Delta + V(Z)]\mathbf{1} \end{bmatrix} + \begin{pmatrix} -\frac{P_{\parallel}^2}{2\tilde{m}} \mathbf{1} & v_{\parallel} \mathbf{P}_{\parallel} \cdot \mathbf{\Sigma}_{\parallel} \\ v_{\parallel} \mathbf{P}_{\parallel} \cdot \mathbf{\Sigma}_{\parallel} & \frac{P_{\parallel}^2}{2\tilde{m}} \mathbf{1} \end{pmatrix} + \sin\theta \begin{bmatrix} -\frac{1}{2\tilde{m}}h(X, Z) & (v_{\parallel} - v_z)\mathcal{H}(X, Z) \\ (v_{\parallel} - v_z)\mathcal{H}(X, Z) & \frac{1}{2\tilde{m}}h(X, Z) \end{bmatrix}. \quad (4)$$

Here, **P** and **Σ** are the momentum and Pauli operators in the new coordinate system (X, Y, Z). The two first terms of Eq. (4) are identical to those of the longitudinal valley. The third term
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Calculated Landau level energies of confined states at the $\bar{\Gamma}$ point for multiquantum well (MQW)-39 for parameters given in Table II. The Landau level from topological interface states (TISs; massive states) are drawn in red (black). Some of the Landau level indices are indicated. (b) Landau level energies of confined states at the $\bar{M}$ point for MQW-39. (c) Expanded view of the region encompassed by the yellow contour shown in (b). Dark and bright blue lines refer the different spin states. The anticrossings due to the two coupling potentials $W_1$ and $W_2$ are highlighted.</figcaption></figure>

includes the anisotropy effects with ($v_\parallel - v_z$) as well as the tilt effect with $\sin\theta$. It is proportional to $\sin\theta$ and thus cancels for the longitudinal valley. Its diagonal elements write as

$$ \pm \frac{1}{2\tilde{m}} h(X, Z) = \pm \frac{1}{2\tilde{m}} [\sin\theta(P_Z^2 - P_X^2) - 2\cos\theta P_Z P_X] \mathbf{1}. \quad (5) $$

These diagonal terms come from the anisotropic far-band correction terms that are $\tilde{m}$ in the $XY$ plane and $\tilde{\mu} \sim \infty$ along $Z$. Similarly, the off-diagonal elements are proportional to ($v_\parallel - v_z$) = $v_z(K-1)$ and account for the anisotropy effect of the valence and conduction bands. They are given by

$$ (v_\parallel - v_z)\mathcal{H}(X, Z) = (v_\parallel - v_z)[(\sin\theta P_Z - \cos\theta P_X)\Sigma_Z - (\sin\theta P_X + \cos\theta P_Z)\Sigma_X]. \quad (6) $$

Solving Eq. (4) in the framework of the perturbation theory, which is detailed in Appendix B, yields the subband dispersions of the oblique valleys as well as the corresponding Landau levels. The intricate behavior of these Landau levels in which many anticrossings appear is revealed in Fig. 3(b). We focus on the first two anticrossings zoomed-in in Fig. 3(c). They occur between Landau levels originating from TIS and E2 or E3 confined states. More generally, all these anticrossings stem from two distinct kinds of couplings $W_1$ and $W_2$. The first coupling $W_1$ between the $i$th and $j$th confined states is given by

$$ |W_1(nB)| = \sqrt{\frac{\hbar e B n}{2} \frac{\sin 2\theta}{2}} \left| \int_{-\infty}^{+\infty} \left\{ \frac{\hbar}{\tilde{m}} \left[ F_+^{(i)} \frac{dF_+^{(j)}}{dZ} + F_-^{(i)} \frac{dF_-^{(j)}}{dZ} \right] - (v_\parallel - v_z) [F_+^{(i)} F_-^{(j)} - F_+^{(j)} F_-^{(i)}] \right\} dZ \right|. \quad (7) $$

The parities of $F_\pm^{(i)}$ and $F_\pm^{(j)}$ impose $j = i \pm 1$ so that this coupling involves two successive confined states. The anticrossings take place between two Landau levels of successive subbands with indices $n$ and $n \pm 1$ and identical spins (see Appendix B). As an example, we highlight the anticrossings that occur for magnetic field $\sim 6$ T and an energy $\sim 40$ meV between subbands TIS and E2 on Figs. 3(b) and 3(c).

The second potential $W_2$ is responsible for an interaction between TIS and E3 due to the parities of $F_\pm^{(i)}$ and $F_\pm^{(j)}$. It involves Landau levels with indexes $n$ and $n \pm 1$ and opposite spins. Such anticrossings correspond, for instance, to the one computed $\sim 14$ T and 60 meV in Figs. 3(b) and 3(c). This second coupling $W_2$ is given by

$$ |W_2| = (v_\parallel - v_z)\hbar \frac{\sin 2\theta}{2} \left| \int_{-\infty}^{+\infty} \left[ F_+^{(i)} \frac{dF_-^{(j)}}{dZ} + F_-^{(i)} \frac{dF_+^{(j)}}{dZ} \right] dZ \right|. \quad (8) $$

We highlight that $W_1$ and $W_2$ are proportional to the tilt $\theta$ of the valley axis with respect to the growth direction and to the valley anisotropy ($v_\parallel - v_z$). Note also that the confinement, whether it occurs in the longitudinal or in the oblique valleys, leads to different energy quantization of the subbands, as shown by the extrapolations to $B=0$ of the distinct Landau level series of Figs. 3(a) and 3(b). This valley splitting is due to the valley-dependent effective mass along the confinement direction [35,59]. For instance, the valley splitting for the confined state E4 is calculated as 3 meV.

As a remark, it is found that the slope of the $n=0$ Landau levels depend on the valleys and are different for TIS and massive states. This slope scales as $\hbar e/\tilde{m}$ at the $\bar{\Gamma}$ point and as $\hbar e \cos^2\theta/\tilde{m}$ for oblique valleys. Therefore, as seen in Fig. 3(b), the slope of the $n=0$ Landau levels are drastically reduced at
<figure><img src="image_6.png" /><figcaption>FIG. 4. (a) Magnetooptical fan chart of multiquantum well (MQW)-39. Dots are experimental transmission minima, and lines are calculated transitions between Landau levels. The red (black) color denotes absorptions that mainly involve topological interface states (massive states). The solid lines correspond to the oblique transitions and the dashed lines to the longitudinal valley transitions. The transition labeling is only efficient for the longitudinal valley (see Fig. 3). (b) Magnetooptical spectra within the yellow rectangle in (a), highlighting the characteristic anticrossings. The two absorptions are marked by arrows whose transparencies mimic absorption intensities.</figcaption></figure>

the $\bar{M}$ points because of the anisotropic far-band contribution. Furthermore, this slope is negative for massive states and positive for TIS. This effect accounts for the band symmetry of the 2D confined states. The massive states stem from the QW material bands, which are inverted, while TIS and TIS' emerge in a trivial ordering in IV–VI. As the symmetry ordering of the remote bands remains identical, their effects on $L_6^\pm$ are constant and tend to push up (down) in magnetic field a $n=0$ Landau level that has a $L_6^-$ ($L_6^+$) symmetry. Thus, the sign of the $n=0$ Landau level slope is a direct indication of the confined state Bloch function symmetries. The far-band contribution, which is sometimes called inversion parameter, is therefore a good marker of band inversion in this system.

## IV. RESULTS AND DISCUSSION

The model derived for bulk $Pb_{1-x}Sn_xSe$ samples (see Sec. III A) is used to interpret the magnetooptical fan chart presented in Fig. 2(a). Only one Dirac velocity is used to fit the data with Eq. (3), meaning that no anisotropy is measured. The absence of valley splitting in the magnetooptical data of the bulk sample implies that indeed $Pb_{0.75}Sn_{0.25}Se$ is a quasi-isotropic system, as we pointed out previously [32]. In the following, we show, however, that in $Pb_{1-x}Sn_xSe/Pb_{1-y}Eu_ySe$ MQWs, the barrier material causes a nonnegligible anisotropy of the matrix elements.

The barrier material $Pb_{1-x}Eu_ySe$ is known to be anisotropic [35,36]. The electron velocities in $Pb_{1-y}Eu_ySe$ were determined as $v_\parallel \sim 6.5 \times 10^5$ m/s and $v_z = 3\sqrt{1.96-3.98 y} \times 10^5$ m/s,

so that $v_z \sim 3.75 \times 10^5$ m/s for MQW-39 ($K = v_\parallel/v_z = 1.73$) and $v_z \sim 4.00 \times 10^5$ m/s for MQW-25 ($K = 1.63$) [35]. The confined states in $Pb_{1-x}Sn_xSe/Pb_{1-y}Eu_ySe$ MQWs, whose probability densities in the barrier are nonnegligible (calculated to be 16 and 10% for MQW-25 and MQW-39, respectively), are sensitive to this anisotropy. TIS and TIS', whose probability densities are peaked at the interface [see Fig. 1(c)], are expected to be significantly influenced by the valley anisotropy of the barriers. The electron velocities of the QW states $v_\parallel$ and $v_z$ depict averaged velocities between those of the QW and barrier materials. They contain the effect of the barrier anisotropy and are treated as fitting parameters.

We use our theoretical model with $v_\parallel \neq v_z$ to interpret the magnetooptical data of MQW-39 shown in Fig. 2(b) and replotted in Fig. 4(a), where the calculated intersubband transitions between Landau levels from longitudinal ($\theta = 0^\circ$) and oblique valleys ($\theta = 70.5^\circ$) are shown by dashed and solid lines, respectively. Magnetooptical transitions occur between electronlike and holelike confined states of opposite parities; thus, the allowed transitions involve, for instance, TIS' and TIS, H2 and E2, ... Otherwise, the selection rules are like that obtained in the bulk case ($n \rightarrow n \pm 1$ and a spin flip) in the Faraday geometry. The absorptions shown by red dots are attributed to transitions between TIS' and TIS, and those in black dots correspond to transitions involving excited states (H2, E2, E3, ...) [see Figs. 3(a) and 3(b)]. The observed magnetooptical absorptions mainly correspond to oblique valley transitions, as they are three times more numerous than the longitudinal valley.
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Derived constant energy contours of topological interface states (TISs) at $\bar{M}$ (in red) and $\bar{\Gamma}$ (in black) for arbitrarily chosen $E = 13$ meV and multiquantum well (MQW)-39. It is directly deduced from the experimental observation. (b) Map of the two-dimensional (2D) Brillouin zone with the calculated dispersions of the MQW-39 TIS and TIS' near the $\bar{M}$ and $\bar{\Gamma}$ points.</figcaption></figure>

The oblique valley transitions nicely fit the nonmonotonic behavior of the Landau level transitions observed in the MQW sample, which is unequivocally due to the tilt and anisotropy of the oblique valleys that induce Landau level anticrossings. The unusual dispersions of the absorptions (red and black dots) vs magnetic field are clearly well-described by considering the interactions between TIS and massive states, as elaborated in Sec. III. For instance, the second red absorption line, which involves TIS' and TIS, has an oscillatory behavior vs magnetic field. This is due to the interactions of the Landau levels of the TIS with the ones of E2, then E3, etc., as the magnetic field is increased. These interactions not only change the transition energies vs magnetic field but also yield a relaxation of the selection rules near the Landau level anticrossings. To illustrate this effect, Fig. 4(b) highlights one of these interactions occurring at around $B = 20$ T and photon energies ~180 meV. The two red arrows account for the swing of the intensity between two absorption lines, which is typical of Landau level anticrossings [51,60]. Note that it is difficult to model quantitatively the absorption intensities as a function of the magnetic field, as the large number of involved Landau levels is strongly mixed between subbands, spins, and their indices. Therefore, our model provides very accurate transition energies, but it is rather difficult to provide a quantitative description of their intensities. These anticrossings are clear experimental manifestations and proofs of the anisotropic and tilted valleys found in the investigated multivalley MQWs.

The parameters used to fit the magnetooptical absorptions in MQW-39 are listed in Table II. The slight discrepancy of the gap values between the well material and Pb$_{0.75}$Sn$_{0.25}$Se thick film can be explained by small Sn content variations or/and residual strain imposed by the barriers in the MQWs. The velocities are found to be slightly different for both types of valleys. Note that the anisotropy induced by the barriers has no reason to be similar for the two types of valleys [59,61]. Indeed, the penetration depth of the envelope functions depends on $v_z$, which is valley dependent [43].

Our analysis and the parameters determined in Table II lead to the accurate determination of the TIS anisotropy in

the XY plane [12,24]. Because of the anisotropy and different tilts of the 3D ellipsoids [see Fig. 1(a)], their projection on the 2D BZ are not equivalent at the $\bar{M}$ and $\bar{\Gamma}$ points, leading to anisotropic constant energy contours, as sketched in Fig. 5(a). Remarkably, this projection due to confinement allows for a different anisotropy factor at $\bar{\Gamma}$ and $\bar{M}$, which is not possible with bulk materials where the anisotropy factor $K$ is identical for oblique and longitudinal valleys. We thus define the two anisotropy factors $K_{\bar{M}}$ and $K_{\bar{\Gamma}}$ associated with the constant energy surfaces of the TIS located at $\bar{M}$ and $\bar{\Gamma}$. They are determined by the ratios between the in-plane velocities in the Y and X directions and are listed in Table II. The coordinate rotation detailed in Appendix A defines $\bar{M}-\bar{\Gamma}$ as the X direction and $\bar{M}-\bar{K}$ as the Y axis [see Fig. 1(a)]. The expression for the TIS dispersions is given in Appendix C. For the longitudinal valley, we obtain an isotropic dispersion in the XY plane, with the velocity $v_\parallel$. Therefore, one gets $K_{\bar{\Gamma}} = 1$ at the $\bar{\Gamma}$ point, which follows directly from the isotropy of the well and barrier materials in the layer plane [see Eq. (1)] for this geometry. However, the TIS dispersions near the $\bar{M}$ points are anisotropic with electron velocities $v_\parallel$ along $\bar{M}-\bar{K}$ and $\cos^2\theta v_\parallel + \sin^2\theta v_z$ in the $\bar{M}-\bar{\Gamma}$ directions (see Appendix C). The magnetooptical determination of $v_\parallel$ and $v_z$ by fitting the anticrossing magnitudes allows the deduction of the topological state dispersion anisotropy, which is given by

$$K_{\bar{M}} = \left( \cos^2\theta + \frac{v_z^o \sin^2\theta}{v_\parallel^o} \right)^{-1}.$$

For the sample MQW-39, we find $K_{\bar{M}} = 1.23$, which yields the constant energy surface shown in red in Fig. 5(a) at 13 meV. For comparison, the isotropic constant energy surface of the longitudinal TIS is drawn in black. The TIS dispersions near the $\bar{M}$ and $\bar{\Gamma}$ points are shown in Fig. 5(b). These data demonstrate that magnetospectroscopy is a powerful technique to determine accurately the anisotropy of the massive or massless Dirac cones (TIS and TIS') at the $\bar{M}$ points.

The magnetooptical spectra up to 34 T obtained on the second MQW sample (MQW-25) are shown in Fig. 6(a), and the
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Magnetooptical spectra of multiquantum well (MQW)-25. Some transmission minima are marked by colored arrows. (b) Magnetooptical fan chart of MQW-25. Dots represent transmission minima pointed out in (a). Those in red represent absorptions mainly involving topological interface state (TIS) and TIS', while the black dots are mainly due to massive states. Dashed (solid) lines are the fit from the Landau level transitions emerging from the longitudinal valley (oblique valleys). The transition labels are given for the longitudinal valley.</figcaption></figure>

interpretation is summarized by Fig. 6(b). The solid lines are transitions between Landau levels of the oblique valleys, and the dashed lines account for the longitudinal one. The fitting parameters are listed in Table II and are found to be like those of MQW-39. Again, our theoretical model remarkably accounts for the observed magnetooptical transitions, including several anticrossings between TIS and massive state Landau levels emerging at the $\tilde{M}$ points. For this sample MQW-25, the anisotropy factor for the TIS dispersion at the $\tilde{M}$ point is deduced as $K_{\tilde{M}} = 1.17$ and is found to be slightly lower than the one of MQW-39. These results confirm our analysis and fully demonstrate the interaction of TIS Landau levels with those of the massive states in the MQWs in the presence of valley tilt and anisotropy.

V. CONCLUSIONS

We have demonstrated significant interactions between topological and massive states in multivalley TCI QW heterostructures. The anisotropy induced by the barrier normal insulator allows us to differentiate between the electronic properties of the topological states emerging from the longitudinal and oblique valleys. The Landau levels associated to the oblique valleys show anticrossings between TISs and massive states near the $\tilde{M}$ points of the 2D BZ, which is not the case for Landau levels emerging from the longitudinal valley, near the $\tilde{\Gamma}$ point. Moreover, our work allows us to accurately determine the anisotropy of the topological states of both kinds of valleys. In this way, we demonstrate an isotropic, i.e., circular dispersion for TIS near $\tilde{\Gamma}$ and a

tunable elliptical anisotropy of TIS near the $\tilde{M}$ points induced by the strong valley anisotropy within the barrier. In this paper, we establish a link between the parent 3D ellipsoids of bulk materials and their projection when a quantum confinement is superimposed and thus opens the door to valley and anisotropy engineering of topological states. The effects demonstrated in this paper (Landau levels avoided crossings and TIS anisotropy) should also be observed in $Pb_{1-x}Sn_xTe$-based heterostructures, where the anisotropy of the TIS should be greatly enhanced, as $Pb_{1-x}Sn_xTe$ is much more anisotropic ($K \sim 3.2$) than $Pb_{1-x}Sn_xSe$ [52,60,62].

In this paper, we thus pave the way to valley engineering in topological insulators and to future investigations showing valley effects that can spontaneously break crystalline symmetry, such as nematic valley ordering. The anisotropy of the TIS dispersions in $Pb_{1-x}Sn_xSe$ QWs is indeed predicted to promote a valley polarization that would drive the formation of quantum Hall ferroelectric states as theoretically anticipated for different material systems [15–17,63,64].

ACKNOWLEDGMENTS

The authors want to thank G. Bastard for enlightening discussions. This paper is supported by a coproject between N° ANR-19-CE30-022-01 for ENS and the Austrian Science Fund FWF, Project I-4493 for JKU. We acknowledge the support of LNCMI-CNRS, a member of the EMFL. M.O. acknowledges the support from ANR No-19-CE30-0032 project Collector.
## APPENDIX A: OBLIQUE VALLEY HAMILTONIAN IN HETEROSTRUCTURE

As pointed out in Sec. III B, for an oblique valley with main axis (x, y, z), the quantum operators have to be expressed in a new coordinate system (X, Y, Z) where Z//B//[111]. The relation between operators in both coordinate systems is encoded in a rotation matrix M, which leads to p = MP and σ = MΣ for the momentum and spin operators, respectively [55]. Using the invariance of a scalar product under rotation, the off-diagonal terms of the Hamiltonian in Eq. (1) transform as

$$v_{\parallel} \mathbf{p}_{\parallel} \boldsymbol{\sigma}_{\parallel} + v_z p_Z \boldsymbol{\sigma}_Z \rightarrow v_{\parallel} \mathbf{P}_{\parallel} \boldsymbol{\Sigma}_{\parallel} + v_z P_Z \boldsymbol{\Sigma}_Z + (v_{\parallel} - v_z) \times (P_Z \boldsymbol{\Sigma}_Z - p_z \boldsymbol{\sigma}_Z),$$

with

$$P_Z \Sigma_Z - p_z \sigma_z = m_{zx} [(m_{zx} P_Z - m_{zz} P_X) \Sigma_Z - (m_{zx} P_X + m_{zz} P_Z) \Sigma_X],$$

$$\left[ \begin{array}{l} -\frac{P_{\parallel}^2}{2\bar{m}} \mathbf{1} - \frac{\sin\theta}{2\bar{m}} h(X, Z) \\ \\ v_{\parallel} \mathbf{P}_{\parallel} \cdot \boldsymbol{\Sigma}_{\parallel} + \sin\theta(v_{\parallel} - v_z) \mathcal{H}(X, Z) \end{array} \right]$$

with $h(X, Z)$ and $\mathcal{H}(X, Z)$ given by Eqs. (5) and (6), respectively. Under a magnetic field along the Z direction, we can write

$$\begin{aligned} P_Z &= -\frac{i\hbar d}{dZ} \\ P_X &= \sqrt{\frac{e\hbar B}{2}}(a + a^+) \\ P_X + iP_Y &= \sqrt{2e\hbar Ba^+} \\ P_X - iP_Y &= \sqrt{2e\hbar Ba}, \end{aligned}$$

where $a$ and $a^+$ are the ladder operators. The perturbative Hamiltonian cannot be projected into a trivial basis of harmonic oscillator functions, as it was the case for the longitudinal valley. The envelope function needs to be expressed in its general form, and for the jth confined state, it is given by

$$\psi^{(j)} = \sum_{n,i} \left\{ \alpha_{n,i} \begin{bmatrix} F_{+}^{(j)} \\ 0 \\ F_{-}^{(j)} \\ 0 \end{bmatrix} + \beta_{n,i} \begin{bmatrix} 0 \\ F_{+}^{(j)} \\ -F_{-}^{(j)} \\ 0 \end{bmatrix} \right\} |n\rangle,$$

where $|n\rangle$ is the harmonic oscillator function with $n = 0, 1, 2, \dots$, and $\alpha$ and $\beta$ are the coefficients in front of the

where $m_{\alpha\beta}$ denote the matrix elements of M. The three oblique valleys are tilted by an angle $\theta$ with respect to Z [see Fig. 1(a)] so that $m_{zz} = \cos\theta$. As $MM^{-1} = 1$, one gets $m_{zx}^2 + m_{zy}^2 + m_{zz}^2 = 1$. The isotropic dispersion in the xy plane allows us to arbitrarily choose $m_{zy} = 0$ (i.e., choosing a rotation with respect to y//[-10]), which leads to $m_{zx} = |\sin\theta|$. One ends up with $P_Z \Sigma_Z - p_z \sigma_z = \sin\theta \mathcal{H}(X, Z)$ as defined in Eq. (6). A similar transformation is done for the diagonal terms $(p_{\parallel}^2/2\bar{m})\mathbf{1}$ of (1):

$$\frac{P_{\parallel}^2}{2\bar{m}} \mathbf{1} \rightarrow \frac{P_{\parallel}^2}{2\bar{m}} \mathbf{1} + \frac{1}{2\bar{m}} (P_Z^2 - p_z^2) \mathbf{1},$$

where $p_z = m_{zx}P_X + m_{zz}P_Z$. The resulting Hamiltonian is given by Eq. (4) as $(P_Z^2 - p_z^2)\mathbf{1} = \sin\theta h(X, Z)$ [Eq. (5)].

## APPENDIX B: OBLIQUE VALLEY LANDAU LEVELS

The perturbative terms for the oblique valleys [second and third terms of Eq. (4)] are given by

$$v_{\parallel} \mathbf{P}_{\parallel} \cdot \boldsymbol{\Sigma}_{\parallel} + \sin\theta(v_{\parallel} - v_z) \mathcal{H}(X, Z) \begin{bmatrix} \frac{P_{\parallel}^2}{2\bar{m}} \mathbf{1} + \frac{\sin\theta}{2\bar{m}} h(X, Z) \end{bmatrix},$$

two spin components. The perturbative theory requires us to calculate matrix elements $\langle\psi^{(i)}|\delta W|\psi^{(j)}\rangle$, where $1 \le i, j \le N$ with $N$ the number of confined states obtained at $k_{\parallel} = 0$. This leads to two equations on $\alpha_{n,j}$ and $\beta_{n,j}$ that are coupled with $\alpha_{n\pm1,j}, \beta_{n\pm1,j}, \alpha_{n\pm2,j}, \dots$. We solved these equations by taking $1 \le j \le N$ and $0 \le n \le 8$, which insures a convergence. For $n = 0$, only the equation on $\beta_{0,j}$ is considered, the equation on $\alpha_{0,j}$ being unphysical in this system. Indeed, the two $n = 0$ Landau levels are spin polarized [50,51]. This resolution leads to the oblique valley Landau levels shown in Fig. 3(b). Moreover, we can extract from these equations the two coupling potentials $W_1$ and $W_2$. They occur between $\beta_{n+1,j}$ and $\beta_{n,j \pm (2q+1)}$ [or between $\alpha_{n+1,j}$ and $\alpha_{n,j \pm (2q+1)}$] for $W_1$ and between $\beta_{n,j}$ and $\alpha_{n+1,j \pm (2q+2)}$ for $W_2$ ($q \in \mathbb{N}$).

## APPENDIX C: ANISOTROPY OF THE TISs

The theoretical model detailed in Sec. III B gives the perturbative matrix for the in-plane dispersion calculations. The resolution is like the one in Appendix B.

As an illustration, in the subspace TIS and TIS' (twofold degenerated due to the spin), the perturbation theory leads to the diagonalization of the following matrix (for the sake of simplicity, small terms like remote band effects are neglected):

$$\left\{ \begin{array}{l} -\delta \Sigma_Z \\ i[\hbar v_{\parallel} k_{+} - \hbar \sin^2 \theta (v_{\parallel} - v_z) k_X] \Sigma_X \\ -i[\hbar v_{\parallel} k_{-} - \hbar \sin^2 \theta (v_{\parallel} - v_z) k_X] \Sigma_X \end{array} \right\}.$$
Here, $\pm\delta$ are the energies of TIS' and TIS calculated at $k_\parallel = 0$ and $k_\perp = k_X \pm ik_Y$. In a good approximation, the TIS dispersions are thus twofold degenerated with

$$E(\theta) = \pm \sqrt{\delta^2 + \hbar^2 v_\parallel^2 k_Y^2 + \hbar^2 [v_\parallel - \sin^2\theta(v_\parallel - v_z)]^2 k_X^2}$$

For the longitudinal valley, we have $\theta = 0^\circ$ so that the intersubband matrix element of momentum between TIS and TIS' is $v_\parallel$. At the $\bar{M}$ point, where $\theta = 70.5^\circ$, one immediately obtains the in-plane velocities and the anisotropy factor $K_{\bar{M}}$ mentioned in Sec. IV.

As a final remark, when a magnetic field along Z is applied, the Hamiltonian above gives, in a good approximation,

the Landau levels of TIS' and TIS for the two types of valley:

$$E_n(\theta) = \pm \sqrt{\delta^2 + 2e\hbar B n \left[ \frac{v_\parallel + v_z}{2} + \frac{\cos^2\theta}{2} (v_\parallel - v_z) \right]^2}$$

Therefore, the in-plane (XY plane) motions of the electrons under magnetic field in both valleys are characterized by the following velocity:

$$\left[ \frac{v_{\parallel} + v_z}{2} + \frac{\cos^2\theta}{2} (v_{\parallel} - v_z) \right]$$

[1] C. E. Nebel, Electrons dance in diamond, Nat. Mater. 12, 690 (2013).

[2] J. R. Schaibley, H. Yu, G. Clark, P. Rivera, J. S. Ross, K. L. Seyler, W. Yao, and X. Xu, Valleytronics in 2D materials, Nat. Rev. Mater. 1, 16055 (2016).

[3] F. Bussolotti, H. Kawai, Z. E. Ooi, V. Chellappan, D. Thian, A. L. C. Pang, and K. E. J. Goh, Roadmap on finding chiral valleys: screening 2D materials for valleytronics, Nano Futur. 2, 32001 (2018).

[4] M. Sun, I. G. Savenko, H. Flayac, and T. C. H. Liew, Multi-valley engineering in semiconductor microcavities, Sci. Rep. 7, 45243 (2017).

[5] J. Isberg, M. Gabrysch, J. Hammersberg, S. Majdi, K. Kovi, and D. J. Twitchen, Generation, transport and detection of valley-polarized electrons in diamond, Nat. Mater. 12, 760 (2013).

[6] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Silicon quantum electronics, Rev. Mod. Phys. 85, 961 (2013).

[7] L. Zhang, J.-W. Luo, A. Saraiva, B. Koiller, and A. Zunger, Genetic design of enhanced valley splitting towards a spin qubit in silicon, Nat. Commun. 4, 2396 (2013).

[8] S. Goswami, K. A. Slinker, M. Friesen, L. M. McGuire, J. L. Truitt, C. Tahan, L. J. Klein, J. O. Chu, P. M. Mooney, D. W. van der Weide, R. Joynt, S. N. Coppersmith, and M. A. Eriksson, Controllable valley splitting in silicon quantum devices, Nat. Phys. 3, 41 (2007).

[9] R. C. C. Leon, C. H. Yang, J. C. C. Hwang, J. C. Lemyre, T. Tanttu, W. Huang, K. W. Chan, K. Y. Tan, F. E. Hudson, K. M. Itoh, A. Morello, A. Laucht, M. Pioro-Ladrière, A. Saraiva, and A. S. Dzurak, Coherent spin control of s-, p-, d- and f-electrons in a silicon quantum dot, Nat. Commun. 11, 797 (2020).

[10] C. H. Yang, R. C. C. Leon, J. C. C. Hwang, A. Saraiva, T. Tanttu, W. Huang, J. Camirand Lemyre, K. W. Chan, K. Y. Tan, F. E. Hudson, K. M. Itoh, A. Morello, M. Pioro-Ladrière, A. Laucht, and A. S. Dzurak, Operation of a silicon quantum processor unit cell above one kelvin, Nature (London) 580, 350 (2020).

[11] A. Rycerz, J. Tworzydło, and C. W. J. Beenakker, Valley filter and valley valve in graphene, Nat. Phys. 3, 172 (2007).

[12] O. Gunawan, Y. P. Shkolnikov, K. Vakili, T. Gokmen, E. P. De Poortere, and M. Shayegan, Valley Susceptibility of an Interacting Two-Dimensional Electron System, Phys. Rev. Lett. 97, 186404 (2006).

[13] G. Y. Wu, N.-Y. Lue, and L. Chang, Graphene quantum dots for valley-based quantum computing: a feasibility study, Phys. Rev. B 84, 195463 (2011).

[14] N. Rohling and G. Burkard, Universal quantum computing with spin and valley states, New J. Phys. 14, 83008 (2012).

[15] I. Sodemann, Z. Zhu, and L. Fu, Quantum Hall Ferroelectrics and Nematics in Multivalley Systems, Phys. Rev. X 7, 041068 (2017).

[16] P. M. Tam, T. Liu, I. Sodemann, and L. Fu, Local probes for quantum Hall ferroelectrics and nematics, Phys. Rev. B 101, 241103 (2020).

[17] X. Li, F. Zhang, and A. H. MacDonald, SU(3) Quantum Hall Ferromagnetism in SnTe, Phys. Rev. Lett. 116, 026803 (2016).

[18] B. E. Feldman, M. T. Randeria, A. Gyenis, F. Wu, H. Ji, R. J. Cava, A. H. MacDonald, and A. Yazdani, Observation of a nematic quantum Hall liquid on the surface of bismuth, Science 354, 316 (2016).

[19] D. L. Mitchell and R. F. Wallis, Theoretical energy-band parameters for the lead salts, Phys. Rev. 151, 581 (1966).

[20] C. M. Polley, P. Dziawa, A. Reszka, A. Szczerbakow, R. Minikayev, J. Z. Domagala, S. Safaei, P. Kacman, R. Buczko, J. Adell, M. H. Berntsen, B. M. Wojek, O. Tjernberg, B. J. Kowalski, T. Story, and T. Balasubramanian, Observation of topological crystalline insulator surface states on (111)-oriented Pb$_{1-x}$Sn$_x$Se films, Phys. Rev. B 89, 075317 (2014).

[21] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).

[22] L. Fu, Topological Crystalline Insulators, Phys. Rev. Lett. 106, 106802 (2011).

[23] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Topological crystalline insulators in the SnTe material class, Nat. Commun. 3, 982 (2012).

[24] P. Dziawa, B. J. Kowalski, K. Dybko, R. Buczko, A. Szczerbakow, M. Szot, E. Łusakowska, T. Balasubramanian, B. M. Wojek, M. H. Berntsen, O. Tjernberg, and T. Story, Topological crystalline insulator states in Pb$_{1-x}$Sn$_x$Se, Nat. Mater. 11, 1023 (2012).

[25] B. M. Wojek, P. Dziawa, B. J. Kowalski, A. Szczerbakow, A. M. Black-Schaffer, M. H. Berntsen, T. Balasubramanian, T. Story, and O. Tjernberg, Band inversion and the topological phase transition in (Pb,Sn)Se, Phys. Rev. B 90, 161202 (2014).

[26] S.-Y. Xu, C. Liu, N. Alidoust, M. Neupane, D. Qian, I. Belopolski, J. D. Denlinger, Y. J. Wang, H. Lin, L. A. Wray, G. Landolt, B. Slomski, J. H. Dil, A. Marcinkova, E. Morosan, Q.
Gibson, R. Sankar, F. C. Chou, R. J. Cava, A. Bansil, and M. Z. Hasan, Observation of a topological crystalline insulator phase and topological phase transition in Pb<sub>1-x</sub>Sn<sub>x</sub>Te, Nat. Commun. **3**, 1192 (2012).

[27] Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T. Takahashi, K. Segawa, and Y. Ando, Experimental realization of a topological crystalline insulator in SnTe, Nat. Phys. **8**, 800 (2012).

[28] J. O. Dimmock and G. B. Wright, Band edge structure of PbS, PbSe, and PbTe, Phys. Rev. **135**, A821 (1964).

[29] G. Martinez, Band inversion in Pb<sub>1-x</sub>Sn<sub>x</sub>Se alloys under hydrostatic pressure. II. Galvanomagnetic properties, Phys. Rev. B **8**, 4686 (1973).

[30] J. O. Dimmock, **k·p** theory for the conduction and valence bands of PbSnTe and PbSnSe alloys, in *Proceedings of the International Conference of the Physics of Semimetals and Narrow Gap Semiconductor* (Pergamon Press, New York, 1971), pp. 319–330.

[31] Y. Ando and L. Fu, Topological crystalline insulators and topological superconductors: from concepts to materials, Annu. Rev. Condens. Matter Phys. **6**, 361 (2015).

[32] G. Krizman, B. A. Assaf, T. Phuphachong, G. Bauer, G. Springholz, L. A. de Vaulchier, and Y. Guldner, Dirac parameters and topological phase diagram of Pb<sub>1-x</sub>Sn<sub>x</sub>Se from magnetospectroscopy, Phys. Rev. B **98**, 245202 (2018).

[33] J. Liu, W. Duan, and L. Fu, Two types of surface states in topological crystalline insulators, Phys. Rev. B **88**, 241303 (2013).

[34] M. Serbyn and L. Fu, Symmetry breaking and Landau quantization in topological crystalline insulators, Phys. Rev. B **90**, 035402 (2014).

[35] M. Simma, G. Bauer, and G. Springholz, Temperature dependent band offsets in PbSe/PbEuSe quantum well heterostructures, Appl. Phys. Lett. **101**, 172106 (2012).

[36] P. Rothlein, H. Pascher, G. Bauer, and M. Tacke, Coherent Raman scattering and magneto-optical interband transitions in diluted magnetic IV–VI compounds, Semicond. Sci. Technol. **5**, S147 (1990).

[37] B. A. Assaf, T. Phuphachong, V. V. Volobuev, A. Inhofer, G. Bauer, G. Springholz, L. A. de Vaulchier, and Y. Guldner, Massive and massless Dirac fermions in Pb<sub>1-x</sub>Sn<sub>x</sub>Te topological crystalline insulator probed by magneto-optical absorption, Sci. Rep. **6**, 20323 (2016).

[38] H. Burkhard, G. Bauer, and W. Zawadzki, Band-population effects and intraband magneto-optical properties of a many-valley semiconductor: PbTe, Phys. Rev. B **19**, 5149 (1979).

[39] Y. Okada, M. Serbyn, H. Lin, D. Walkup, W. Zhou, C. Dhital, M. Neupane, S. Xu, Y. J. Wang, R. Sankar, F. Chou, A. Bansil, M. Z. Hasan, S. D. Wilson, L. Fu, and V. Madhavan, Observation of Dirac node formation and mass acquisition in a topological crystalline insulator, Science **341**, 1496 (2013).

[40] G. Bauer, H. Pascher, and W. Zawadzki, Magneto-optical properties of semimagnetic lead chalcogenides, Semicond. Sci. Technol. **7**, 703 (1992).

[41] M. S. Adler, C. R. Hewes, and S. D. Senturia, **k·p** model for the magnetic energy levels in PbTe and Pb<sub>1-x</sub>Sn<sub>x</sub>Te, Phys. Rev. B **7**, 5186 (1973).

[42] M. Simma, T. Fromherz, G. Bauer, and G. Springholz, Type I/type II band alignment transition in strained PbSe/PbEuSeTe multiquantum wells, Appl. Phys. Lett. **95**, 212103 (2009).

[43] G. Krizman, B. A. Assaf, T. Phuphachong, G. Bauer, G. Springholz, G. Bastard, R. Ferreira, L. A. de Vaulchier, and Y. Guldner, Tunable Dirac interface states in topological superlattices, Phys. Rev. B **98**, 075303 (2018).

[44] B. Zhou, H.-Z. Lu, R.-L. Chu, S.-Q. Shen, and Q. Niu, Finite Size Effects on Helical Edge States in a Quantum Spin-Hall System, Phys. Rev. Lett. **101**, 246807 (2008).

[45] Y. Zhang, K. He, C. Z. Chang, C. L. Song, L. L. Wang, X. Chen, J. F. Jia, Z. Fang, X. Dai, W. Y. Shan, S. Q. Shen, Q. Niu, X. L. Qi, S. C. Zhang, X. C. Ma, and Q. K. Xue, Crossover of the three-dimensional topological insulator Bi<sub>2</sub>Se<sub>3</sub> to the two-dimensional limit, Nat. Phys. **6**, 584 (2010).

[46] S. Safaei, P. Pacman, and R. Buczko, Topological crystalline insulator (Pb,Sn)Te: surface states and their spin polarization, Phys. Rev. B **88**, 045305 (2013).

[47] B. A. Assaf, T. Phuphachong, V. V. Volobuev, G. Bauer, G. Springholz, L. A. de Vaulchier, and Y. Guldner, Magnetooptical determination of a topological index, npj Quantum Mater. **2**, 26 (2017).

[48] A. J. Strauss, Inversion of conduction and valence bands in Pb<sub>1-x</sub>Sn<sub>x</sub>Se alloys, Phys. Rev. **157**, 608 (1967).

[49] A. R. Calawa, J. O. Dimmock, T. C. Harman, and I. Melngailis, Magnetic Field Dependence of Laser Emission in Pb<sub>1-x</sub>Sn<sub>x</sub>Se Diodes, Phys. Rev. Lett. **23**, 7 (1969).

[50] B. A. Assaf, T. Phuphachong, E. Kampert, V. V. Volobuev, P. S. Mandal, J. Sánchez-Barriga, O. Rader, G. Bauer, G. Springholz, L. A. de Vaulchier, and Y. Guldner, Negative Longitudinal Magnetoresistance from the Anomalous $N = 0$ Landau Level in Topological Materials, Phys. Rev. Lett. **119**, 106602 (2017).

[51] G. Krizman, B. A. Assaf, M. Orlita, T. Phuphachong, G. Bauer, G. Springholz, G. Bastard, R. Ferreira, L. A. de Vaulchier, and Y. Guldner, Avoided level crossing at the magnetic field induced topological phase transition due to spin-orbital mixing, Phys. Rev. B **98**, 161202 (2018).

[52] G. Bauer, Magnetooptical properties of IV–VI compounds, in *Narrow Gap Semiconductors Physics and Applications*, Lecture Notes in Physics, edited by W. Zawadzki (Springer, Berlin, Heidelberg, 1980), Vol. 133, pp. 427–446.

[53] K. K. Tikuišis, J. Wyzula, L. Ohnoutek, P. Cejpek, K. Uhlířová, M. Hakl, C. Faugeras, K. Výborný, A. Ishida, M. Veis, and M. Orlita, Landau Level spectroscopy of the PbSnSe topological crystalline insulator, Phys. Rev. B **103**, 155304 (2021).

[54] H. Pascher, G. Bauer, and R. Grisar, Magnetooptical investigations and four-wave-mixing spectroscopy of PbSe, Phys. Rev. B **38**, 3383 (1988).

[55] A. C. Bartnik, A. L. Efros, W.-K. Koh, C. B. Murray, and F. W. Wise, Electronic states and optical properties of PbSe nanorods and nanowires, Phys. Rev. B **82**, 195313 (2010).

[56] G. Bastard, Theoretical investigations of superlattice band structure in the envelope-function approximation, Phys. Rev. B **25**, 7584 (1982).

[57] G. Bastard, Superlattice Band Structure in the Envelope-Function Approximation, Phys. Rev. B **24**, 5693 (1981).

[58] G. Bastard, *Wave Mechanics Applied to Semiconductor Heterostructures* (Les Editions de Physique, Les Ulis, 1988).

[59] E. A. de Andrada e Silva, Optical transition energies for leadsalt semiconductor quantum wells, Phys. Rev. B **60**, 8859 (1999).
[60] D. C. Tsui, G. Kaminsky, and P. H. Schmidt, Tunneling study of surface quantization in PbTe, Phys. Rev. B 9, 3524 (1974).

[61] R. Buczko and Ł. Cywiński, PbTe/PbSnTe heterostructures as analogs of topological insulators, Phys. Rev. B 85, 205319 (2012).

[62] G. E. Tudury, M. V. Marquezini, L. G. Ferreira, L. C. Barbosa, and C. L. Cesar, Effect of band anisotropy on electronic structure of PbS, PbSe, and PbTe quantum dots, Phys. Rev. B 62, 7357 (2000).

[63] D. A. Abanin, S. A. Parameswaran, S. A. Kivelson, and S. L. Sondhi, Nematic valley ordering in quantum Hall systems, Phys. Rev. B 82, 035428 (2010).

[64] S. A. Parameswaran and B. E. Feldman, Quantum Hall valley nematics, J. Phys.: Condens. Matter 31, 273001 (2019).