  We report on microscopic tight-binding modeling of surface states in
Bi$_2$Se$_3$ three-dimensional topological insulator, based on a sp$^3$
Slater-Koster Hamiltonian, with parameters calculated from density functional
theory. The effect of spin-orbit interaction on the electronic structure of the
bulk and of a slab with finite thickness is investigated. In particular, a
phenomenological criterion of band inversion is formulated for both bulk and
slab, based on the calculated atomic- and orbital-projections of the
wavefunctions, associated with valence and conduction band extrema at the
center of the Brillouin zone. We carry out a thorough analysis of the
calculated bandstructures of slabs with varying thickness, where surface states
are identified using a quantitative criterion according to their spatial
distribution. The thickness-dependent energy gap, attributed to inter-surface
interaction, and the emergence of gapless surface states for slabs above a
critical thickness are investigated. We map out the transition to the
infinite-thickness limit by calculating explicitly the modifications in the
spatial distribution and spin-character of the surface states wavefunction with
increasing the slab thickness. Our numerical analysis shows that the system
must be approximately forty quintuple-layers thick to exhibit completely
decoupled surface states, localized on the opposite surfaces. These results
have implications on the effect of external perturbations on the surface states
near the Dirac point.
