THEORY



\(
\def\vb{{\bf}}
\def\bold#1{{\bf #1}}
\)

1. $ k\cdot p $ effective Hamiltonian

When we only care about wave vector $ {\vb K} $ around a specific wave vector $ {\vb k}_0 $ in the Brillouin zone, it can be well described by using the $ k\cdot p $ effective model. Suppose that $ \psi_{n{\vb K}}({\vb r}) $ is the wavefunction of $ n $-th band which satisfies Schrödinger equation with spin-orbit coupling (SOC) $ \hat{H}_B\psi_{n{\vb K}}({\vb r})=\epsilon_{n}({\vb K})\psi_{n{\vb K}}({\vb r}) $, where
\begin{equation}
\hat{H}_B=\frac{\hat{p}^2}{2m}+V({\vb r})+\frac{1}{2m^2c^2}\left(\vb{\hat{ s}}\times\nabla V({\vb r})\right)\cdot\hat{{\vb p}}
\end{equation}
is the Bloch Hamiltonian operator with SOC, and $ \hat{\vb p} $ is the momentum operator, $ V({\vb r}) $ is the potential in crystal, $ \hat{{\vb s}} $ is the spin momentum operator, $ m $ is the electron mass and $ c $ is the light speed. To expand the Hamiltonian at $ \vb k_0 $, introducing a transformation $ \phi_{n{\vb k}}({\vb r})=\psi_{n{\vb K}}({\vb r})e^{-i\vb k\cdot\vb r} $, where $ {\vb k}={\vb K}-{\vb k}_0 $ is the deviation from $ {\vb k}_0 $, the Schrödinger equation which $ \phi_{\vb k}({\vb r}) $ obeys is
\begin{equation}
\left(\hat{H}_B+\frac{\hbar^2k^2}{2m}+\hat{H}'({\vb k})\right)\phi_{n{\vb k}}({\vb r})=\epsilon_n({\vb K})\phi_{n{\vb k}}({\vb r})
\label{eq2}
\end{equation}
where $ \hat{H}'({\vb k})=\frac{\hbar}{m}{\vb k}\cdot{\vb \pi} $ is the first-order term, and $ \hat{\vb \pi}=\hat{\vb p}+\frac{1}{2mc^2}\left(\hat{{\vb s}}\times\nabla V({\vb r})\right) $ is the generalized momentum operator with SOC. We can take $ \hat{H}^{kp}=\hat{H}_B+\frac{\hbar^2k^2}{2m}+\hat{H}'({\vb k}) $ as the equivalent Hamiltonian of $ \hat{H}_B $ since the eigenvalues of them are all equal.

Suppose that we have got the eigenenergies $\{ \epsilon_n({\vb k}_0) \}$ and the eigenstates $ \{|{n({\vb k}_0)}\rangle\} $ of $ \hat{H}_B $. Then we can obtain the generalized momentum elements by $ {\vb \pi}_{mn}=\langle{m({\vb k}_0)}|{\hat{\vb \pi}}|{n({\vb k}_0)}\rangle $, thus the matrix elements of $ \hat{H}^{kp} $ can be obtained by
\begin{equation}
H^{kp-\text{all}}_{mn}({\vb k})=\left(\epsilon_n({\vb k}_0)+\frac{\hbar^2k^2}{2m}\right)\delta_{mn}+\frac{\hbar}{m}{\vb \pi}_{mn}\cdot{\vb k}
\label{kp-initial}
\end{equation}
Usually, we are aimed at several low-energy bands (the set of which is denoted as $\mathcal{A}$). The set of other bands is denoted as $ \mathcal{B} $. Then we can fold down the $H^{kp-\text{all}}_{mn}({\vb k})$ Hamiltonian into a subspace $\mathcal{A}$ via Löwdin partitioning. After two-order Löwdin partitioning, the effective $ k\cdot p $ Hamiltonian is transformed as
\begin{equation}
\begin{aligned}
H^{kp}_{\alpha\beta}({\vb k})=&\left(\epsilon_\alpha({\vb k}_0)+\frac{\hbar^2k^2}{2m}\right)\delta_{\alpha\beta}+\frac{\hbar}{m}{\vb \pi}_{\alpha\beta}\cdot{\vb k}+\frac{\hbar^2}{2m^2}\sum_{l\in \mathcal{B}}\sum_{ij}\pi^{i}_{\alpha l}\pi^{j}_{l\beta}k^ik^j\times \left(\frac{1}{\epsilon_\alpha({\vb k}_0)-\epsilon_l({\vb k}_0)}+\frac{1}{\epsilon_\beta({\vb k}_0)-\epsilon_l({\vb k}_0)}\right)
\end{aligned}
\label{eq-kp-origin}
\end{equation}
with $ i,j\in \{x,y,z\} $. Hereafter we use $\alpha, \beta\in\mathcal{A}$, $l\in\mathcal{B}$ and $m, n\in \mathcal{A}\cup \mathcal{B}$.

2. Zeeman's coupling

When a magnetic field is applied to the system, the momenta $ \hbar k^i $ should be replaced by $ (-i\hbar \partial^i+eA^i) $ according to Peierls substitution, where $ {\vb A} $ is the vector potential and $e$ (positively valued) is the elementary charge. Therefore, $ \hbar^2k^ik^j $ in the summation will be replaced by the sum of the gauge dependent term $ \frac{1}{2}\left\{-i\hbar \partial^i+eA^i,-i\hbar \partial^j+eA^j\right\} $ and gauge independent term $ \frac{1}{2}\left[-i\hbar \partial^i+eA^i,-i\hbar \partial^j+eA^j\right]=-\frac{i\hbar e}{2}\sum_k\epsilon^{ijk}B_k $, where $ {\vb B} $ is the magnetic induction intensity. In this case, the total Hamiltonian can be written by the summation of the $ k\cdot p $ effective Hamiltonian $ H^{kp}_{\alpha\beta} $ and Zeeman's coupling $ H^{Z}_{\alpha\beta} $, which are gauge dependent and gauge invariant, respectively. The $ k\cdot p $ effective Hamiltonian $ H^{kp}_{\alpha\beta} $ can be expressed by
\begin{equation}
\begin{aligned}
H^{kp}_{\alpha\beta}=&\epsilon_\alpha({\vb k}_0)\delta_{\alpha\beta}+\frac{\hbar}{m}{\vb \pi}_{\alpha\beta}\cdot\left(-i\nabla+\frac{e}{\hbar}{\vb A}\right)+\sum_{ij}M_{\alpha\beta}^{ij}\left(-i\partial^i+\frac{e}{\hbar}A^i\right)\left(-i\partial^j+\frac{e}{\hbar}A^j\right)
\end{aligned}
\label{2kp-1}
\end{equation}
where
\begin{equation}
\begin{aligned}
M_{\alpha\beta}^{ij}=&\frac{\hbar^2}{2m}\delta_{\alpha\beta}\delta_{ij}+\frac{\hbar^2}{4m^2}\sum_{l\in\mathcal{B}}\left(\pi^{i}_{\alpha l}\pi^{j}_{l\beta}+\pi^{j}_{\alpha l}\pi^{i}_{l\beta}\right)\times\left(\frac{1}{\epsilon_\alpha({\vb k}_0)-\epsilon_l({\vb k}_0)}+\frac{1}{\epsilon_\beta({\vb k}_0)-\epsilon_l({\vb k}_0)}\right)
\end{aligned}
\label{2kp-2}
\end{equation}
The Zeeman's coupling term can be expressed by
\begin{equation}
H^{Z}_{\alpha\beta}=\frac{\mu_B}{\hbar}\left({\vb L}_{\alpha\beta}+2{\vb s}_{\alpha\beta}\right)\cdot{\vb B}
\label{Zeeman}
\end{equation}
where
\begin{equation}
\begin{aligned}
L_{\alpha\beta}^{k}=& - \frac{i\hbar}{2m}\sum_{l\in\mathcal{B}}\sum_{ij}\epsilon^{ijk}\pi^{i}_{\alpha l}\pi^{j}_{l\beta}\times\left(\frac{1}{\epsilon_\alpha({\vb k}_0)-\epsilon_l({\vb k}_0)}+\frac{1}{\epsilon_\beta({\vb k}_0)-\epsilon_l({\vb k}_0)}\right)
\end{aligned}
\label{Zeeman-L}
\end{equation}
can be considered as the orbital contribution, $ {\vb s}_{\alpha\beta}=\langle{\alpha({\vb k}_0)}|{\hat{\vb s}}|{\beta({\vb k}_0)}\rangle $ are the spin elements, and $ \mu_B $ is the Bohr magneton.

3. $ {vasp2mat} $ : to compute matrix elements from DFT wavefunctions

The matrix elements of $\hat{\vb \pi}$, $\hat{\vb s}$, $\hat{T}$ and $\hat{R}$ are required in constructing the $ k\cdot p $ Hamiltonian and Zeeman's coupling under the DFT wavefunctions as shown in Sec. 1. In VASP, the projector augmented wave (PAW) potential is used, which is a combination and generalization of the linear augmented plane wave method as well as the pseudopotential. Therefore, we derive these matrix elements under the PAW wavefunctions, which are implemented in $ vasp2mat $ to compute the matrix elements. The relation of all electron wavefunction ($ |{n(\vb k_0)}\rangle $) and pseudo wavefunction ($ |{\widetilde{n}(\vb k_0)}\rangle $) in PAW is
\begin{equation}
|{n(\vb k_0)}\rangle=\mathcal{T}|{\widetilde{n}(\vb k_0)}\rangle,
\end{equation} where the linear transformation $ \mathcal{T} $ can be expressed as
\begin{equation}
\mathcal{T}=1+\sum_{a\mu\zeta}\left(|{\phi^a_{\mu}\zeta}\rangle-|{\widetilde{\phi}^a_{\mu}\zeta}\rangle\right)\langle{\widetilde{p}^a_{\mu}\zeta}|
\label{T-trans}
\end{equation}
The $ |{\phi^a_{\mu}\zeta}\rangle $ is the direct product of the real space all electron partial wavefunction $ \phi^a_{\mu}(\vb r) $ at the $ a $-th atom with a spinor wavefunction $ |{\zeta}\rangle (\zeta=\uparrow \text{or}\downarrow) $, $ |{\widetilde{\phi}^a_{\mu}\zeta}\rangle $ is the direct product of the real space pseudo partial wavefunction $ \widetilde{\phi}^a_{\mu}(\vb r) $ with a spinor wavefunction $ |{\zeta}\rangle $, and $ |{\widetilde{p}^a_{\mu}\zeta}\rangle $ is a projector wavefunction comprised of the direct product of a real space projector wavefunction $ \widetilde{p}^a_{\mu}({\vb r}) $ and a spinor wavefunction $ |{\zeta}\rangle$. $ \phi^{a}_{\mu}(\vb r) $'s are obtained by all-electron calculation for the reference atom. The pseudo partial wavefunctions $ \widetilde{\phi}^a_{\mu}(\vb r) $ are identical to $ \phi^{a}_{\mu}(\vb r) $ outside the augmentation sphere of the corresponding atom and are much softer than $ \phi^{a}_{\mu}(\vb r) $ inside the augmentation sphere. $ \widetilde{\phi}^{a}_{\mu}(\vb r) $'s provide a complete basis for the pseudo wavefunction $|{\widetilde{n}(\boldsymbol{k}_0)}\rangle$ inside the augmentation sphere.
The projector wavefunctions are defined in such a way, ${\it i.e.}$, $\langle \widetilde{p}^a_{\mu} | \widetilde{\phi}^{a}_{\mu'} \rangle=\delta_{\mu\mu'}$, that $\langle \widetilde{p}_{\mu}^a\zeta | \widetilde{n}(\boldsymbol{k}_0)\rangle$ gives the expanding coefficients of $| \widetilde{n}(\boldsymbol{k}_0)\rangle$ on $ \widetilde{\phi}^a_{\mu}(\vb r) $. Usually the projector wavefunctions are chosen to be zero outside the core radius. Therefore, outside the augmentation sphere the third and second terms in $\mathcal{T}$ cancel each other exactly, and inside the augmentation sphere the first and third terms cancel each other exactly. $\mathcal{T}$ leaves $| \widetilde{n}(\boldsymbol{k}_0)\rangle$ unchanged outside the augmentation sphere and maps it to the all-electron wavefunction inside the augmentation sphere. Both ${\phi}^a_{\mu}(\vb r) $ and $ \widetilde{\phi}^a_{\mu}(\vb r) $ are stored as a radial part times an angular part, which can be expressed as
\begin{equation}
\left\{
\begin{aligned}
\phi^{a}_{\mu}(\vb r)=Y^{m \mu}_{l \mu}(\widehat{\vb r-{\vb \tau}_a})R^a_{\mu}(\left|\vb r-{\vb \tau}_a\right|)\\
\widetilde{\phi}^{a}_{\mu}(\vb r)=Y^{m \mu}_{l \mu}(\widehat{\vb r-{\vb \tau}_a})\widetilde{R}^a_{\mu}(\left|\vb r-{\vb \tau}_a\right|)
\end{aligned}\right.
\end{equation}
where $ {\vb \tau}_a $ is the site of the $ a $-th atom, $Y^{m}_l$ are sphere harmonics, and $R^{a}_\mu$ are real functions. Here the hatted vectors represent the unit vector along the corresponding directions.
According to Eq. (\ref{T-trans}), the PAW matrix form of a local operator $ \hat{F} $ can be expressed by
\begin{equation}
\begin{aligned}
F_{mn}=&\langle{m(\vb k_0)}|{\hat{F}}|{n(\vb k_0)}\rangle=\langle{\widetilde{m}(\vb k_0)}|{\mathcal{T}^{\dagger}\hat{F}\mathcal{T}}|{\widetilde{n}(\vb k_0)}\rangle\\
=&\langle{\widetilde{m}(\vb k_0)}|{\hat{F}}|{\widetilde{n}(\vb k_0)}\rangle+\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle F^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\vb k_0)}\rangle
\end{aligned}
\label{matrix-f}
\end{equation}
where $ F^{a}_{\mu\zeta,\nu\zeta'} $ is the projection matrix of the operator $ \hat{F} $ in the $ a $-th atom's augmentation sphere, which is defined by
\begin{equation}
F^{a}_{\mu\zeta,\nu\zeta'}=\langle{\phi^a_{\mu}\zeta}|{\hat{F}}|{\phi^{a}_{\nu}\zeta'}\rangle-\langle{\widetilde{\phi}^a_{\mu}\zeta}|{\hat{F}}|{\widetilde{\phi}^{a}_{\nu}\zeta'}\rangle
\label{ele-f}
\end{equation}
The first term in $F^{a}_{\mu\zeta,\mu'\zeta'}$ gives the contribution from the all-electron wavefunction in the augmentation sphere, and the second term cancels the contribution from the pseudo wavefunction in the augmentation sphere that are counted in the first term of Eq. (\ref{matrix-f}).

In PAW, the plane wave is used to span the pseudo wavefunction $ \widetilde{\psi}_{n \vb k_0}(\vb r,\zeta)=\langle{\vb r,\zeta}|{\widetilde{n}(\vb k_0)}\rangle $, which is expressed by
\begin{equation}
\widetilde{\psi}_{n\vb k_0}(\vb r, \zeta)=\sum_{\vb G}c^{n\vb k_0}_{\zeta\vb G}e^{i({\vb k_0+\vb G})\cdot{\vb r}}
\label{plane-wave}
\end{equation}
where $ c^{n\vb k_0}_{\zeta\vb G} $ are the plane wave coefficients. The projection coefficients $ \langle{\widetilde{p}^a_{\mu}\zeta}|{\widetilde{n}(\vb k_0)}\rangle $ in the second term of Eq. (\ref{matrix-f}) have already been calculated in VASP. Therefore, to calculate the matrix of the local operator $ \hat{F} $, we need to calculate $\langle{\widetilde{m}(\vb k_0)}|{\hat{F}}|{\widetilde{n}(\vb k_0)}\rangle$ and $F^{a}_{\mu\zeta,\nu\zeta'}$ in Eq. (\ref{matrix-f}).

3.1 Generalized momentum matrix

The matrix of the generalized momentum $ \hat{\vb \pi} $ can be expressed by substituting $ \hat{\vb \pi} $ into $ \hat{F} $ in Eq. (\ref{matrix-f})
\begin{equation}
\begin{aligned}
&\vb \pi_{mn}=\langle{\widetilde{m}(\vb k_0)}|{\hat{\vb p}}|{\widetilde{n}(\vb k_0)}\rangle+\frac{1}{2mc^2}\langle{\widetilde{m}(\vb k_0)}|{\hat{\vb s}\times \nabla V}|{\widetilde{n}(\vb k_0)}\rangle+\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\vb p^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\vb k_0)}\rangle\\
+&\frac{1}{2mc^2}\sum_{a\mu\nu\zeta\zeta'}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\langle{\phi^a_{\mu}\zeta}|{\hat{\vb s}\times\nabla V}|{\phi^{a}_{\nu}\zeta'}\rangle\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\vb k_0)}\rangle-\frac{1}{2mc^2}\sum_{a\mu\nu\zeta\zeta'}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\langle{\widetilde{\phi}^a_{\mu}\zeta}|{\hat{\vb s}\times\nabla V}|{\widetilde{\phi}^{a}_{\nu}\zeta'}\rangle\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\vb k_0)}\rangle
\label{pi-1}
\end{aligned}
\end{equation}
Since the SOC effect is considered only within the augmentation spheres in VASP, where $ |{\widetilde{n}(\vb k_0)\rangle}=\sum_{a\mu\zeta}|{\widetilde{\phi}^a_{\mu}\zeta}\rangle\langle{\widetilde{p}^a_\mu\zeta}|{\widetilde{n}(\vb k_0)}\rangle $, thus making the second term and the fifth term in Eq. (\ref{pi-1}) cancel out. Therefore, the generalized momentum matrix can be simplified as
\begin{equation}
\begin{aligned}
\vb \pi_{mn}=&\langle{\widetilde{m}(\vb k_0)}|{\hat{\vb p}}|{\widetilde{n}(\vb k_0)}\rangle+\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\vb \pi'^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\vb k_0)}\rangle
\end{aligned}
\label{pi-simp}
\end{equation}
where
\begin{equation}
\begin{aligned}
\vb \pi'^{a}_{\mu\zeta,\nu\zeta'}=&\delta_{\zeta\zeta'}\left(\langle{\phi^a_\mu}|{\hat{{\vb p}}}|{\phi^a_\nu}\rangle-\langle{\widetilde{\phi}^a_\mu}|{\hat{{\vb p}}}|{\widetilde{\phi}^a_\nu}\rangle\right)+\frac{\hbar^2}{4mc^2}\vb \sigma_{\zeta\zeta'}\times\langle{\phi^a_\mu}|{\nabla V}|{\phi^a_\nu}\rangle
\end{aligned}
\label{pi'-matrix}
\end{equation}
The first term of $ \vb \pi_{mn} $ in Eq. (\ref{pi-simp}) can be calculated by
\begin{equation}
\langle{\widetilde{m}(\vb k_0)}|{\hat{\vb p}}|{\widetilde{n}(\vb k_0)}\rangle=\sum_{\vb G}\hbar(\vb k_0+\vb G)c^{*m\vb k_0}_{\zeta\vb G}c^{n\vb k_0}_{\zeta\vb G}
\label{p-matrix}
\end{equation}
where $ \hat{\vb p}=-i\hbar\nabla $. The integrals in $ \vb \pi'^a_{\mu\zeta,\nu\zeta'} $ can be calculated by separating the radial part and angular part, which are expressed as
\begin{equation}
\left\{\begin{aligned}
&\begin{aligned}
\langle{\phi^a_\mu}|{\hat{{\vb p}}}|{\phi^a_\nu}\rangle&=-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\nabla Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2R^{a*}_\mu R^a_\nu-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\vb r}{r}Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2R^{a*}_{\mu}\partial_rR^{a}_\nu
\end{aligned}\\
&\begin{aligned}
\langle{\widetilde{\phi}^a_\mu}|{\hat{{\vb p}}}|{\widetilde{\phi}^a_\nu}\rangle=&-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\nabla Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2\widetilde{R}^{a*}_\mu \widetilde{R}^a_\nu-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\vb r}{r}Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2\widetilde{R}^{a*}_{\mu}\partial_r\widetilde{R}^{a}_\nu
\end{aligned}\\
&\langle{\phi^a_\mu}|{\nabla V}|{\phi^a_\nu}\rangle\approx\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\vb r}{r}Y^{m_\nu}_{l_\mu}\int\mathrm{d}rr^2R^{a*}_{\mu}\partial_rVR^{a}_\nu
\end{aligned}\right.
\label{pi'-matrix-ele}
\end{equation}
Finally, by substituting Eqs. (\ref{pi'-matrix}-\ref{pi'-matrix-ele}) into Eq. (\ref{pi-simp}), the matrices of the generalized momentum $ \hat{\vb \pi} $ can be obtained.

3.2 Spin matrices

By substituting the local operator $ \hat{\vb s} $ for $ \hat{F} $ in Eq. (\ref{matrix-f}), the corresponding matrix elements are given explicitly. The first term in Eq. (\ref{matrix-f}) can be calculated by
\begin{equation}
\begin{aligned}
\langle{\widetilde{\alpha}(\vb k_0)}|{\hat{\vb s}}|{\widetilde{\beta}(\vb k_0)}\rangle=&\frac{\hbar}{2}\sum_{\vb G\vb G'\zeta\zeta'}\vb \sigma_{\zeta\zeta'}c^{*\alpha\vb k_0}_{\zeta\vb G}c^{\beta\vb k_0}_{\zeta'\vb G'} \delta_{\vb G,\vb G'}
\label{spin-T-R-1}
\end{aligned}
\end{equation}
where $ \hat{\vb s}=\frac{\hbar}{2}\hat{\vb \sigma} $ and the projection matrix can be obtained by
\begin{equation}
\begin{aligned} \vb s^{a}_{\mu\zeta,\nu\zeta'}=&\frac{\hbar}{2}\delta_{l_\mu l_\nu}\delta_{m_\mu m_\nu}\vb \sigma_{\zeta\zeta'}\int\mathrm{d}rr^2\left(R^{a*}_\mu R^{a}_\nu-\widetilde{R}^{a*}_\mu \widetilde{R}^{a}_\nu\right)
\end{aligned}
\label{spin-T-R-2}
\end{equation}
Finally, by substituting Eqs. (\ref{spin-T-R-1},\ref{spin-T-R-2}) into Eq. (\ref{matrix-f}), the matrices of the spin $ \hat{\vb s} $ can be obtained in the PAW wavefunctions.

3.3 Space group operator matrices

The general (magnetic) space group operator (SGO) is expressed by $F=\{R|\vb v\}$ or $F=\{TR|\vb v\}$, which is usually a nonlocal (NL) operator. One should notice that the SGO commutes with $\mathcal{T}$ in Eq. (\ref{T-trans}). Thus the matrix element can be written as
\begin{align}
F^{NL}_{mn}=& \langle{m(\vb k_0)}|{\hat{F}}|{n(\vb k_0)}\rangle=\langle{\widetilde{m}(\vb k_0)}|{\mathcal{T}^{\dagger} \mathcal{T}}|{\hat{F} \cdot \widetilde{n}(\vb k_0)}\rangle \nonumber\\
=& \langle {\widetilde{m}(\vb k_0)} | \hat{F} \cdot \widetilde{n}(\vb k_0)\rangle+ \sum_{a\mu\nu\zeta}\langle{\widetilde{m}(\vb k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle C^{a}_{\mu,\nu}\langle{\widetilde{p}^{a}_{\nu}\zeta}|{ \hat{F} \cdot \widetilde{n}(\vb k_0)}\rangle,
\end{align}
where
\begin{equation}
\begin{aligned}
&C^a_{\mu,\nu}=\delta_{l_\mu l_\nu}\delta_{m_\mu m_\nu}\int\mathrm{d}rr^2\left(R^{a*}_\mu R^{a}_\nu-\widetilde{R}^{a*}_\mu \widetilde{R}^{a}_\nu\right)\\
&\langle{\vb r,\zeta}|{\left\{R|\vb v\right\}}|{\widetilde{\beta}(\vb k_0)}\rangle=\sum_Gc^{\beta\vb k_0}_{\zeta\vb G} e^{i({\vb k_0}+{\vb G})\cdot R^{-1}({\vb r}-{\vb v})}\\
&\langle{\vb r,\zeta}|{\left\{TR|\vb v\right\}}|{\widetilde{\beta}(\vb k_0)}\rangle=\eta_\zeta\sum_Gc^{\beta\vb k_0*}_{\overline{\zeta}\vb G} e^{-i({\vb k_0}+{\vb G})\cdot R^{-1}({\vb r}-{\vb v})}
\end{aligned}
\end{equation}
Here, $\eta_\uparrow=1$, $\eta_\downarrow=-1$, $\bar\zeta$ indicates the opposite of spin $\zeta$. In our convention, the translation {$\vb L$} is expressed by a phase factor of $e^{-i \vb k_0\cdot \vb L}$.