\section{Introduction}

Over the last decade, the effective model Hamiltonian for topological surface states developed in Refs.~[\onlinecite{Zhang_NATPHYS_2009}] and~[\onlinecite{Liu_PRB_2010}] is commonly accepted as a tool to include in a simple manner their remarkable features: the linear energy-momentum dispersion and helical in-plane spin structure. This model has been applied to a variety of topologically non-trivial materials in the spirit of the classical Rashba model, which, since the seminal paper by LaShell et al.~[\onlinecite{LaShell_PRL_1996}], has been used to fit the two-dimensional (2D) spin-orbit-split states at trivial surfaces. The spin-orbit splitting $k_{\pm}$-linear term is the same for trivial and for topological surface states, and it yields an orthogonal spin-momentum locking commonly considered a hallmark of a strong spin-orbit interaction (SOI).

The simplified model of Refs.~[\onlinecite{Zhang_NATPHYS_2009}] and [\onlinecite{Liu_PRB_2010}] needs to be extended in order to describe the non-orthogonality between spin and momentum in realistic systems, as, e.g., observed in photoemission from Bi$_2$Se$_3$~[\onlinecite{Wang_PRL_2011}]. While it naturally arises in \textit{ab initio} calculations, in $\mathbf{k}\cdot\mathbf{p}$ theory, in order to yield a deviation from orthogonality, an effective Hamiltonian must include higher-order in $\mathbf{k}$ spin-orbit terms. In Ref.~[\onlinecite{Basak_PRB_2011}], for structures with the $C_{3v}$ crystal symmetry and time-reversal symmetry it was suggested to include a $k_{\pm}^5$ term to allow for the non-orthogonality, and a minimal fifth-order $\mathbf{k}\cdot\mathbf{p}$ model was applied to Bi$_2$Te$_3$. Following Ref.~[\onlinecite{Basak_PRB_2011}], in Ref.~[\onlinecite{Hopfner_PRL_2012}] this model was also used to analyze the spin structure of the Au/Ge(111) surface state, which is rather far from being Rashba-like. Since the fifth-order Hamiltonian was constructed based on symmetry arguments rather than derived directly from \textit{ab initio} spinor wave functions, its parameters were found by fitting to the \textit{ab initio} band structure, which is an approximate procedure sensitive to the choice of the energy interval of interest and to the order of the $\mathbf{k}\cdot\mathbf{p}$ expansion: with each successive order the complexity and ambiguity grow, and the parameters become increasingly less physically meaningful.

Here, we study the angle between the spin and momentum in the surface states of Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se within our \textit{ab initio} relativistic $\mathbf{k}\cdot\mathbf{p}$ approach introduced in Refs.~[\onlinecite{Nechaev_PRBR_2016, Nechaev_PRB_2018, Nechaev_PRB_2019}]. This approach has been successfully applied to different materials~[\onlinecite{Nechaev_SciRep_2017, Nechaev_PRB_2019, Susanne2019, Usachov_arXiv_2020}] and established as a reliable theoretical tool for deriving few-band $\mathbf{k}\cdot\mathbf{p}$ Hamiltonians capable of comprehensive description of the surface spin structure. We take Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se as vivid examples of the topological insulators (TIs) with a rather wide absolute bulk band gap bridged by the partly occupied topological surface state and a local projected gap well above the Fermi level hosting the so-called ``second topological surface state''~[\onlinecite{Niesner_PRB_2012, Sobota_PRL_2013, Niesner_JESRP_2014, Datzer_PRB_2017, Aguilera_PRB_2019}]. The wide gap is favorable for minimizing the effect of the proximity of bulk states on the surface-state spin structure. The presence of the second surface state makes it possible to derive a \textit{two-band seventh-order} Hamiltonian by applying the L\"{o}wdin partitioning to a \textit{four-band third-order} Hamiltonian generated for the two surface states within our \textit{ab initio} approach.

For Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se, we derive the seventh-order Hamiltonian allowing for the non-orthogonal locking between momentum and spin. We show that there is an oscillation of the spin around the momentum-perpendicular direction with a $\pi/3$ periodicity as a function of the polar angle of $\mathbf{k}$ due to the presence of $\mathbf{k}$-dependent effective spin-orbit magnetic fields of the six-fold symmetry. The effective Hamiltonian facilitates the inclusion of the non-orthogonality in the description of spin-related properties of the TI surfaces and their interpretation within $\mathbf{k}\cdot\mathbf{p}$ theory. Thus, our study can also be considered as an \textit{ab initio} substantiation of the fifth-order Hamiltonian proposed in Ref.~[\onlinecite{Wang_PRL_2011}], based on an unambiguous algorithm for its parameters.

\section{Computational details}

The \textit{ab initio} band structure is obtained with the extended linear augmented plane waves method~[\onlinecite{Krasovskii_PRB_1997}] using the full potential scheme of Ref.~[\onlinecite{Krasovskii_PRB_1999}] within the local density approximation (LDA). The spin-orbit interaction was treated as a second variation~[\onlinecite{Koelling_1977}]. The surfaces of the TIs are simulated by bulk-truncated centrosymmetric six-QL (quintuple layer) slabs of space group $P\bar{3}m1$ (no.~164). The experimental crystal lattice parameters were taken from Ref.~[\onlinecite{Wyckoff_RWG}]. In the case of Bi$_2$Te$_2$Se, the experimental atomic positions of Ref.~[\onlinecite{Wyckoff_RWG}] were used, while for Bi$_2$Se$_3$ we took the LDA relaxed atomic positions of Ref.~[\onlinecite{Nechaev_PRB_2013_BISE}].

\section{Spin-momentum locking angle}

Figure~\ref{fig1} shows the calculated LDA band structure of the Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se surfaces along $\bar{\Gamma}$-$\bar{K}$. Two topological surface states are clearly identified in the spectra of both TIs. The two states are numbered $n=1$ and 2 in order of increasing energy. The spin-resolved constant energy contours (CECs) at $\sim0.2$~eV above the Dirac points of the low-energy surface states ($n=1$) are shown in Fig.~\ref{fig2} together with the respective angles of deviation from the orthogonal spin-momentum coupling $\delta$. As a function of the polar angle $\varphi_{{\mathbf k}}$ of the momentum $\mathbf{k}$, the deviation angle demonstrates an oscillating behavior with a $\pi/3$ periodicity and an amplitude close to 1.5$^{\circ}$ for Bi$_2$Se$_3$, which is in good agreement with the experiment of Ref.~[\onlinecite{Wang_PRL_2011}], and about 3.0$^{\circ}$ for Bi$_2$Te$_2$Se.

\begin{figure}[tbp]
\centering
\includegraphics[width=\columnwidth]{Fig_1.png}
\caption{Band structure of the surface of Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se along $\bar{\Gamma}$-$\bar{K}$ by the full LDA Hamiltonian $H^{\rtm{LDA}}_{{\mathbf k}}$ (black lines), the $8\times8$ Hamiltonian~(\ref{HamTwoSurfaces}) (light green lines), the $4\times4$ Hamiltonian~(\ref{TWOSSHam}) (green lines), and the $2\times2$ Hamiltonian~(\ref{UPTO7}) with two different values of $\widetilde{\gamma}^{(7)}$ (dark green and orange lines, see text). The upper (lower) border of the horizontal stripes corresponds to the energy at which the model (\textit{ab initio}) CECs shown in Fig.~\ref{fig2} are calculated. The $\mathbf{k}$-points on these CECs along $\bar{\Gamma}$-$\bar{K}$ are marked by vertical blue lines.}
\label{fig1}
\end{figure}

We start with a ${\mathbf k}\cdot{\mathbf p}$ model comprising both Dirac surface states $n=1$ and 2. Being eigenfunctions of a centrosymmetric slab Hamiltonian at $\bar{\Gamma}$, these states form four Kramers-degenerate pairs with spinor wave functions $\Psi_{m\mu}$, which we group into two twin pairs with two members, $m = 2n-1$ and $2n$. Here, $\mu=\uparrow$ or $\downarrow$ indicates the sign $+$ or $-$ of the expectation value $\langle\Psi_{m\mu}| \widehat{J}_{z} |\Psi_{m'\mu'}\rangle_{\tau}=\langle J_z\rangle_{m\mu}\delta_{mm'}\delta_{\mu\mu'}$ of the $z$ projection of the total angular momentum $\widehat{\mathbf{J}}$ at the (symmetry equivalent) atomic sites of type $\tau$, which has the largest weight $\langle\Psi_{m\mu}|\Psi_{m\mu}\rangle_{\tau}$, see Ref.~[\onlinecite{Nechaev_PRBR_2016}]. The integration is over the muffin-tin spheres of this type, and the positive value is $\langle J_z\rangle_{m\uparrow} = -\langle J_z\rangle_{m\downarrow}$. With this basis set, we first microscopically derive an eight-band $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian $H_{\rtm{\mathbf{kp}}}$ from an \textit{ab initio} relativistic ${\mathbf k}\cdot{\mathbf p}$ perturbation expansion around the $\Gamma$ point. The expansion is carried out up to the third order in \textbf{k} by applying  the L\"{o}wdin partitioning~[\onlinecite{Leowdin_JCP_1951, Schrieffer_PR_1966, Winkler_KP}] to the original Hilbert space of the $\Gamma$-projected LDA Hamiltonian $H^{\rtm{LDA}}_{{\mathbf k}}$, see Appendix~\ref{A_expansion}.

Because $\Psi_{m\mu}$ are slab eigenfunctions representing the surface states, each twin-pair is characterized by two doubly degenerate slab levels $E_{2n-1}$ and $E_{2n}$ separated by $\Delta_n = E_{2n}-E_{2n-1}$ of a few meV due to the bonding-antibonding interaction. Since the $\bar{\Gamma}$ point is a TRIM (time reversal invariant momentum), the spinors $\Psi_{m\mu}$ are also parity eigenfunctions, and the two pairs for a given $n=1$ or 2 have different parity. We now transfer to a new basis $|\Phi^{\pm}_{n\mu}\rangle =\frac{1}{\sqrt{2}} \left[|\Psi_{2n-1\mu}\rangle \pm |\Psi_{2n\mu}\rangle\right]$, where the new basis functions $|\Phi^{\pm}_{n\mu}\rangle$ are no longer parity eigenfunctions but are localized at one of the two surfaces of the 6QL slab, ``$+$'' or ``$-$''. In this surface-resolved basis, the original $8\times8$ Hamiltonian reads
\begin{equation}\label{HamTwoSurfaces}
H_{\rtm{\mathbf{kp}}}\longrightarrow H^{\rtm{Film}}_{\rtm{\mathbf{kp}}}=\left(
\begin{array}{cc}
H_{\rtm{Surf}}^{+} & H_{\rtm{int}}  \\
H^{\dag}_{\rtm{int}} & H_{\rtm{Surf}}^{-}
\end{array}
\right).
\end{equation}
In Fig.~\ref{fig1}, the bands obtained by diagonalizing this Hamiltonian are shown by light green lines for both TIs.

\begin{figure*}[tbp]
\centering
\includegraphics[width=\textwidth]{Fig_2.png}
\caption{Spin-resolved constant energy contours (upper panels) and the non-orthogonality $\delta$ as a function of the polar angle $\varphi_{\mathbf{k}}$ (lower panels) for Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se. In the upper panels, the colored areas highlight the deviation of the in-plane spin direction from the classical Rashba in-plane spin at a given $\mathbf{k}$-point in the contour: the border of the areas is given by the length of the relevant $\mathbf{k}$ plus a scaled $\mathbf{k}$-projection of the in-plane spin: $k+R\sin\delta$, $R$ being the scaling factor. The green areas and lines correspond to the LDA and ${\mathbf k}\cdot{\mathbf p}$ calculations, while the orange ones -- to the ${\mathbf k}\cdot{\mathbf p}$ calculations with the magnitude of $\gamma^{(5)}$ manually increased by a factor of 4.5 for Bi$_2$Se$_3$ and of 3.8 for Bi$_2$Te$_2$Se, see text. The contours are calculated at the energies marked in Fig.~\ref{fig1}.}
\label{fig2}
\end{figure*}

Further we neglect the coupling of the surfaces due to the overlap between the $+$ and $-$ new basis functions, $H_{\rtm{int}}\rightarrow0$, and in the following we will consider only the $-$ surface, so we omit the superscript $-$. In a compact form, the resulting $4\times4$ Hamiltonian, which is just the term $H_{\rtm{Surf}}^{-}$ of the Hamiltonian~(\ref{HamTwoSurfaces}), reads
\begin{equation}\label{TWOSSHam}
H_{\rtm{\mathbf{kp}}}^{4\times4}=\left(
\begin{array}{cc}
E_1+H_1+H^{\rtm{R}}_1      & H_0 +\widetilde{H} \\
H_0 +\widetilde{H}^{\dag}  & E_2+H_2+H^{\rtm{R}}_2
\end{array}
\right).
\end{equation}
Here, each term of the diagonal and non-diagonal blocks is a $2\times2$ matrix, whose implicit form directly follows from the \textit{ab initio} ${\mathbf k}\cdot{\mathbf p}$ expansion up to the third order in \textbf{k}: $E_n=\epsilon_n\rtm{\mathbb{I}}_{2\times2}$, with $\rtm{\mathbb{I}}_{2\times2}$ being the $2\times2$ identity matrix, $H_n=M_{n}k^2\rtm{\mathbb{I}}_{2\times2}$, $k=\sqrt{k_x^2+k_y^2}$, and
\begin{equation}\label{Rashba2}
H^{\rtm{R}}_n=\left(
\begin{array}{cc}
-iW_n (k_+^3-k_-^3)      & i\tilde{\alpha}_{n}k_- \\
-i\tilde{\alpha}_{n}k_+ & iW_n(k_+^3-k_-^3)
\end{array}
\right),
\end{equation}
where $\tilde{\alpha}_{n}=\alpha_n^{(1)}+\alpha_n^{(3)}k^2$ and $k_{\pm}=k_x\pm k_y$. The well-known $2\times2$ Rashba term $H^{\rtm{R}}_n$ is responsible for the out-of-plane and in-plane spin structure typical of hexagonal structures, see, e.g, Ref.~[\onlinecite{Nechaev_PRB_2019}] and references therein. The interaction between the states $n=1$ and 2 is realized through the term
\begin{equation}\label{HintTWOSS}
\widetilde{H}=\left(
\begin{array}{cc}
i\theta k_+^3+i\eta k_-^3      & i\tilde{\zeta}k_- +Dk_+^2 \\
-i\tilde{\zeta}k_+ -Dk_-^2 & i\theta k_-^3+i\eta k_+^3
\end{array}
\right)
\end{equation}
with $\tilde{\zeta}=\zeta^{(1)}+\zeta^{(3)}k^2$.

\begin{table}
\caption{\label{tab:table1} Parameters of the four-band Hamiltonians~(\ref{TWOSSHam}) (based on calculations for 6QL-layer slabs with the lattice parameter $a=7.8187$~a.u.  for Bi$_2$Se$_3$ and $a=8.0880$~a.u. for Bi$_2$Te$_2$Se). All parameters are in Rydberg atomic units except for $\epsilon_1$ and $\epsilon_2$ presented in eV.}
\begin{ruledtabular}
\begin{tabular}{ldd}
                                                            &       \mbox{Bi$_2$Se$_3$}     &  \mbox{Bi$_2$Te$_2$Se}  \\
\hline
$\epsilon_1$                                    &   -0.118                   &   -0.278   \\
$\epsilon_2$                                    &    1.429                    &   1.053  \\
$\alpha_1^{(1)}$                           &    0.174                     &   0.187  \\
$\alpha_2^{(1)}$                          &    -0.265                     &   -0.150  \\
$\alpha_1^{(3)}$                          &     28.30                    &  -29.03     \\
$\alpha_2^{(3)}$                          &    141.80                  &   -22.35     \\
$\theta$                                           &   8.33                       &   22.39  \\
$\eta$                                              &    -1.02                        &  -3.00 \\
$\zeta^{(1)}$                                 &    -0.048                     &    -0.078 \\
$\zeta^{(3)}$                                 &    -52.83                     &    4.51   \\
$D$                                                   &   -2.64                       &   -3.40 \\
$M_1$                                               &   7.97                        &    15.67 \\
$M_2$                                               &  -2.56                        &    -1.76 \\
$M_0$                                               &  -0.37                        &    0.61 \\
$W_1$                                               &  -4.71                        &  -17.32 \\
$W_2$                                               &   5.53                       &    11.88  \\
$s^{\shortparallel}_1$                                    &    0.70                       &   0.63  \\
$s^{\shortparallel}_2$                                    &    0.42                       &   0.39  \\
$ s^{z}_1$                                        &     0.40                       &   0.26 \\
$ s^{z}_2$                                        &   -0.16                       &   -0.21 \\
$\tilde{s}^{\shortparallel}$                          &  -0.21                        &  -0.10  \\
$\tilde{s}^{z}$                                &  -0.41                       &   -0.20 \\
\end{tabular}
\end{ruledtabular}
\end{table}

In the new basis, the spin matrix that yields the spin structure of the states under study is defined as
\begin{equation}\label{SpinTWOSS}
\rtm{\mathbf{S}}^{4\times4}_{\rtm{\mathbf{kp}}}=\left(
\begin{array}{cc}
\rtm{\mathbf{S}}_1 & \widetilde{\rtm{\mathbf{S}}}  \\
\widetilde{\rtm{\mathbf{S}}} & \rtm{\mathbf{S}}_2
\end{array}
\right)
\end{equation}
with $\rtm{\mathbf{S}}_n=(s^{\shortparallel}_n\bm{\sigma}_{\shortparallel}, s^{z}_n\sigma_z)$ and $\widetilde{\rtm{\mathbf{S}}}=(\tilde{ s}^{\shortparallel}\bm{\sigma}_{\shortparallel}, \tilde{s}^{z}\sigma_z)$, where  $\bm{\sigma}_{\shortparallel} = (\sigma_x, \sigma_y)$ and $\sigma_x$, $\sigma_y$, and $\sigma_z$ are the Pauli matrices. The elements of the spin matrix $\left[\rtm{\mathbf{S}}^{4\times4}_{\rtm{\mathbf{kp}}} \right]^{n\mu}_{l\nu} = \langle\Phi_{n\mu}|\bm{\sigma}|\Phi_{l\nu}\rangle$ enter into the expression for the spin expectation value
\begin{equation}\label{modelRealSpin}
\langle \mathbf{S}_{\mathbf{k}\lambda}\rangle = \frac{1}{2} \langle \widetilde{\Phi}^{\lambda}_{\mathbf{k}}|\bm{\sigma}|\widetilde{\Phi}^{\lambda}_{\mathbf{k}}\rangle
= \frac{1}{2}\sum\limits_{n\mu l\nu} C_{{\mathbf{k}}n\mu}^{\lambda\ast}C_{{\mathbf{k}}l\nu}^{\lambda}\left[\rtm{\mathbf{S}}^{4\times4}_{\rtm{\mathbf{kp}}} \right]^{n\mu}_{l\nu}
\end{equation}
in the state $|\widetilde{\Phi}^{\lambda}_{\mathbf{k}}\rangle =\sum\limits_{n\mu}C_{\mathbf{k}n\mu}^{\lambda} |\Phi_{n\mu}\rangle$ of the reduced Hilbert space of the Hamiltonian $H^{4\times4}_{\rtm{\mathbf{kp}}}$. The four-dimensional vectors $\mathbf{C}^{\lambda}_{\mathbf{k}}$ diagonalize this Hamiltonian  $H^{4\times4}_{\rtm{\mathbf{kp}}} \mathbf{C}^{\lambda}_{\mathbf{k}} = E^{\lambda}_{\mathbf{k}} \mathbf{C}^{\lambda}_{\mathbf{k}}$. The parameters in Eqs.~(\ref{TWOSSHam}) and (\ref{SpinTWOSS}) are listed in Table~\ref{tab:table1}. The bands obtained with these parameters are shown in Fig.~\ref{fig1} by green lines.

Next, we \textit{analytically} transform the Hamiltonian~(\ref{TWOSSHam}) by means of the L\"{o}wdin partitioning, retaining terms up to seventh-order in $\mathbf{k}$ for the block $E_1+H_1+H^{\rtm{R}}_1$ of this Hamiltonian. As a result, we arrive at the $2\times2$ Hamiltonian that describes the low-energy Dirac surface state:
\begin{widetext}
\begin{equation}\label{UPTO7}
H_{\rtm{\mathbf{kp}}}^{2\times2}=\left(
\begin{array}{cc}
\epsilon_1+\mathcal{M}k^2-i\mathcal{W}(k_+^3-k_-^3)+\mathcal{N}(k_+^6+k_-^6)  & i\tilde{\alpha}k_--i\tilde{\gamma}k_+^5+i\xi k_-^7 \\
-i\tilde{\alpha}k_++i\tilde{\gamma}k_-^5-i\xi k_+^7  & \epsilon_1+\mathcal{M}k^2+i\mathcal{W}(k_+^3-k_-^3)+\mathcal{N}(k_+^6+k_-^6)
\end{array}
\right),
\end{equation}
\end{widetext}
where $\mathcal{M}=\sum_{m=0}^{2}M^{(2m)}k^{2m}$, $\mathcal{W}=\sum_{m=0}^{2}W^{(2m)}k^{2m}$, $\tilde{\alpha}=\alpha_1^{(1)}+\sum_{m=1}^{3}\alpha^{(2m+1)}k^{2m}$, and $\tilde{\gamma}=\gamma^{(5)}+\gamma^{(7)}k^2$, see Appendix~\ref{A_params}. All the parameters are listed in Table~\ref{tab:table2}. With these parameters, the diagonalization of the Hamiltonian~(\ref{UPTO7}) yields the bands shown by dark green lines in Fig.~\ref{fig1}.

Up to fifth order in \textbf{k}, the Hamiltonian~(\ref{UPTO7}) of our \textit{ab initio} ${\mathbf k}\cdot{\mathbf p}$ theory is in accord with the form of  the two-band Hamiltonian constructed in Ref.~[\onlinecite{Basak_PRB_2011}] for Bi$_2$Te$_3$ considering the $C_{3v}$ crystal symmetry and time-reversal symmetry. In Ref.~[\onlinecite{Hopfner_PRL_2012}], the Hamiltonian of Ref.~[\onlinecite{Basak_PRB_2011}] was modified by adding a sixth-order term $k_+^6+k_-^6$ in order to reproduce the hexagonal warping of the Au/Ge(111) surface state not related to the spin-orbit effect. Obviously, this term is naturally present in our theory. Note that in Refs.~[\onlinecite{Basak_PRB_2011}] and [\onlinecite{Hopfner_PRL_2012}] the values of  the parameters were found by fitting the model Hamiltonian to \textit{ab initio} results, and, for example, in the case of the surface state of Bi$_2$Te$_3$~[\onlinecite{Basak_PRB_2011}] the values for the lower-order terms differ strongly from those obtained in Refs.~[\onlinecite{Liu_PRB_2010}] and [\onlinecite{Nechaev_PRBR_2016}], which are currently commonly accepted. In contrast to a fitting method, in our $\mathbf{k}\cdot\mathbf{p}$ theory the shape and the value of a given order term are independent on whether or not we include higher-order terms (and it does not affect the lower-order terms), since our $\mathbf{k}\cdot\mathbf{p}$ expansion uniquely follows from the basis set---the eigenfunctions of the original {\it ab initio} Hamiltonian.

\begin{table}
\caption{\label{tab:table2} Parameters of the two-band Hamiltonians~(\ref{UPTO7}) in Rydberg atomic units. The parameters $\epsilon_1$, $\alpha_1^{(1)}$, $s^{\shortparallel}_1$, $ s^{z}_1$, and the lattice parameters $a$ are listed in Table~\ref{tab:table1}. }
\begin{ruledtabular}
\begin{tabular}{ldd}
                                                            &       \mbox{Bi$_2$Se$_3$}     &  \mbox{Bi$_2$Te$_2$Se}  \\
\hline
$\alpha^{(3)}$                               & 27.91                      &    -28.27    \\
$\alpha^{(5)}$                               &  -575.48                    &  115.57     \\
$\alpha^{(7)}$                               &  -92731.72               &   63344.71  \\
$\gamma^{(5)}$                            & 529.40                   &   1735.90    \\
$\gamma^{(7)}$                            &  -39773.06              &   319791.93  \\
$\xi$                                                  &  51.10                    &    589.48  \\
$\mathcal{N}$                               &  -157.85                   &     -789.83\\
$\mathcal{M}^{(0)}$                   &   7.95                         &     15.61 \\
$\mathcal{M}^{(2)}$                   &   -110.14                  &    -122.40 \\
$\mathcal{M}^{(4)}$                   & -36925.53                &    -27778.63 \\
$\mathcal{W}^{(0)}$                   &  -5.82                         &     -20.03 \\
$\mathcal{W}^{(2)}$                   &  -1389.00                  &     -162.47 \\
$\mathcal{W}^{(4)}$                   &   -111467.25            &     65344.75 \\
\end{tabular}
\end{ruledtabular}
\end{table}

The spin-resolved CECs and the non-orthogonality by our ${\mathbf k}\cdot{\mathbf p}$ model are shown in Fig.~\ref{fig2}. As seen in the figure, the effective model underestimates the non-orthogonality for the lower-energy Dirac surface states. A better agreement with the respective \textit{ab initio} results is achieved by increasing the magnitude of $\gamma^{(5)}$  by a factor of 4.5 ($\gamma^{(5)}=2382.3$~a.u.) and 3.8 ($\gamma^{(5)}=6596.4$~a.u.) for Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se, respectively, see the orange areas and curves in Fig.~\ref{fig2} (the respective energy bands are shown by orange lines in Fig.~\ref{fig1}). Note that by manually correcting this parameter we reproduce more accurately not only the non-orthogonality, but also the hexagonal warping of the contours.

\section{Effective fields and multiple winding}

\begin{figure}[tbp]
\centering
\includegraphics[width=\columnwidth]{Fig_3.png}
\caption{Effective spin-orbit magnetic fields in Eq.~(\ref{TwoBandEMF}) as a function of polar angle $\varphi_{{\mathbf k}}$.  The in-plane fields ${\bm{\mathcal{B}}}_{\rtm{R}}^{(1)}$,  ${\bm{\mathcal{B}}}^{(5)}$, and ${\bm{\mathcal{B}}}^{(7)}$ are represented by arrows showing the direction of the field vectors at a given ${\mathbf k}$. The out-of-plane field ${\bm{\mathcal{B}}}_{\rtm{Z}}^{(3)}$ is illustrated by the blurred sixty-degree sectors of green and orange color for the negative and positive sign of its $z$ projection, respectively.}
\label{fig3}
\end{figure}

We focus now on the terms of the $2\times2$ Hamiltonian that cause the non-orthogonality in the in-plane-spin structure. We rewrite the Hamiltonian~(\ref{UPTO7}) in terms of the Pauli matrices:
\begin{equation}\label{TwoBandHam}
  H^{2\times2}_{\rtm{\mathbf{kp}}} = \mathcal{E}(\mathbf{k})\sigma_0 + {\bm{\mathcal{B}}}\cdot{\bm\sigma},
\end{equation}
where $\mathcal{E}(\mathbf{k}) = \epsilon_1 + \mathcal{M}k^2+2\mathcal{N}k^6\cos6\varphi_{{\mathbf k}}$ represents the dispersion of the doubly degenerate bands with the hexagonal warping of their CECs. The SOI-induced splitting of the bands
\begin{equation}\label{HamEig}
E^{\pm}(\mathbf{k}) = \mathcal{E}(\mathbf{k}) \pm |\bm{\mathcal{B}}|
\end{equation}
is due to the Zeeman-like term with the effective (spin-orbit) magnetic field
\begin{equation}\label{TwoBandEMF}
{\bm{\mathcal{B}}} = \widetilde{\alpha}{\bm{\mathcal{B}}}_{\rtm{R}}^{(1)} + 2\mathcal{W}{\bm{\mathcal{B}}}_{\rtm{Z}}^{(3)} + \widetilde{\gamma}{\bm{\mathcal{B}}}^{(5)} +\xi{\bm{\mathcal{B}}}^{(7)}.
\end{equation}
This field consists of the classical (linear) Rashba magnetic field ${\bm{\mathcal{B}}}_{\rtm{R}}^{(1)} =  k(\sin\varphi_{{\mathbf k}},-\cos\varphi_{{\mathbf k}},0)$, the cubic field ${\bm{\mathcal{B}}}_{\rtm{Z}}^{(3)} =  k^3(0,0,\sin3\varphi_{{\mathbf k}})$ responsible for the well-known three-fold  symmetric pattern of the spin $z$ component and contributing to the hexagonal warping of the CECs~[\onlinecite{Fu_PRL_2009}], and two higher-order six-fold symmetric fields ${\bm{\mathcal{B}}}^{(5)} =  k^5(\sin5\varphi_{{\mathbf k}},\cos5\varphi_{{\mathbf k}},0)$ and ${\bm{\mathcal{B}}}^{(7)}  =  k^7(\sin7\varphi_{{\mathbf k}}, - \cos7\varphi_{{\mathbf k}},0)$, see Fig.~\ref{fig3}. Note that since the spin matrix $(s^{\shortparallel}_1\bm{\sigma}_{\shortparallel}, s^{z}_1\sigma_z)$ of our two-band ${\mathbf k}\cdot{\mathbf p}$ model [the upper-left $2\times2$ block of the spin matrix~(\ref{SpinTWOSS})] differs from the matrix $(\bm{\sigma}_{\shortparallel}, \sigma_z)$ of a model built on a scalar-relativistic basis only by the non-unity coefficients $s^{\shortparallel}_1$ and $s^{z}_1$, one is tempted to treat the Pauli matrices in Eq.~(\ref{TwoBandHam}) as if they were spin matrices. Then, the spin expectation value is $\mathbf{S}^{\pm}({\mathbf{k}})=\pm \frac{1}{2}\bm{\mathcal{B}}/ |\bm{\mathcal{B}}|$. However, irrespective of the interpretation of $\bm{\sigma}$ in Eq.~(\ref{TwoBandHam}), in our model the non-orthogonality is characterized by the deviation angle $\delta^{\pm}$ found from the dot product
\begin{equation}\label{delta_nonorth}
\sin\delta^{\pm} = \frac{\mathbf{S}^{\pm}_{\shortparallel}({\mathbf{k}})\cdot \mathbf{k}} {|\mathbf{S}^{\pm}_{\shortparallel}({\mathbf{k}})|k}  = \pm \frac{k^5}{ |\bm{\mathcal{B}}_{\shortparallel}|}(\widetilde{\gamma}+\xi k^2)\sin6\varphi_{\mathbf{k}},
\end{equation}
where the parallel (in-plane) component of the effective field~(\ref{TwoBandEMF}) is
\begin{eqnarray}\label{EMF_parallel}
|\bm{\mathcal{B}}_{\shortparallel}|^2&=&(\widetilde{\alpha}^2+\widetilde{\gamma}^2k^8+ \xi^2k^{12})k^2 \nonumber \\
&-& 2\widetilde{\alpha} k^6(\widetilde{\gamma}-\xi k^2)\cos6\varphi_{\mathbf{k}} \\
&-&2\widetilde{\gamma}\xi k^{12} \cos12\varphi_{\mathbf{k}}. \nonumber
\end{eqnarray}
We see that the angle $\delta^{\pm}$ has a non-trivial dependence on the polar angle with a $\pi/3$ periodicity due to the presence of the fields ${\bm{\mathcal{B}}}^{(5)}$ and ${\bm{\mathcal{B}}}^{(7)}$. Acting separately, these fields may cause a quintuple or a septuple winding of the in-plane spin, respectively, in contrast to the Rashba field ${\bm{\mathcal{B}}}_{\rtm{R}}^{(1)}$ yielding a single winding, Fig.~\ref{fig3}.

At a given $\mathbf{k}$, the importance of each contribution to the effective magnetic field of Eq.~(\ref{TwoBandEMF}) depends on the respective parameter of the Hamiltonian~(\ref{UPTO7}):  $\widetilde{\alpha}$, $\mathcal{W}$, $\widetilde{\gamma}$ or $\xi$. According to their values in Table~\ref{tab:table2}, the effect of the field ${\bm{\mathcal{B}}}^{(7)}$ is expected to be negligible, because $\xi$ is much smaller than $\alpha^{(7)}$ and $\gamma^{(7)}$. At the same time, the contribution of ${\bm{\mathcal{B}}}^{(5)}$ depends on the parameters $\gamma^{(5)}$ and $\gamma^{(7)}$, which are comparable to or even larger than $\alpha^{(5)}$ and $\alpha^{(7)}$, respectively. However, because of the dominant contribution of the linear Rashba field ${\bm{\mathcal{B}}}_{\rtm{R}}^{(1)}$ due to the rather large $\alpha^{(1)}_1$ in $\widetilde{\alpha}$ of Eq.~(\ref{TwoBandEMF}), the superposition of all the in-plane fields produces a single winding of the in-plane spin~[\onlinecite{Note1}].

The in-plane-filed contribution~(\ref{EMF_parallel}) as well as the out-of-plane contribution $|{\bm{\mathcal{B}}}_{Z}|^2 = 2 \mathcal{W}^2 k^6 (1-\cos6\varphi_{\mathbf{k}})$ of the effective field~(\ref{TwoBandEMF}) affects the eigenvalues~(\ref{HamEig}) of the Hamiltonian~(\ref{TwoBandHam}) though the splitting term $\pm|{\bm{\mathcal{B}}}|$. This means that the SOI-induced hexagonal warping of the CECs is due not only to the cubic field ${\bm{\mathcal{B}}}_{Z}$ as, e.g., in Ref.~[\onlinecite{Fu_PRL_2009}], but also to the fields ${\bm{\mathcal{B}}}^{(5)}$ and ${\bm{\mathcal{B}}}^{(7)}$, which contribute to the warping through the scalar products ${\bm{\mathcal{B}}}^{(5)}\cdot{\bm{\mathcal{B}}}_{\rtm{R}}^{(1)}$ and ${\bm{\mathcal{B}}}^{(7)}\cdot{\bm{\mathcal{B}}}_{\rtm{R}}^{(1)}$ [the terms proportional to $\widetilde{\alpha}\widetilde{\gamma}$ and $\widetilde{\alpha}\xi$ in Eq.~(\ref{EMF_parallel}), respectively]. As follows form Tables~\ref{tab:table1} and \ref{tab:table2}, the fields ${\bm{\mathcal{B}}}_{\shortparallel}$ and ${\bm{\mathcal{B}}}_{\rtm{Z}}$ are equally important for the $\cos6\varphi_{\mathbf{k}}$ distortion of the CEC. In addition, the hexagonal warping due to ${\bm{\mathcal{B}}}_{\rtm{Z}}$ gives rise to the $z$ spin component, so if one neglects the contribution of ${\bm{\mathcal{B}}}_{\shortparallel}$ and fits only the cubic field to calculated or measured CECs, one may arrive at a large out-of-plane spin polarization with the spin-momentum locking unaffected by the warping. In contrast, the CEC warping caused by the fields ${\bm{\mathcal{B}}}^{(5)}$ and ${\bm{\mathcal{B}}}^{(7)}$ is accompanied by a change of the locking angle between spin and momentum. This explains why a stronger  warping may imply a larger non-orthogonality.

\section{Conclusions}

To summarize, within a fully \textit{ab initio} ${\mathbf k}\cdot{\mathbf p}$ perturbation approach we have developed a two-band effective ${\mathbf k}\cdot{\mathbf p}$ model for the surface states of the three-dimensional topological insulators Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se. The model includes terms to seventh order in $\mathbf{k}$ and provides a comprehensive description of the surface spin structure characterized by a non-orthogonality between the surface-electron spin and its momentum. In the ${\mathbf k}\cdot{\mathbf p}$ theory, the non-orthogonality that arises naturally in the \textit{ab initio} calculations is included in the effective models through the higher-order terms in \textbf{k}. Our $\mathbf{k}\cdot\mathbf{p}$ expansion builds on the eigenfunctions of the {\it ab initio} Hamiltonian, and, therefore, a term of a given order is unambiguously determined by the \textit{ab initio} spinor wave functions and, in contrast to a fitting method, does not depend on the presence of other terms. We have shown that the $k_{\pm}^5$- and $k_{\pm}^7$-terms represent effective spin-orbit magnetic fields with six-fold symmetric patterns on the two-dimensional momentum plane and, thereby, can lead to a non-orthogonality with the $\pi/3$ periodicity as a function of the polar angle of $\mathbf{k}$. For Bi$_2$Se$_3$ and Bi$_2$Te$_2$Se, we have found that the contribution of the $k_{\pm}^7$-term is rather small, and it is the $k_{\pm}^5$-term that causes the few-degree deviation of the actual spin direction from the classical orthogonality.

Finally, we would like to note that the derived two-band Hamiltonian is fully applicable to classical Rashba systems such as the Au(111) surface state. Here, the non-orthogonality appears to be negligibly small, albeit nonzero. A similar study for the giant Rashba spin-split conduction state of a single BiTeI trilayer reveals a substantial non-orthogonality, rather different for the inner and outer constant energy contour. In fact, the two-band Hamiltonian (\ref{UPTO7}) can be considered typical of hexagonal structures. Thus, the simplified picture that the in-plane spin and momentum are locked perpendicular to each other by spin-orbit interaction might overlook important features inherent in the spin-related phenomena at the surfaces and in 2D structures.

\begin{acknowledgments}
This work was supported by funding from the Department of Education of the Basque Government under Grant No.~IT1164-19 and by the Spanish Ministry of Science, Innovation and Universities (Project No.~FIS2016-76617-P).
\end{acknowledgments}


