ψ$_{αβ}$$_{1}$ =    D$_{+}$λ 2 α − L$_{2}$ + E − iA$_{1}$ ( βλ$_{α}$ ) 0 A$_{2}$k$_{+}$    , (7) ψ$_{αβ}$$_{2}$ =    0 iA$_{1}$ ( βλ$_{α}$ ) D$_{−}$λ 2 α − L$_{1}$ + E    . (8) Ψ( E, k, z ) = ∑ α =1 , 2 ∑ β = ± ∑ γ =1 , 2 C$_{αβγ}$ψ$_{αβγ}$e βλ α $^{z}$, (9) F = A 2 1 + D$_{+}$ ( E − L$_{1}$ ) + D$_{−}$ ( E − L$_{2}$ ) , R = F 2 − 4 D$_{+}$D$_{−}$ [( E − L$_{1}$ )( E L$_{2}$ ) − A 2 $_{2}$k$_{+}$ k$_{−}$ ] , D$_{±}$ = D$_{1}$ ± B$_{1}$, L$_{1}$ = C + M + ( D$_{2}$ − B$_{2}$ ) k $^{2}$, L$_{2}$ = C − M + ( D$_{2}$ + B$_{2}$ ) k $^{2}$. (6) λ$_{α}$ ( E ) = [ − F 2 D$_{+}$D$_{−}$ + ( − 1) α − 1 √$_{R}$ 2 D$_{+}$D$_{−}$ ] 1 2 , (5) det | H ( k, − iλ ) − E | = 0 (4) H ( k, − i∂$_{z}$ ) ψ = Eψ, (3) ψ = ψ$_{λ}$e λz (2) 2 II. MODEL AND GENERAL SOLUTIONS FOR 3DTI thin films of Bi$_{2}$Se$_{3}$ are in the QSH regime. Finally, a conclusion is presented in Sec. VI. putting a four-component trial solution into the Schr¨ odinger equation ( E is the eigenvalue of en- ergy) the secular equation gives four solutions of λ ( E ), denoted as βλ$_{α}$ ( E ), with α ∈ { 1 , 2 } , β ∈ { + , −} , and where for convenience we have defined A. Model for 3DTI FIG. 1: Schematic of a topological insulator film grown on substrate. The grown direction is defined as z axis. The thickness of the film is L . As shown in Fig. 1, we will consider a thin film grown along z direction. The thickness of the film is L . We assume translational symmetry in x - y plane so that the wave numbers k$_{x}$ and k$_{y}$ are good quantum numbers. We start with the effective model proposed to describe the bulk states near the Γ point for the bulk Bi$_{2}$Se$_{3}$ $^{19}$. The states are mainly contributed by four hybridized states of Se and Bi p$_{z}$ orbitals, denoted as {| p 1 + z , ↑〉 , | p 2 − z , ↑〉 , | p 1 + z , ↓〉 , | p 2 − z , ↓〉} , where + ( − ) stands for the even (odd) parity. The Hamiltonian is given by Because of double degeneracy, each of the four βλ$_{α}$ ( E ) corresponds to two linearly independent four-component vectors, found as H ( k ) = ϵ$_{0}$ ( k ) I$_{4}$$_{×}$$_{4}$ +    M ( k ) − iA$_{1}$∂$_{z}$ 0 A$_{2}$k$_{−}$ − iA$_{1}$∂$_{z}$ −M ( k ) A$_{2}$k$_{−}$ 0 0 A$_{2}$k$_{+}$ M ( k ) iA$_{1}$∂$_{z}$ A$_{2}$k$_{+}$ 0 iA$_{1}$∂$_{z}$ −M ( k )    , (1) The general solution should be a linear combination of these eight functions where k$_{±}$ = k$_{x}$ ± ik$_{y}$ , ϵ$_{0}$ ( k ) = C − D$_{1}$∂ 2 z + D$_{2}$k $^{2}$, M ( k ) = M + B$_{1}$∂ 2 z − B$_{2}$k $^{2}$, and k 2 = k 2 $_{x}$+ k 2 $_{y}$, with A$_{1}$ , A$_{2}$ , B$_{1}$ , B$_{2}$ , C , D$_{1}$ , D$_{2}$ , and M the model parameters. This model has the time reversal symmetry and the inversion symmetry. Though we start with a concrete model, the conclusion in this paper should be applicable to other topological insulator films. We shall demonstrate that this model for the bulk states can produce the surface states with appropriate boundary condition. with the superposition coefficients C$_{αβγ}$ to be determined by boundary conditions. In the following, we will con- sider two different boundary conditions: one is semi- infinite focusing on only one surface at z = 0; the other includes two opposite surfaces at z = ± L/ 2. In both cases we assume open boundary conditions (Ψ = 0) for the surface states at the surfaces. Following the method by Zhou et al. $^{33}$, the general solution for either the bulk states or the surface states can be derived analytically. Despite the existence of time- reversal symmetry, the A$_{2}$k$_{±}$ term couples opposite spins in Hamiltonian (1), and one has to solve a 4 × 4 matrix, instead of the simplified 2 × 2 one in the 2D case $^{33}$. By The surface states have a finite distribution near the boundary. For a film thick enough that the states at op- posite surfaces barely couple to each other, we can focus C. Solutions for the surface states with semi-infinite boundary conditions B. General solutions of the surface states