Explicit expressions of the density matrices were obtained (see SI for detailed forms) by using the python-based SecondQuantizationAlgebra library and Wick's theorem. Finally, the matrix elements of the SOMF Hamiltonian can be obtained by the contraction of the transformed density matrices with the effective one-electron integrals in Eq. 10. The complex-valued Hermitian matrix is then diagonalized to give eigenvalues and eigenstates of the SOC-included Hamiltonian operator.