Inversion symmetry breaking in the probability density by surface-bulk hybridization in topological insulators

Topological insulators (TIs) have emerged as fascinating materials with unique electronic properties that hold great promise for potential applications in electronic transport, spintronics Hasan and Kane (2010) and quantum computation Fu and Kane (2008). Characterized by a large band-gap in the bulk, they possess topologically protected surface states with pseudo-relativistic (Dirac) dispersion that confer them non-trivial conducting properties Hasan and Kane (2010). The interplay between the bulk and surface states in TIs plays a pivotal role in shaping their behavior and functionality Moore (2010).

Previous studies have focused on the properties of surface states in topological insulators, providing valuable insights into their behavior. Extensive exploration has been carried out on confined electronic states and optical transitions in nanoparticles with different geometries Castaño-Yepes and Muñoz (2022); Imura et al. (2012); Tian et al. (2013); Governale et al. (2020). From both theoretical and experimental perspectives, the spin properties of surface states have garnered attention in the field of the so-called quantum spin-Hall effect Zhou et al. (2008); Lu et al. (2010). Furthermore, there is evidence suggesting the existence of residual bulk carriers that couple with surface states through disorder Saha and Garate (2014); Black-Schaffer and Balatsky (2012); Mandal et al. (2023); Yang et al. (2020). On the other hand, it has been observed that charge carriers in TI nanostructures consist of both surface and bulk electronic states Checkelsky et al. (2011); Chiatti et al. (2016); Liao et al. (2017); Steinberg et al. (2011); Zhao et al. (2013). These studies point to the hybridization and interaction between topologically trivial and non-trivial states. The precise mechanism governing their coupling is not yet fully understood, but its comprehension is vital for elucidating the intricate nature of transport and photoemission experiments Li et al. (2015); Rosenbach et al. (2022); Singh et al. (2017). Earlier studies have described the surface to bulk hybridization mechanism via a Fano model, with a phenomenological coupling constant Hsu et al. (2014), combined with ab-initio calculations in a slab geometry. Those results show that, despite the hybridization with the bulk state, the spin polarization of the surface state tends to be preserved by this mechanism Hsu et al. (2014).

In this study, we investigate the hybridization between surface and bulk states in TIs, calculated from a paradigmatic continuous Hamiltonian model Liu et al. (2010); Zhang et al. (2009) in the geometry of a flat slab of finite thickness $d$. First principles calculations Reid et al. (2020) in this particular geometry point towards the relevance of the slab thickness on the size of the gap and for the existence of the surface state. Our focus lies in the observation of a significant asymmetry in the electronic density, resulting from the hybridization between bulk and surface states, and we highlight the role of parity and spin in this simple but rather surprising effect. Interestingly, this analytical result may provide an interpretation for the symmetry breaking effect observed in magnetotransport experiments in $\text{Bi}_{2}\text{Se}_{3}$ nanoplates on $\text{Fe}_{3}\text{O}_{4}$ substrates Buchenau et al. (2017).

To analyze the possibility for hybridization between surface and bulk states in a topological insulator, along with its physical consequences, we shall start by considering a generic, low-energy continuum model obtained and reported in the literature Liu et al. (2010); Zhang et al. (2009) by $(\mathbf{k}\cdot\mathbf{p})$-perturbation theory, that represents the band structure of a typical TI near the $\Gamma$-point Liu et al. (2010); Zhang et al. (2009):

	$\displaystyle\hat{\mathcal{H}}$	$\displaystyle=$	$\displaystyle\left(\frac{\Delta_{0}}{2}+\frac{\gamma}{k_{\Delta}}\mathbf{k}^{2}\right)\hat{\tau}_{z}\otimes\hat{\sigma}_{0}$		(1)
		$\displaystyle+$	$\displaystyle\gamma\left(k\hat{\tau}_{x}\otimes\hat{\sigma}_{z}+k_{-}\hat{\tau}_{x}\otimes\hat{\sigma}_{+}+k_{+}\hat{\tau}_{x}\otimes\hat{\sigma}_{-}\right),$

where the spin ($\hat{\sigma}$) and pseudo-spin ($\hat{\tau}$) degrees of freedom are incorporated. Here, $\mathbf{k}=(k_{x},k_{y},k_{z})=-\mathrm{i}\nabla$ represents the momentum operator, with $k_{\pm}=k_{x}\pm\mathrm{i}k_{y}$ combinations of its components, while $\hat{\sigma}_{\pm}=(\hat{\sigma}_{x}\pm\mathrm{i}\hat{\sigma}_{y})/2$ are the ladder spin operators. As a reference, representative numerical values for the parameters $\Delta_{0}$, $\gamma$, and $k_{\Delta}$ for three different TI materials are presented in Table 1.

In order to simplify the mathematical notation, and further algebraic manipulations, we define the “mass" differential operator

$\displaystyle\mathbb{M}(\mathbf{k})\equiv\frac{\Delta_{0}}{2}+\frac{\gamma}{k_{\Delta}}\mathbf{k}^{2},$

(2)

so that the Hamiltonian Eq. (1) can be written in the alternative and more compact form

$\displaystyle\hat{\mathcal{H}}=\mathbb{M}(\mathbf{k})\hat{\beta}+\gamma\left(\bm{\alpha}\cdot\mathbf{k}\right),$

(3)

where we defined the Dirac $4\times 4$ matrices in the standard representation

$\displaystyle\bm{\alpha}=\begin{pmatrix}0&\bm{\hat{\sigma}}\\ \bm{\hat{\sigma}}&0\end{pmatrix},~{}~{}~{}\hat{\beta}=\begin{pmatrix}\mathbb{I}&0\\ 0&-\mathbb{I}\end{pmatrix},$

(4)

with $\bm{\hat{\sigma}}=\left(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}\right)$. The representation Eq. (3) also provides the conceptual advantage to possess the explicit form of a pseudo-relativistic Dirac Hamiltonian with a momentum-dependent "mass" operator as defined by Eq. (2). One can thus make the appropriate analogies to interpret the physical meaning of the parameters involved, where $\Delta_{0}$ clearly represents the "mass" that defines the finite gap of the bulk spectrum at the $\Gamma$-point, while $\gamma=\hbar v_{F}$ is proportional to the Fermi velocity, and $k_{\Delta}$ controls the curvature of the small quadratic correction to the linear dispersion relation at such point. We remark that an effective continuum model, such as Eq. (1) or its equivalent form Eq. (3), constitutes a fair representation of the physics at length-scales larger than the lattice constant of the material. As shown in Table 1, the lattice constants for the two metal trichalcogenides Bi₂Se₃ and Bi₂Te₃ are on the order of $\sim 4$Å. Based on these reference values, a continuum model representation is fair for confinement length scales on the order of $d\gtrsim 5R_{0}$, where $R_{0}\equiv 2\gamma/|\Delta_{0}|$ is a material-dependent parameter whose pseudo-relativistic interpretation in the context of the Dirac Hamiltonian is the "Compton wavelength", i.e. a characteristic length-scale for the support of the spinor eigenstates.

We assume that our system is a planar slab of a topological insulator with thickness $d$, whose Hamiltonian in dimensionless form can be written as

$\displaystyle\widetilde{\mathcal{H}}=\left(1+\frac{s}{\mathcal{R}}\widetilde{\mathbf{k}}^{2}\right)\hat{\beta}+s(\mathbf{\widetilde{k}}\cdot\hat{\bm{\alpha}}),$

(5)

where $\widetilde{\mathcal{H}}\equiv 2\hat{\mathcal{H}}/\Delta_{0}$, $\widetilde{k}\equiv R_{0}\hat{k}$, $\mathcal{R}\equiv k_{\Delta}R_{0}$, and $s=\text{sign}(\Delta_{0})$. In this dimensionless representation of the parameters, the corresponding eigenvalue problem is stated by the equation

$\displaystyle\widetilde{\mathcal{H}}\hat{\Psi}(\mathbf{r}_{\perp},z)=\epsilon\hat{\Psi}(\mathbf{r}_{\perp},z).$

(6)

From the geometry of our system, the translational invariance on the plane $(x,y)=\mathbf{r}_{\perp}$ normal to the interfaces at $z=0$ and $z=d$, respectively, allows us for separation of variables in the eigenfunction

$\displaystyle\hat{\Psi}(\mathbf{r}_{\perp},z)=e^{\mathrm{i}\mathbf{k}_{\perp}\cdot\mathbf{r}_{\perp}}\Psi(z),$

(7)

where $\mathbf{k}_{\perp}=(k_{x},k_{y})$, and $\Psi(z)$ is a four-component spinor. This separation of variables leads to the secular equation

$\displaystyle\left[1+\frac{s}{\mathcal{R}}\left(\widetilde{\mathbf{k}}_{\perp}^{2}+k_{z}^{2}\right)\hat{\beta}+s(\mathbf{\widetilde{k}_{\perp}}\cdot\hat{\bm{\alpha}}_{\perp}+\widetilde{k}_{z}\hat{\alpha}_{z})\right]\Psi(z)=\epsilon\Psi(z),$

(8)

where we defined the dimensionless energy eigenvalues

$\displaystyle\epsilon\equiv\frac{E}{\Delta_{0}/2}.$

(9)

The 4-component spinor eigenstates can be expressed in terms of its bi-spinor components $\chi$ and $\eta$, i.e. $\Psi(z)=e^{\mathrm{i}k_{z}z}(\eta,\chi)^{T}$, such that the eigenvalue problem Eq. (6) spans a system of coupled algebraic equations

	$\displaystyle\left[1+\frac{s}{\mathcal{R}}\left(\widetilde{\mathbf{k}}_{\perp}^{2}+k_{z}^{2}\right)-\epsilon\right]\eta+s\left(\widetilde{\mathbf{k}}_{\perp}\cdot\hat{\bm{\sigma}}_{\perp}+\widetilde{k}_{z}\hat{\sigma}_{z}\right)\chi$	$\displaystyle=$	$\displaystyle 0$
	$\displaystyle\left[1+\frac{s}{\mathcal{R}}\left(\widetilde{\mathbf{k}}_{\perp}^{2}+k_{z}^{2}\right)+\epsilon\right]\chi-s\left(\widetilde{\mathbf{k}}_{\perp}\cdot\hat{\bm{\sigma}}_{\perp}+\widetilde{k}_{z}\hat{\sigma}_{z}\right)\eta$	$\displaystyle=$	$\displaystyle 0.$

The latter system of equations implies that the bi-spinors $\eta$ and $\chi$ are not independent. Indeed, for a generic $(\tilde{\mathbf{k}}_{\perp},\tilde{k}_{z})$, it is straightforward to obtain from the second equation in the system Eq. (LABEL:eq_Dirac_system)

$\displaystyle\chi=\frac{s\left(\widetilde{k}_{z}\hat{\sigma}_{z}+\widetilde{\mathbf{k}}_{\perp}\cdot\hat{\bm{\sigma}}_{\perp}\right)}{1+\frac{s}{\mathcal{R}}\left(\widetilde{k}_{z}^{2}+\widetilde{\mathbf{k}}_{\perp}^{2}\right)+\epsilon}\eta,$

(11)

and the corresponding energy eigenvalues are functions of $\tilde{k}_{z}^{2}+\widetilde{\mathbf{k}}_{\perp}^{2}$, as follows

$\displaystyle\epsilon^{2}=\widetilde{k}_{z}^{2}+\widetilde{\mathbf{k}}_{\perp}^{2}+\left[1+\frac{s}{\mathcal{R}}\left(\widetilde{k}_{z}^{2}+\widetilde{\mathbf{k}}_{\perp}^{2}\right)\right]^{2}.$

(12)

Clearly then, the contribution to the energy eigenvalues arising from the transverse degrees of freedom is a monotonically increasing function of $\mathbf{\widetilde{k}_{\perp}}^{2}$. Therefore, since we are exploring the consequences of spontaneous surface-bulk hybridization, we shall be interested in the lowest energy eigenvalues, and hence we choose $\mathbf{\widetilde{k}_{\perp}}=0$, with the corresponding consequence on the transverse plane-wave phase of the spinor function in Eq. (7). Therefore, by setting $\mathbf{\widetilde{k}_{\perp}}=0$ in Eq. (12), the general energy dispersion relation that we shall use in what follows is

$\displaystyle\epsilon^{2}=\widetilde{k}_{z}^{2}+\left(1+\frac{s}{\mathcal{R}}\widetilde{k}_{z}^{2}\right)^{2}.$

(13)

The surface states are characterized by an exponential localization near the boundaries $z=0$ and $z=d$, respectively. This behaviour is captured by choosing $\tilde{k}_{z}=\mathrm{i}\tilde{\lambda}$ a complex number. Under this ansatz, Eq. (13) admits four possible solutions for $\tilde{k}_{z}$, given by $\pm\tilde{\lambda}_{1}$, and $\pm\tilde{\lambda}_{2}$, so that:

$\displaystyle\widetilde{\lambda}^{2}_{\alpha}=\frac{\mathcal{R}^{2}}{2}\left[1+\frac{2s}{\mathcal{R}}+(-1)^{1+\alpha}\sqrt{1+\frac{4}{\mathcal{R}}\left(s+\frac{\epsilon^{2}}{\mathcal{R}}\right)}\right],$

(14)

where we used the dimensionless form $\widetilde{\lambda}=R_{0}\lambda$. Figure 1 shows how each $\lambda$ changes as a function of the material-dependent parameter $\mathcal{R}$. In particular, as can be trivially verified from Eq. (14), for surface states ($\epsilon\simeq 0$) the imaginary part of $\lambda_{1,2}$ vanishes when $\mathcal{R}>4$, while it remains present for smaller values.

Then, by solving the linear system Eq. (LABEL:eq_Dirac_system), we obtain four independent solutions (for $\alpha=1,2$ and $\widetilde{z}=z/R_{0}$)

	$\displaystyle\varphi_{\pm}^{(\alpha,\uparrow)}(\widetilde{z},\epsilon)$	$\displaystyle=$	$\displaystyle e^{\pm\widetilde{\lambda}_{\alpha}\widetilde{z}}\left(1,0,\pm\mathrm{i}A_{\alpha},0\right)^{T},$
	$\displaystyle\varphi_{\pm}^{(\alpha,\downarrow)}(\widetilde{z},\epsilon)$	$\displaystyle=$	$\displaystyle e^{\pm\widetilde{\lambda}_{\alpha}\widetilde{z}}\left(0,1,0,\mp\mathrm{i}A_{\alpha}\right)^{T},$		(15)

where we defined the parameters

$\displaystyle A_{\alpha}=\frac{\widetilde{\lambda}_{\alpha}}{1-\frac{s}{\mathcal{R}}\widetilde{\lambda}^{2}_{\alpha}+\epsilon}.$

(16)

The exponential decay length in Eq. (15) is thus a non-linear function of the energy eigenvalues after Eq. (14). In particular, for $\mathcal{R}\gg 1$ after Eq. (14) it is clear that $\tilde{\lambda}_{\alpha}\sim\mathcal{R}/\sqrt{2}$, i.e. the exponential decay of the surface states from each side of the slab is approximately determined by the magnitude of this parameter. Concerning the materials we present in Table 1, realistic values are $\mathcal{R}=0.169$ for ${\rm{Bi}}_{2}{\rm{Se}}_{3}$ and $\mathcal{R}=0.021$ for ${\rm{Bi}}_{s}{\rm{Te}}_{3}$, respectively, and hence the imaginary part of the $\lambda_{1,2}$ does not vanish in both cases. Moreover, based on these material specific examples, we shall focus on the analysis for representative values $\mathcal{R}\leq 1$.

We remark that these independent solutions are trivially orthogonal for opposite spin components, i.e. (for $\beta,\,\beta^{\prime}=\pm$, $\alpha,\,\alpha^{\prime}=1,2$)

$$\langle\alpha\uparrow;\beta|\alpha^{\prime}\downarrow;\beta^{\prime}\rangle=\int_{0}^{\tilde{d}}d\tilde{z}\left[\varphi^{(\alpha,\uparrow)}_{\beta}(\tilde{z})\right]^{\dagger}\varphi^{(\alpha^{\prime},\downarrow)}_{\beta^{\prime}}(\tilde{z})=0,$$

(17)

while those with parallel spin components satisfy the inner product result

	$\displaystyle\langle\alpha\downarrow;\beta|\alpha^{\prime}\downarrow;\beta^{\prime}\rangle$	$\displaystyle=$	$\displaystyle\langle\alpha\uparrow;\beta|\alpha^{\prime}\uparrow;\beta^{\prime}\rangle$		(18)
		$\displaystyle=$	$\displaystyle\tilde{d}\left(1+\beta\beta^{\prime}A_{\alpha}A_{\alpha^{\prime}}\right)e^{\left(\beta\tilde{\lambda}_{\alpha}+\beta^{\prime}\tilde{\lambda}_{\alpha^{\prime}}\right)\tilde{d}/2}$
			$\displaystyle\times\frac{\sinh\left[\left(\beta\tilde{\lambda}_{\alpha}+\beta^{\prime}\tilde{\lambda}_{\alpha^{\prime}}\right)\tilde{d}/2\right]}{\left(\beta\tilde{\lambda}_{\alpha}+\beta^{\prime}\tilde{\lambda}_{\alpha^{\prime}}\right)\tilde{d}/2}.$

Therefore, it is natural to define independent surface states for each fixed spin projection $\sigma=\uparrow$ or $\sigma=\downarrow$ as follows

$\displaystyle\Psi_{\text{surface}}^{(\sigma)}(\widetilde{z})=\sum_{\alpha=1,2}\sum_{\beta=\pm}a_{\beta}^{(\alpha)}\varphi_{\beta}^{(\alpha,\sigma)}(\widetilde{z}),$

(19)

where $a_{\beta}^{(\alpha)}$ represent arbitrary coefficients in the general linear combination. Those coefficients are further defined upon imposing, for each spin projection state $\sigma=\uparrow$ or $\sigma=\downarrow$, the boundary conditions at each surface of the slab,

$\displaystyle\Psi_{\text{surface}}^{(\sigma)}(\widetilde{z}=0)=\Psi_{\text{surface}}^{(\sigma)}(\widetilde{z}=d/R_{0})=0.$

(20)

Those lead to the secular linear system

$\displaystyle\mathbb{A}~{}\bm{a}=0,$

(21)

where $\bm{a}\equiv\left(a_{+}^{(1)},a_{-}^{(1)},a_{+}^{(2)},a_{-}^{(2)}\right)^{T}$ is the vector of coefficients, and the matrix $\mathbb{A}$ is defined as (for $\widetilde{d}\equiv d/R_{0}$)

$\displaystyle\mathbb{A}\equiv\begin{pmatrix}1&1&1&1\\ \mathrm{i}A_{1}&-\mathrm{i}A_{1}&\mathrm{i}A_{2}&-\mathrm{i}A_{2}\\ e^{\widetilde{\lambda}_{1}\widetilde{d}}&e^{-\widetilde{\lambda}_{1}\widetilde{d}}&e^{\widetilde{\lambda}_{2}\widetilde{d}}&e^{-\widetilde{\lambda}_{2}\widetilde{d}}\\ \mathrm{i}A_{1}e^{\widetilde{\lambda}_{1}\widetilde{d}}&-\mathrm{i}A_{1}e^{-\widetilde{\lambda}_{1}\widetilde{d}}&\mathrm{i}A_{2}e^{\widetilde{\lambda}_{2}\widetilde{d}}&-\mathrm{i}A_{2}e^{-\widetilde{\lambda}_{2}\widetilde{d}}\end{pmatrix}.$

(22)

The existence of non-trivial solutions for the coefficients $\bm{a}$ that satisfy the boundary conditions Eq. (20) is granted provided $\det\mathbb{A}=0$. This equation can be expressed in a closed algebraic form as

	$\displaystyle\frac{\widetilde{\lambda}_{1}}{\widetilde{\lambda}_{2}}\frac{1+\epsilon-\widetilde{\lambda}_{2}^{2}/\mathcal{R}}{1+\epsilon-\widetilde{\lambda}_{1}^{2}/\mathcal{R}}v(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})=u(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})$
	$\displaystyle-\sqrt{u^{2}(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})-v^{2}(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})},$		(23)

with the definitions $u(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})\equiv\cosh(\widetilde{d}~{}\widetilde{\lambda}_{1})\cosh(\widetilde{d}~{}\widetilde{\lambda}_{2})-1$ and $v(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2})\equiv\sinh(\widetilde{d}~{}\widetilde{\lambda}_{1})\sinh(\widetilde{d}~{}\widetilde{\lambda}_{2})$, respectively.

Furthermore, from the first three rows of Eq. (21) we obtain

	$\displaystyle\bar{a}_{-}^{(1)}$	$\displaystyle\equiv$	$\displaystyle a_{-}^{(1)}/a_{+}^{(1)}=\frac{\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{2}}-1\right)e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}A_{1}+\left(1+e^{2\widetilde{d}~{}\widetilde{\lambda}_{2}}-2e^{\widetilde{d}(\widetilde{\lambda}_{1}+\widetilde{\lambda}_{2})}\right)e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}A_{2}}{e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{2}}-1\right)A_{1}-\left(e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}-2e^{\widetilde{d}~{}\widetilde{\lambda}_{2}}+e^{\widetilde{d}(\widetilde{\lambda}_{1}+2\widetilde{\lambda}_{2})}\right)A_{2}}$
	$\displaystyle\bar{a}_{+}^{(2)}$	$\displaystyle\equiv$	$\displaystyle a_{+}^{(2)}/a_{+}^{(1)}=\frac{\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{1}}-1\right)e^{\widetilde{d}~{}\widetilde{\lambda}_{2}}A_{2}-\left(e^{\widetilde{d}~{}\widetilde{\lambda}_{2}}-2e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}+e^{\widetilde{d}(2\widetilde{\lambda}_{1}+\widetilde{\lambda}_{2})}\right)A_{1}}{e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{2}}-1\right)A_{1}-\left(e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}-2e^{\widetilde{d}~{}\widetilde{\lambda}_{2}}+e^{\widetilde{d}(\widetilde{\lambda}_{1}+2\widetilde{\lambda}_{2})}\right)A_{2}}$
	$\displaystyle\bar{a}_{-}^{(2)}$	$\displaystyle\equiv$	$\displaystyle a_{-}^{(2)}/a_{+}^{(1)}=\frac{\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{1}}-1\right)e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}A_{2}+\left(1+e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}-2e^{\widetilde{d}(\widetilde{\lambda}_{1}+\widetilde{\lambda}_{2})}\right)e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}A_{1}}{e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}\left(e^{2\widetilde{d}~{}\widetilde{\lambda}_{2}}-1\right)A_{1}-\left(e^{\widetilde{d}~{}\widetilde{\lambda}_{1}}-2e^{\widetilde{d}~{}\widetilde{\lambda}_{2}}+e^{\widetilde{d}(\widetilde{\lambda}_{1}+2\widetilde{\lambda}_{2})}\right)A_{2}}.$		(24)

Note that the constants in Eq. (19) are the same for both spin projections, so we can finally write

	$\displaystyle\Psi_{\text{surface}}^{(\uparrow)}(\widetilde{z},\epsilon)$	$\displaystyle=$	$\displaystyle C_{\uparrow}\sum_{\alpha=1,2}\sum_{\beta=\pm}\bar{a}_{\beta}^{(\alpha)}\varphi_{\beta}^{(\alpha,\uparrow)}(\widetilde{z})$
	$\displaystyle\Psi_{\text{surface}}^{(\downarrow)}(\widetilde{z}.\epsilon)$	$\displaystyle=$	$\displaystyle C_{\downarrow}\sum_{\alpha=1,2}\sum_{\beta=\pm}\bar{a}_{\beta}^{(\alpha)}\varphi_{\beta}^{(\alpha,\downarrow)}(\widetilde{z}).$		(25)

Here, we defined $\bar{a}_{+}^{(1)}=1$, and the coefficients $C_{\uparrow}$, $C_{\downarrow}$ are determined by the normalization conditions as follows

	$\displaystyle|C_{\uparrow}|$	$\displaystyle=$	$\displaystyle\left[\sum_{\alpha=1,2}\sum_{\alpha^{\prime}=\pm}\sum_{\beta=1,2}\sum_{\beta^{\prime}=1,2}\bar{a}_{\beta}^{(\alpha)}\bar{a}_{\beta^{\prime}}^{(\alpha^{\prime})}\langle\alpha\uparrow;\beta|\alpha^{\prime}\uparrow;\beta^{\prime}\rangle\right]^{-1/2}$
	$\displaystyle|C_{\downarrow}|$	$\displaystyle=$	$\displaystyle\left[\sum_{\alpha=1,2}\sum_{\alpha^{\prime}=\pm}\sum_{\beta=1,2}\sum_{\beta^{\prime}=1,2}\bar{a}_{\beta}^{(\alpha)}\bar{a}_{\beta^{\prime}}^{(\alpha^{\prime})}\langle\alpha\downarrow;\beta|\alpha^{\prime}\downarrow;\beta^{\prime}\rangle\right]^{-1/2},$		(26)

where the analytical expression for the inner product coefficients are presented in Eq. (18). Moreover, from Eq. (17) it is clear that surface states with opposite spin components are mutually orthogonal

			$\displaystyle\langle\uparrow;\text{surface},\epsilon|\downarrow;\text{surface},\epsilon\rangle$		(27)
		$\displaystyle=$	$\displaystyle\int_{0}^{\tilde{d}}d\tilde{z}\left[\Psi_{\text{surface}}^{(\uparrow)}(\tilde{z},\epsilon)\right]^{\dagger}\Psi_{\text{surface}}^{(\downarrow)}(\tilde{z},\epsilon)=0.$

The system of equations Eq. (14) and Eq. (23) must be solved in favor of $\epsilon$. By inspection, it is evident that a solution is found when

$\displaystyle\widetilde{d}~{}\widetilde{\lambda}_{2}=0~{}~{}\text{and}~{}~{}\epsilon=-1$

(28)

so that, this condition implies, after Eq. (14)

$\displaystyle\widetilde{\lambda}_{1}=\mathcal{R}\left(1+\frac{2s}{\mathcal{R}}\right).$

(29)

On the other hand, nontrivial solutions corresponding to true surface states can exist when $\epsilon\simeq 0$, and those can only be found numerically as a function of the system thickness $d$.

An important property of these surface states, determined from the numerical solutions of the algebraic system of equations detailed in this section, is the well-defined parity of their components with respect to the center of the slab at $z=d/2$, i.e.

$$\Psi^{(\text{surface},\sigma)}_{j}(d/2-\zeta,\epsilon)=(-1)^{P}\Psi^{(\text{surface},\sigma)}_{j}(d/2+\zeta,\epsilon)$$

(30)

where $j$ labels the component ($j=1,\ldots,4$), and $\zeta\in[-d/2,d/2]$ represents the distance from the center of the slab. Here, the parity index $P$ is odd ($P=1$) or even ($P=2$), depending on the spin polarization of the surface state and the magnitude of $\mathcal{R}$. As shown in Fig. 2, for the $\sigma=\uparrow$ spin polarization, the only non-vanishing components are $j=1,3$. When $\mathcal{R}=0.1$ and $\mathcal{R}=1$, the component $j=1$ has even parity $P=2$, while the $j=3$ component possesses odd parity $P=1$. In contrast, we observe that the parities are reversed when $\mathcal{R}=0.5$.

A similar behavior is observed in Fig. 3 for the $\sigma=\downarrow$ spin polarization, when in this case it is the $j=2$ and $j=4$ components which are non-zero and alternate their well defined but opposite parities as a function of $\mathcal{R}$.

Nevertheless, despite both components always have opposite parities, the corresponding probability density $\rho_{{\rm{surface}}}^{(\sigma)}(z)=\Psi_{{\rm{surface}}}^{\sigma\dagger}(z)\Psi_{{\rm{surface}}}^{\sigma}(z)$ for either spin polarization $\sigma$ remains symmetric with respect to the center of the slab

$$\rho_{{\rm{surface}}}^{(\sigma)}(d/2-\zeta)=\rho_{{\rm{surface}}}^{(\sigma)}(d/2+\zeta).$$

(31)

The bulk states are obtained from the same system of equations Eq. (LABEL:eq_Dirac_system), but with the $z$-component of momentum $\tilde{k}_{z}$ a real number. From Eq. (11), we have that the general solution is given by the linear combination of spinors

	$\displaystyle\Phi^{\pm}_{\text{bulk}}(z)$	$\displaystyle=$	$\displaystyle\begin{pmatrix}\eta_{\text{bulk}}^{\pm}\\ \\ \chi_{\text{bulk}}^{\pm}\end{pmatrix}$		(32)
		$\displaystyle=$	$\displaystyle ue^{\pm\mathrm{i}\tilde{k}_{z}\widetilde{z}}\begin{pmatrix}1\\ 0\\ B_{z}\\ 0\end{pmatrix}+ve^{\pm\mathrm{i}\tilde{k}_{z}\widetilde{z}}\begin{pmatrix}0\\ 1\\ 0\\ -B_{z}\end{pmatrix},$

with

$\displaystyle B_{z}=\frac{s\tilde{k}_{z}}{1+\frac{s}{\mathcal{R}}\tilde{k}_{z}^{2}+\epsilon_{\text{bulk}}}.$

(33)

Now, by applying the confinement boundary conditions

$\displaystyle\phi_{\text{bulk}}(\widetilde{z}=0)=\phi_{\text{bulk}}(\widetilde{z}=\tilde{d})=0,$

(34)

we conclude that the appropriate linear combination (modulo a normalization constant) is

$$\phi_{\text{bulk}}(\widetilde{z})=\Phi^{+}_{\text{bulk}}(\widetilde{z})-\Phi^{-}_{\text{bulk}}(\widetilde{z})\propto\sin(\tilde{k}_{z}\widetilde{z}),$$

(35)

provided

$\displaystyle\sin(\tilde{k}_{z}\tilde{d})=0,$

(36)

which leads to the quantization condition for $\tilde{k}_{z}$

$\displaystyle\widetilde{k}_{z}=\frac{n\pi}{\widetilde{d}},~{}\text{for}~{}n=\pm 1,\pm 2,\ldots$

(37)

Therefore, from Eq. (13) the quantized energy eigenvalues due to spatial confinement in the region $0\leq z\leq d$ are given by

$\displaystyle\epsilon_{\text{bulk}}^{2}=\frac{n^{2}\pi^{2}}{\widetilde{d}^{2}}+\left(1+\frac{s}{\mathcal{R}}\frac{n^{2}\pi^{2}}{\widetilde{d}^{2}}\right)^{2}.$

(38)

The corresponding normalized bulk eigenstates for each spin projection are

	$\displaystyle\phi_{\text{bulk}}^{(n,\uparrow)}(\widetilde{z})$	$\displaystyle=$	$\displaystyle\sqrt{\frac{2}{(1+B_{n}^{2})\tilde{d}}}\sin\left(\frac{n\pi\widetilde{z}}{\widetilde{d}}\right)\left(1,0,B_{n},0\right)^{T},$
	$\displaystyle\phi_{\text{bulk}}^{(n,\downarrow)}(\widetilde{z})$	$\displaystyle=$	$\displaystyle\sqrt{\frac{2}{(1+B_{n}^{2})\tilde{d}}}\sin\left(\frac{n\pi\widetilde{z}}{\widetilde{d}}\right)\left(0,1,0,-B_{n}\right)^{T},$

with $\widetilde{d}\equiv d/R_{0}$, and the coefficients $B_{n}$ defined as

$$B_{n}=B_{z}\left(\tilde{k}_{z}=\frac{n\pi}{\widetilde{d}}\right)=\frac{sn\pi/\widetilde{d}}{1+\frac{s}{\mathcal{R}}\left(\frac{n\pi}{\widetilde{d}}\right)^{2}+\epsilon_{\text{bulk}}}.$$

(40)

The normalization constants are defined such that (for $\sigma=\uparrow,\downarrow$)

$\displaystyle\int_{0}^{\tilde{d}}d\widetilde{z}\left[\phi_{\text{bulk}}^{n,\sigma}(\widetilde{z})\right]^{\dagger}\phi_{\text{bulk}}^{n,\sigma}(\widetilde{z})=1.$

(41)

A noteworthy property of these spin projection eigenstates is their even/odd parity with respect to the center of the slab $z=d/2$, since for $\zeta\in[-d/2,d/2]$ the distance from the center of the slab, we have

$$\sin\left(\frac{n\pi}{d}\left(d/2-\zeta\right)\right)=(-1)^{n+1}\sin\left(\frac{n\pi}{d}\left(d/2+\zeta\right)\right).$$

(42)

Therefore, the same even/odd parity is inherited by the spinor functions

$$\phi_{\text{bulk}}^{n,\sigma}\left(d/2-\zeta\right)=(-1)^{n+1}\phi_{\text{bulk}}^{n,\sigma}\left(d/2+\zeta\right).$$

(43)

As discussed in Eq. (30), depending on the magnitude of $\mathcal{R}$, the spinor components of the surface states possess well defined but opposite parities $P$ with respect to a mirror plane located at the center of the slab $z=d/2$. On the other hand, the bulk states exhibit a global parity determined by the number of nodes within the finite domain $z\in[0,d]$ imposed by the boundaries of the slab, which depends on the quantum number $n\geq 1$, as stated in Eq. (43).

Let us now consider the existence of a hybrid surface-bulk state described by the linear combination

$\displaystyle\Psi_{n,\epsilon}^{\left(\sigma,\sigma^{\prime}\right)}(\widetilde{z})=\alpha\Psi^{(\sigma)}_{\text{surface}}(\widetilde{z},\epsilon)+\beta\phi_{\text{bulk}}^{(n,\sigma^{\prime})}(\widetilde{z}).$

(44)

The probability distribution corresponding to this hybrid state is given by the expression

$\displaystyle\rho^{(\sigma,\sigma^{\prime})}_{n,\epsilon}(z)=\rho_{\text{surface}}^{(\sigma)}(z)+\rho_{\text{bulk}}^{(\sigma^{\prime})}(z)+\rho_{\text{mix}}^{(\sigma,\sigma^{\prime})}(z),$

(45)

where

$\displaystyle\rho_{\text{mix}}^{(\sigma,\sigma^{\prime})}=2\text{Re}\left\{\alpha^{*}\beta\braket{\text{surface},\sigma}{\text{bulk},\sigma^{\prime}}\right\}$

(46)

represents the contribution due to the quantum interference effect between surface and bulk in the hybrid state in Eq. (44). As a trivial consequence of Eq. (30) and Eq. (43), the probability density contributions $\rho_{\text{surface}}^{(\sigma)}(z)=|\alpha|^{2}|\Psi_{\text{surface}}^{(\sigma)}(z,\epsilon)|^{2}$ and $\rho_{\text{bulk}}^{(\sigma)}(z)=|\beta|^{2}|\phi_{\text{bulk}}^{(n,\sigma)}(z)|^{2}$ arising independently from surface and bulk, respectively, are symmetric with respect to the center of the slab, regardless of their spin or parity. This feature is clearly illustrated in the central row of the panel Fig. 4, by the solid (bulk) and dashed (surface) lines, respectively.

Clearly, for opposite spins $\sigma\neq\sigma^{\prime}$ the components of the surface and bulk spinor states are mutually orthogonal such that $\rho_{\text{mix}}^{(\sigma\neq\sigma^{\prime})}=0$, and hence the probability density is indeed trivially symmetric with respect to a mirror plane at the center of the slab $z=d/2$, i.e.

$$\rho_{n,\epsilon}^{(\sigma\neq\sigma^{\prime})}(d/2-\zeta)=\rho_{n,\epsilon}^{(\sigma\neq\sigma^{\prime})}(d/2+\zeta),$$

(47)

for $\zeta\in[-d/2,d/2]$ the distance from the center. This is clearly shown in the first row of the panel Fig. 4, where the dotted line represents the total probability density $\rho^{(\uparrow,\downarrow)}$ for opposite spin components in the hybrid state when $|\alpha|=|\beta|$ in the superposition Eq. (44).

In contrast, a very non-trivial behavior is observed when surface and bulk in the hybrid state Eq. (44) possess parallel spins ($\sigma=\sigma^{\prime}$). In this last case, the mixture term does not vanish, as illustrated in the bottom row of the panel Fig. 4, for the choice $|\alpha|=|\beta|$ in the combination of the hybrid state in Eq. (44). Interestingly, given that the surface spinor components possess alternate parities (see Eq. (30)), while the bulk spinor possesses a well defined global parity (see Eq. (43)), the mixture term $\rho_{{\rm{mix}}}^{(\sigma,\sigma)}$ does not exhibit, in general, a well defined parity for arbitrary values of $n$ and $\mathcal{R}$. However, as can be appreciated in Fig. 2 for $\sigma=\uparrow$ and Fig. 3 for $\sigma=\downarrow$, respectively, when $\mathcal{R}\lesssim 1$ we have that

$\displaystyle\left|\text{Re}\Psi^{(\text{surface},\sigma)}_{j}\right|\ll\left|\text{Re}\Psi^{(\text{surface},\sigma)}_{j+2}\right|,$

(48)

and in such case the mix term acquires an approximate symmetry imposed by the parity of the bulk state. This property of the mixture density term is appreciated in Fig. 5, for $\mathcal{R}=0.1$ and different values of $n$, the later determining the parity of the bulk state according to Eq. (43). In consequence, the mirror symmetry of the total probability density for the hybrid surface-bulk state defined by Eq. (45) may be broken for parallel spin components, depending on the combination of values of $\mathcal{R}$ and $n$. These occurrences are material-specific because $\mathcal{R}$ depends on the Compton decay length. Consequently, the subtleties of symmetry breaking remain distinctive for each material due to their unique characteristics. A graphical illustration of this effect is observed in the first row of the panel Fig. 4, where the solid line represents the total probability density $\rho^{(\uparrow,\uparrow)}$ for parallel spin components in the hybrid state when $|\alpha|=|\beta|$ in the superposition Eq. (44).

This rather counter-intuitive result, given the symmetric boundary conditions in the slab geometry, is clearly a pure quantum mechanical effect due to interference contained in the mixture term.

Figure 6 highlights the effect of the coefficients $\alpha$ and $\beta$ in the linear combination described by Eq. (44). The general considerations regarding inversion symmetry remain unaffected. As expected, larger values of $|\alpha|>|\beta|$ imply an increasingly important contribution from the surface state in the hybrid combination Eq. (44), and correspondingly the weight of the probability distribution is more concentrated towards the edges. Conversely, when $|\beta|>|\alpha|$ the probability density of the hybrid state displays more weight towards the bulk. A similar effect is observed as a function of $n$, since for $n\gg 1$, the hybrid state displays more statistical weight in the bulk.

By means of a standard continuum effective model for a TI in a slab geometry, we explicitly obtained the surface and confined bulk states, in order to investigate the potential effects of their hybridization in the probability density, an important physical parameter that determines the electronic transport properties of such materials near the Fermi level. We assumed symmetric boundary conditions that consistently lead to probability density distributions for the bulk and surface states that possess a mirror symmetry with respect to a plane at the center of the slab. Surprisingly, our analytical and numerical results reveal that, when the spin components of surface and bulk in the hybrid state are parallel to each other, the symmetry in the probability density is spontaneously broken under certain conditions, depending on the combination of $n$ (the quantum number in the bulk solution) and the material-dependent parameter $\mathcal{R}$. We can explain this effect in simple terms, due to the spatially asymmetric interference pattern arising between surface and bulk, as captured in the $\rho_{{\text{mix}}}^{(\sigma,\sigma)}$ contribution for parallel spin components.

In this study, we are not delving into the detailed microscopic mechanism leading to the surface-bulk hybridization, which has been argued to be triggered, for instance, by disorder effects Saha and Garate (2014); Black-Schaffer and Balatsky (2012); Mandal et al. (2023); Yang et al. (2020); Hikami et al. (1980), or analyzed by means of a Fano model with a phenomenological coupling constant Hsu et al. (2014) in previous works. In contrast, we are here focusing on its possible consequences on the resulting probability density distribution. We found that under certain configurations, when the surface and bulk spin polarizations in the hybrid state are parallel, it may display an asymmetric probability density distribution, a theoretical result not previously reported in the literature Hsu et al. (2014); Saha and Garate (2014); Black-Schaffer and Balatsky (2012); Mandal et al. (2023); Yang et al. (2020); Hikami et al. (1980), even though the importance of spin polarization for the emergence of the surface state was already recognized by first-principles calculations in slab geometries Reid et al. (2020).

An experimental scenario where this interesting effect could be observed involves the presence of a constant, weak magnetic field $\mathbf{B}=\hat{z}B$ acting upon the topological insulator slab, as depicted in Fig. 7. In this scenario, the external field acts as a perturbation that does not alter the system’s ground state but breaks its spin degeneracy via the Zeeman coupling. This coupling, involving a magnetic moment $\mu$, results in a linear shift in the corresponding energy eigenvalues described by Eq. (9) as $\epsilon\to\epsilon-s\mu B/2$, where $s=-1$ denotes the $(\uparrow,\downarrow)$ state and $s=+1$ represents the $(\uparrow,\uparrow)$ state. Consequently, under this condition the system tends to favor a hybrid quantum state with parallel spin components for surface and bulk, where we anticipate the emergence of this novel phenomenon. Supporting experimental evidence for this conjecture is indeed presented in the literature Buchenau et al. (2017). For instance, a study on magnetotransport in Bi₂Se₃ nanoplate devices reveals that the proximity of a magnetic substrate to a surface channel disrupts the symmetry between the top and bottom electronic densities, resulting in the decoupling of the quantum Hall effects on each surface Buchenau et al. (2017).

JDCY and EM acknowledge financial support from ANID PIA ANILLO ACT/192023. EM also acknowledges ANID FONDECYT No 1230440.

The authors declare no conflict of interest.

All data generated or analysed during this study are included in this published article. Any additional information is available from the corresponding author on reasonable request.

Topological insulators, surface states, bulk states, surface-bulk hybridization, symmetry breaking, electronic density.