E ﬀ ect of Rashba Impurities on Surface State of a Topological Kondo Insulator

: In this communication, we report surface state, with Rashba impurities, of a generic topological Kondo insulator (TKI) system by performing a mean-ﬁeld theoretic (MFT) calculation within the framework of slave-boson protocol. The surface metallicity together with bulk insulation is found to require very strong f -electron localization. The possibility of intra-band as well as inter-band unconventional plasmons exists for the surface state spectrum. The paramountcy of the bulk metallicity, and, in the presence of the Rashba impurities, the TKI surface comprising of ‘helical liquids’ are the important outcomes of the present communication. The access to the gapless Dirac spectrum leads to spin-plasmons with the usual wave vector dependence q 1 / 2 . The Rashba coupling does not impair the Kondo screening and does not a ﬀ ect the quantum critical point (QCP) for the bulk.


Introduction
The Periodic Anderson model (PAM) [1][2][3][4][5][6] basically involves highly correlated electrons (localized magnetic moments) in one f -orbital, which are screened by weakly correlated electrons in a second d orbital. Although this model has been thoroughly examined for several decades, the model itself and its extensions [5,6] are still relevant for the theoretical condensed-matter physics. A captivating development in recent years is the discovery that Kondo insulators can develop topological order to form a topological Kondo insulator (TKI). The theoretical description of a Kondo insulator [7][8][9][10][11] is usually based on PAM. In this paper we shall focus on a generic topological Kondo insulator with the simplest band structure. Upon including the effect of a Rashba system (such as Au(111), Ag(111), or Ni(111)) on the TKI surface, which bolster the Rashba spin-orbit interaction (RSOI) between f -electrons in the surface state Hamiltonian, we get access to the gapless Dirac spectrum and the spin-plasmons with the usual wave vector dependence q 1/2 . The Kondo screening in the bulk remains unaffected by RSOI.
It may be mentioned that the first example of correlated topological Kondo insulator (TKI) [1,7-10] is SmB 6 . It has a cubic crystal structure of the cesium-chloride type with a lattice constant of a ≈ 4.13 Å. At high temperature, the material behaves similar to a metal, but when reducing the temperature below~50 K, it exhibits insulating behavior [12]. The hybridization gap (∆) has been measured at 15 meV [12,13]. This is a peculiarity exhibited by a Kondo insulator. The compound SmB 6 exhibits yet another oddity. It has been found [14] that its resistivity at temperatures above 5 K increases upon decreasing the temperature. The saturation value in the resistivity, unaffected by the presence of non-magnetic impurities due to protection accorded to the surface states by the time reversal symmetry (TRS) [15], is achieved at temperatures below 5 K. A TRS breaking perturbations, such as the presence of the magnetic impurities, are expected to destroy topological surface state. Quite surprisingly, contrary to this expectation, the surface state in SmB 6 is found to be insensitive to the magnetic field or magnetic impurities by Cooley et al. [16]. It had been shown by them that the application of large magnetic fields up to 60 T is insufficient to destroy the low temperature saturation of resistance.

Periodic Anderson Model
The widely accepted model of a topological Kondo insulator (TKI) [7][8][9][10][11] involves, alongside a strong spin-orbit coupling, hybridization between an odd-parity nearly localized band and an even-parity delocalized conduction band. In the case of SmB 6 with a cubic crystal structure, these bands correspond to 4f and 5d electrons, respectively. It is then imperative that we start with the periodic Anderson model (PAM) where there are two different species of electrons, namely conduction electrons and localized electrons, often originating from d and f orbitals, respectively. We shall call them the conduction electrons and the valence electrons, respectively. The model [1][2][3][4][5][6]8,24] given below ignores the complicated multiplet structure of the d and f orbitals usually encountered in real TKIs such as SmB 6 . This has no major implication as most topological properties of cubic Kondo insulators do not depend on the precise form of the hopping and hybridization matrix elements or, on the particular shape of the orbitals.

Model
We consider below a well-known model [1][2][3][4][5][6]8] for a (topological) Kondo insulator on a simple cubic lattice with one spin-degenerate orbital per lattice site each for d and f -electrons. In momentum-space, we represent them by creation (annihilation) operators d † k ζ (d k ζ ) and f † k ζ (f k ζ ), respectively. Here, the index ζ (= ↑,↓) represents the spin or pseudo-spin of the electrons, where the latter is relevant for the localized electrons with generally large SOC. The Hamiltonian consists of two parts, namely, the bare hopping of the individual orbitals plus the hybridization between d and f orbitals (ℵ), and an onsite repulsion or f -electrons (ℵ int ). We have where d k = [2t d1 c 1 (k) + 4t d2 c 2 (k) + 8t d3 c 3 (k)] and f k = − f + 2t f 1 c 1 (k) + 4t f 2 c 2 (k) + 8t f 3 c 3 (k) . The t d1 , t f 1 , t d2 , t f 2 , and t d3 , t f 3 , respectively, are the NN, NNN, and NNNN hopping parameters. In addition, c 1 (k) = (cos k x a + cos k y a + cos k z a), c 2 (k) = (cos k x a cos k y a + cos k y a cos k z a + cos k z a cos k x a), c 3 (k) = (cos k x a cos k y a cos k z a) With a as the lattice constant. The sums run over all values for the crystal momentum (k) and the index ζ = (↑,↓) for all (cubic) lattice sites. The dispersion of the d and f -electrons is described by the first and the second term, respectively. The momentum-dependent form-factor matrix Γ ζ = ↑,↓ (k), in the third term, corresponds to the Dirac-type spin-orbit coupling associated with the hybridization between the conduction and the f -electrons [1,7,21]. The matrix is given by Γ ζ = ↑,↓ (k) = −2V(s(k).ζ), where ζ α are the Pauli matrices in physical spin space, and s(k) = sin k x a, sin k y a, sin k z a . What is the source of the momentum dependence of the matrix Γ ζ = ↑,↓ ? Due to spin-orbit coupling the f -states are eigenstates of the total angular momentum J, and hence hybridize with conduction band states with the same symmetry. This gives rise to the momentum-dependence of the form-factor Γ ζ = ↑,↓ (k). Furthermore, the hybridization is the source of non-trivial topology of the emergent bands. Since the f -and d-states have different parities, for the hybridization we must have odd parity: The negative sign of t f 1 is necessary for the band inversion, which induces the topological state [25,26]. The system shows the insulating as well as the metallic phases. The difference between the metallic and insulating phase is the sign of t f1 : It is positive for the metallic and negative for the insulating phase. The band width of the f -electrons is much smaller than the band width of the conduction electrons and we therefore assume that |t f 1 | << |t d1 | and similar relations hold for second-and third-neighbor hopping amplitudes. The hybridization is characterized by the parameter V for which we typically use |V| < |t d1 |. Throughout the whole paper, we choose t d1 to be the unit of energy, t d1 = 1. We note that the hybridization is an odd function of k in order to preserve time reversal symmetry (TRS), as it couples to the physical spin of the electron. The interaction Here we have assumed that the f -electrons locally interact via a Hubbard-U repulsion while the d electrons are non-interacting.
The total Hamiltonian (without the interaction term U) yields the single-particle spectrum ∈ (ζ) for upper band (lower band), ζ = ±1 labels the eigenstates ( ↑, ↓ ) of ζ z , and ∈ h = −2V s 2 x + s 2 y + s 2 z 1/2 . The pictorial depiction of the spectrum is shown in Figure 1. It is clear from the figure that the bulk metallicity is accessible when t f 1 is positive, as the conduction bands in this case are half empty. On the other hand, when t f 1 is negative, the band gap exists between the valence and the conduction bands leading to the bulk insulation.
Surfaces 2020, 3 FOR PEER REVIEW 4 interacting. The total Hamiltonian (without the interaction term U) yields the single-particle spectrum ∈ ( ) ( ) = − + + ∈ where α = 1 (−1) for upper band (lower band), ζ = ±1 labels the eigenstates (↑, ↓) ofζ z, and∈ = −2 ( + + ) 1/2 . The pictorial depiction of the spectrum is shown in Figure 1. It is clear from the figure that the bulk metallicity is accessible when is positive, as the conduction bands in this case are half empty. On the other hand, when is negative, the band gap exists between the valence and the conduction bands leading to the bulk insulation. Our main aim in this paper is to capture physics associated with the inclusion of the RSOI generating impurities in the system, modelled by appropriate surface Hamiltonian derived from the bulk Hamiltonian ℵ + ℎ + ℵ . The interaction term ℵ can be studied nonperturbatively by various methods, including dynamical mean-field theory, Gutzwiller-projected variational wave-functions, or slave-particle representations [1,27,28]. In what follows we shall use the slave-particle protocol. The f-electron correlations, being larger than all other energy scales in the problem, effectively enforce a no-double-occupancy constraint on each site, and therefore contribute through virtual processes at low temperatures only. The original Hilbert space thus get projected onto a smaller subspace, where the double occupancy is excluded.

Slave-Boson Protocol
One utilizes the well-known slave-boson protocol [1,20,21,27,28] to do the projection onto the Hilbert space. In this protocol, the electron operator is expressed in terms of psuedo-fermions and slave-bosons. The operators a † (heavy slave-boson creation operator) and b † (light slave-boson creation operator), respectively, creates doubly-occupied and empty bosonic impurity states out of vacuum. The bosons created by b † are supposed to carry the electron's charge. The single site fermionic occupation operator is denoted by s † . This corresponds to a fermionic operator and carry the electron's spin, but no charge. For U → ∞ , the double occupancy is prohibited, and hence the  Figure 1. The two-band energy spectra of a bulk TKI. The parameter values are, ak y = π, t d1 = 1, t d2 = 0.01, t d3 = 0.001, t f 2 = 0.01, t f 3 = 0.001, u = 0.0, f = −0.02, V = 0.10. In (a), t f 1 = 0.3. The conduction band is partially empty which implies the bulk metallicity. In (b), t f 1 = −0.3. In case this case, the Fermi level (chemical potential represented by solid, horizontal line) is located in the hybridization gap, which implies bulk insulation.
Our main aim in this paper is to capture physics associated with the inclusion of the RSOI generating impurities in the system, modelled by appropriate surface Hamiltonian derived from the bulk Hamiltonian ℵ + Rashba coupling + ℵ int . The interaction term ℵ int can be studied non-perturbatively by various methods, including dynamical mean-field theory, Gutzwiller-projected variational wave-functions, or slave-particle representations [1,27,28]. In what follows we shall use the slave-particle protocol. The f -electron correlations, being larger than all other energy scales in the problem, effectively enforce a no-double-occupancy constraint on each site, and therefore contribute through virtual processes at low temperatures only. The original Hilbert space thus get projected onto a smaller subspace, where the double occupancy is excluded.

Slave-Boson Protocol
One utilizes the well-known slave-boson protocol [1,20,21,27,28] to do the projection onto the Hilbert space. In this protocol, the electron operator is expressed in terms of psuedo-fermions and slave-bosons. The operators a † (heavy slave-boson creation operator) and b † (light slave-boson creation operator), respectively, creates doubly-occupied and empty bosonic impurity states out of vacuum. The bosons created by b † are supposed to carry the electron's charge. The single site fermionic occupation operator is denoted by s † . This corresponds to a fermionic operator and carry the electron's spin, but no charge. For U → ∞ , the double occupancy is prohibited, and hence the operators a † and a will not be under consideration. The connection of the remaining auxiliary operators (b, s) to the physical f -electron operator is f † ζ = s † ζ b. In the slave-boson mean-field theory of the infinite-U model (which consists of replacing slave-boson field (b) at each lattice site by the modulus of its expectation value-a c number), the anti-commutation relation To take care of the conservation of auxiliary particle number one needs to impose the restriction ζ s † It is now straight-forward to write the thermal average of the TKI slave-boson mean-field Hamiltonian where the additional terms, in comparison with (1) The first term is the constraint which fixes the total number of particles N (N = N d + N s ), the second term enforces the fact that there are equal d and s fermions (N d = N s ), and the third term describes the constraint on the pseudo-particles due to the infinite Coulomb repulsion. In order to have a Kondo insulator, formation of singlet states between d and s fermions is needed at each lattice site. This means that the number of d and s fermions are equal on average. Here λ is a Lagrange multiplier, and N s is the number of lattice sites for s electrons (similarly, N d corresponds to d-electrons). The dispersion of the f -electron is renormalized by λ and its hopping amplitude by b 2 . Moreover, the hybridization amplitude is also renormalized by the c-number b. The total parameters are slave-boson field b, auxiliary chemical potentials ξ and µ (µ is a free parameter), and the Lagrange multiplier λ. Though these facts are explained clearly in [1], we have never-the-less found it necessary to include them to make this paper self-contained. One obtains equations for the parameters (b,λ,ξ) minimizing the thermodynamic potential per unit volume O sb = −(βV) −1 lnTr exp(−βℵ sb (b, λ, ξ)): ∂O sb ∂b = 0, ∂O sb /∂λ = 0, and ∂O sb /∂ξ = 0. Here, β denotes the inverse of the product of temperature T and Boltzmann constant k B . The thermodynamic potential has been calculated by the method outlined in [29] and [30] (see Appendix A). These equations are: The averages d † k.ζ d k,ζ , bs † k.ζ bs k,ζ , etc., have been calculated in Appendix A below. Their expressions show that in the zero-temperature and the long-wavelength limits, the contribution of the averages d † k↑ bs k,↓ and bs † k,↑ d k↓ , etc., to the derivatives in Equations (4) . This observation allows us to approximate the equations as in these limits. As we already have noted, for the conservation of auxiliary particle number, one needs s k e ik.r s k,ζ where N s is the number of s-fermion, equivalently, one may also write this as dr ζ bs † ζ (r)bs ζ (r) N s b 2 1 − b 2 . This is the fourth equation, with (9)-(11) as the first three, and we have four unknowns, viz. (b,λ,ξ,ak F ), where a ≈ 4.13 A is the latticeconstant and k F is the Fermi wave number. In view of the fact that the chemical potential is a free parameter and could be somewhere between the valence and the conduction bands once again it is easy to see that for the low-lying states in the zero-temperature limit one may write the imposed restriction as 2 N −1 s k,ζ 1 = b 2 − b 4 . This is an equation for b 2 in terms of (ak F ). After a little algebra, we find that whereas (9) and (10) We estimate (ak F ) in the following manner : Since the Fermi velocity v * F of the low-lying states is known to be less than 0.2 eV/Å [24], taking the effective fermion mass (m * = k F v * F ) a 100 times that of an electron [25], we find that (ak F ) ∼ 0.01. This is consistent with the long wavelength limit we have assumed. The two values of b 2 obtained from here are close to but less than one (1 − ) and close to zero (0 + ). The fact to remember is when b is nonzero, the system is in Kondo state, and when b vanishes, the system is in normal gas state. The admissible value of b 2 will be thus be 1 − . With this a reasonable spectral gap is obtained in the Kondo insulating state (see Figure 2b) in the long wavelength limit. We see in the figure that in order to have a gapped Kondo state it is required that the tunneling of the localized fermions must have an opposite sign compared with the conduction fermions. This is difficult to achieve, because the signs of the tunneling for the lowest Bloch bands are usually the same. However, it has been taken to be positive for the metallic phase and negative for the insulating phase for the following reason: For the metallic phase the leading term for f -electrons in the Hamiltonian, viz. −2t f 1 c 1 (k) bs † k.ζ bs k,ζ , corresponds to a minimum at (0,0,0) (electron-like band) and in the insulating phase it corresponds to a maximum (hole-like band). Furthermore, it must be mentioned that in their quantum simulation of the topological Kondo insulator in ultra-cold atoms, Zheng et al. [26] have suggested a transformation leading to a staggered Kondo coupling-positive for the metallic phase and negative for the insulating phase. Abiding by this specification, we find that in the metallic phase ξ is −2.7 and, in the insulating phase it is −3.3. All the unknown parameters now stand determined.
leading to a staggered Kondo coupling-positive for the metallic phase and negative for the insulating phase. Abiding by this specification, we find that in the metallic phase ξ is −2.7 and, in the insulating phase it is −3.3. All the unknown parameters now stand determined.
An investigation on Kondo system is incomplete without the possibility of the Kondo screening being explicitly shown. We, therefore, calculate the Kondo singlet density which may be defined as The calculation details is given in Appendix A. K singlet (k,µ) = _( , , ,µ, ) ( (∈ ( ) ( ) ) + 1 − ( (∈ ( ) ( ) ) + 1) ] Where , , ∈ ( ) , etc., are given by Equations. (A7)-(A10).This average is the ultimate signature of the Kondo insulating state, where there is precisely one conduction electron paired with an impurity spin. The point we wish to make is that unless we have an anti-ferromagnetic exchange field in the bulk Hamiltonian this Kondo screening term will not be non-zero. Assuming(V/ ) ≪ 1, we model the interaction between an impurity moment and the itinerant (conduction) electrons in the system with coupling term −| | ∑m Sm. Sm, where Sm is the m th-site impurity spin, s = , is the fermion annihilation operator at site-m and spin-state ζ (=↑,↓) and ζ z is the z-component of the Pauli matrices. We make the approximation of treating the impurity spins as classical vectors. The latter is valid for| | > 1. Absorbing the magnitude of the impurity spin into the coupling constant (M = | || |/ )it follows that the exchange field term, in the ( k↑ f k↑ d k↓ f k↓ ) T basis, appears as {ζ z ⊗M (τ0+τz)/2}, where τ0 and τz, respectively, are the identity and the z-component of An investigation on Kondo system is incomplete without the possibility of the Kondo screening being explicitly shown. We, therefore, calculate the Kondo singlet density which may be defined as The calculation details is given in Appendix A.
where ε + , ε − , ∈ (ζ) α , etc., are given by Equations (A7)-(A10). This average is the ultimate signature of the Kondo insulating state, where there is precisely one conduction electron paired with an impurity spin. The point we wish to make is that unless we have an anti-ferromagnetic exchange field in the bulk Hamiltonian this Kondo screening term will not be non-zero. Assuming(V/t d1 ) 1, we model the interaction between an impurity moment and the itinerant (conduction) electrons in the system with coupling term −| J | m S m . S m , where S m is the m th-site impurity spin, s m = 1 2 d † mζ ζ z d mζ , d mζ is the fermion annihilation operator at site-m and spin-state ζ (=↑,↓) and ζ z is the z-component of the Pauli matrices. We make the approximation of treating the impurity spins as classical vectors. The latter is valid for |S| > 1. Absorbing the magnitude of the impurity spin into the coupling constant J(M = | J ||S|/t d1 ) it follows that the exchange field term, in the (d k↑ f k↑ d k↓ f k↓ ) T basis, appears as {ζ z ⊗M (τ 0 +τ z )/2}, where τ 0 and τ z , respectively, are the identity and the z-component of Pauli matrix for the pseudo-spin orbital indices. We thus obtain the dimensionless contribution [M k,ζ sgn(ζ) d † k,ζ d k,ζ ] to the momentum space Hamiltonian above. It is not difficult to see that the terms ε + , ε − , and ∈ (ζ) α are to be redefined now in the following manner: (16) Upon substituting these re-defined terms in Equation (12) we obtain non-zero values of K singlet . In Figure 2  . This is typical Doniach-like phase diagram [31]. At T = 0 K M > M QCP corresponds to the heavy-fermion liquid, while M < M QCP corresponds to anti-ferromagnetic liquid. From the plots we notice that the location of QCP depends on µ and t f 1 . In fact, QCP increases with increase in µ and decreases with increase in |t f 1 |.
In order to ascertain whether there is any effect of Rashba coupling (RC) (λ R ) on the Kondo screening, we add the term [2λ R k k y + ik x bs † k↑ bs k,↓ + H.C.] to the bulk Hamiltonian. We have assumed to have deposited particles with high Rashba spin-orbit (RSO) interactions λ R , such as Au(111), on the bulk of material. We once again calculate the average S = [2Vb s x − is y ( d † k↑ bs k,↓ + bs † k,↑ d k↓ ) + H.C.] corresponding to the Kondo singlet following the method outlined in the Appendix A. We find that the RC does not impair the Kondo screening and does not affect the QCP for the bulk. Even in the surface state Hamiltonian in the next section, upon adding the Rashba coupling, one can show that the average similar to S is non-zero. The stage is now set to investigate the surface state.

Surface States
We choose a representation involving the states (| c k,↑,± , | c k,↓,± ) , where the operators c k,↑,± and c k,↓,± , respectively, are the spin-up and spin-down annihilation operators corresponding to the upper and lower bands with spin-splitting (∈ (↑) ± (k) and ∈ ( ↓) ± (k)), as our basis. In the basis chosen above our bulk mean-field Hamiltonian matrix in slave-boson protocol will appear as The eigenvalues of this matrix are given by The energy eigenvectors corresponding to these eigenvalues ∈ (ζ) α (k) are given by α (k). Our aim in this section is to examine the surface plasmons of the system under consideration. The plasmons are defined as longitudinal in-phase oscillation of all the carriers driven by the self-consistent electric field generated by the local variation in charge density on the surface. We require surface state single particle excitation spectrum for this purpose. To study the surface states, we first consider a thick slab limited in z ∈ [−d/2, d/2], with open boundary conditions, where d is the thickness of the slab in z direction. In this case k z is not a good quantum number, and should be replaced by −i∂ z . The Hamiltonian ℵ slab k x , k y for the slab geometry under consideration can be obtained from the Hamiltonian ℵ bulk k x , k y , − i ∂ z considering the orthonormal function | ϕ n (z) = ψ_n(z), where With this ansatz, the matrix elements may be written as The eigenstates of ℵ mn slab k x , k y are the so called edge states. Equation (17) yields ℵ mn slab as where Since two blocks involve zeros in Equation (23), we definitely require a better version of surface state Hamiltonian by guessing an improved ansatz. The work is under way in this direction. Now upon noting that λ = −6t f 1 + 6b 2 t f 1 and ξ ≈ −3t d1 , the eigenvalues of this matrix is given by where α = +1 (−1) for conduction (valence) band, and K = k x , k y . The term v * F (ζ) will be interpreted as the Fermi velocity where v * 3d and 96a 7d , respectively, for the up and down quasi-particle spins. The graphical representation of Equation (25) is shown in Figure 3. The thin curves (solid and dashed) and the thick curves in both the figures correspond to spin-split conduction and valence bands. The chemical potential is taken to be zero and represented by a thick horizontal straight line in this figure.
The eigenvectors corresponding to the eigenvalue ∈ and those corresponding to the eigenvalue ∈ (ζ = −1) α,sur f ace (K) are where (N 1 , N 2 ) correspond to the normalization terms. We notice that if the f -electrons somehow ) small compared to t d1 − b 2 t f 1 × (aK) 2 , one obtains gapless surface state bands at (0,0) wavevector for the insulating bulk. For this to happen, the severe restrictions are with b 2 = 0.9890 i.e., | t f 1 | should be less than | ε f | by one order of magnitude an the high symmetry point (0,0) is unapproachable as 1. Here , and | ε f | = 50 meV, we find the estimated value v * F ≈ 10 4 m/s. The hybridization parameter V cannot contribute here as we have chosen an ortho-normal function to obtain the Hamiltonian for the slab structure fromthe bulk Hamiltonian. It must be noted that Legner et al. [27] have found an expression of the Fermi velocity which with the same set of parameter values as above and V = 100 meV yields v * F nearly the same value as obtained by us. (In fact, they derived the surface state dispersion and found that the Fermi velocity for the electrons is v F = 4|V | (|t f1 t d1 |(t d1 − t f1 ) 2 ) 1/2 . This result implies that the effective mass of the surface electrons m * = p F /v F is quite heavy since the hybridization amplitude is small compared to other relevant energy scale, while the expression under the square root is of the order O(1)).
The Kondo screening length ξ (=¯h v * F k B T K,s ) for the surface states will be, therefore, be one order of magnitude higher than the lattice constant(a) for the surface Kondo temperature T K,s ∼ 25 K. We note that, with the exchange field introduced in the previous section, the eigenvalues of the matrix (23) is given by Since the conduction bands are partially empty, the surface state will be metallic in both (a) and (b). In (a), the hopping integral t f 1 = 0.31 is positive and therefore it corresponds to metallic bulk (Figure 1a). In (b), t f 1 = −0.31 is negative and therefore corresponds to insulating bulk (Figure 1b). The thin lines (solid and dashed lines) and thick lines in both the figures correspond to spin-split conduction and valence bands (see notes below Equation (25)).
We note that, with the exchange field introduced in the previous section, the eigenvalues of the matrix (23) is given by Surfaces 2020, 3

495
where v * F (ζ) = t d1 − b 2 t f 1 3t f 1 + ζM 2 − b 2 2 ε f will now be interpreted as the Fermi velocity. It is thus demonstrated that the exchange field can be used to open an additional gap at the surface state dispersion.
We now consider the problem of the TKI surface plasmonics. The local variation in charge density in materials gives rise to an electric field. The quanta of in-phase longitudinal density oscillation of charge carriers at the surface of materials driven by this electric field are defined as the surface plasmons (SPs). Inside the bulk, SPs evanesce severely owing to the heavy energy loss. These collective density oscillations can be excited in the conventional metal surfaces. It is clear from the plot of Equation (25) in Figure 3b that the surface could be metallic, even when the bulk is insulating, owing to the conduction bands being partially empty. This is, therefore, an appropriate case to investigate the possibility of the SPs. Equation (27) is expected to yield unconventional q 1/2 plasmons [32,33] obtainable under very stringent conditions given by Equation (28). Since the calculation of surface metallicity is in the zero-temperature limit, a prediction of the previous works [7 -10] stating that "when temperature is lowered a Kondo insulator may turn into a topological insulator with a metallic surface state" remains unverified. However, our finding of surface states with gapless Dirac dispersion albeit under very stringent conditions corroborates an important experimental finding of Xiang et al. [27]. They have found that, in the case of the prototypical TKI SmB 6 , there is a broken rotating symmetry in the amplitude of the main de Haas-van Alphen oscillation branch consistent with Lifshitz-Kosevich theory confirming a 2D nature of the surface electronic state. The finding is similar to quantum oscillation experiments [34,35] for the conventional TIs to probe the surface Dirac fermions. The transport measurements [35] in the past have also demonstrated the insulating bulk and metallic surface separation. For potential applications toward scalable quantum information processing [36], this bulk and surface separation is especially important.

Plasmon Frequency
The plasmons are defined as longitudinal in-phase oscillation of all the carriers driven by the self-consistent electric field generated by the local variation in induced charge density ρ (r,ω). In a linear-response approximation, we have ρ(r,ω) = e 2 d 2 r χ(r, r , ω)Φ(r , ω), where Φ is induced local potential and χ is the fermion response function or the dynamic polarization. This is a quantity of interest for many physical properties, since it determines e.g., the plasmon and phonon spectra. Assuming plasmon oscillation for the 3D system under consideration entirely a surface phenomenon, in the random phase approximation (RPA) [37,38], we write the dynamical polarization function χ (aq, ω) in the momentum space, as The symbol ∈ ζ,α,sur f ace (aK) stands for the surface state single-particle excitation spectrum given by Equation (25) and Ψ ζ,α (a(K − q)) Ψ ζ ,α (aK) 2 for the band-overlap of wave functions. This part of the paper leans on the previous investigations [32,33]. Since we are presently interested in intra-band plasmons only, we write down the explicit expression for the intra-band overlap. In view of (26), this is given by F α,α,ζ,ζ (K,q) = 1 2 1 + ζζ cos θ α,α,K,q , where cos θ α,α,K,q, = In the long wavelength limit and ζζ = +1. The Fermi function with the band index α is given by n ζ,α (aK) = 1/[exp(β (∈ ζ,α,sur f ace (aK)) − βµ) + 1]. Upon using the Sokhotski-Plemelj identity (x ± iη) −1 = P x −1 ∓ iπ(x) with P as the principal part, the real-part of the polarization function appears as The imaginary part is given by Since we have chosen t d1 to be the unit of energy before, the quantity ω is also in the same unit. In the long-wavelength limit, the band structure in Equation (24) yields where Within the random phase approximation (RPA), the plasmon dispersion is obtained by finding zeros of the dynamical dielectric function, which is expressed in terms of Coulomb's potential as e ζ,α a q , ω = 1 − V (q)χ ζ ,α a q , ω where ω = ω − iγ, γ is the decay rate of plasmons, and the expression V(q) is the Fourier transform of the Coulomb potential in two dimensions. For weak damping, the equation Re e ω, a q = 0 yields the plasmon frequency ω pl . In the intra-band case, from Equation (32), in the high frequency limit we obtain where V 0 = e 2 2ε 0 ε r , ε 0 is the vacuum permittivity, and ε r is the relative permittivity of the surrounding medium. The denominator of the summand in Equation (36) involves a scalar product. We make use of the standard integral Surfaces 2020, 3

497
The denominator of the summand in Equation (38) can be expanded using the Binomial theorem. We obtain then the following implicit equation as the Plasmon dispersion for the system: where dK = d 2 (aK) (2π) 2 and It may be noted that we have put ζζ = +1 in F α,α,ζ,ζ (K,q) = 1 2 1 + ζζ cos θ α,α,K,q and written it as ( 1 + (aq)FF K ) in Equation (39), as we are presently considering the intra-band case.
The integrals in Equation (39) where πV 0 (I 1 +(aq)I 2 ), in the long wavelength limit. To the leading order, I 1 and I 2 are positive as I n ∼ t d1 (aK F ) n+2 n+2 . It is to be noted that, though the relation is linear, the group velocity (v g ) and the phase velocity (v p ) are not equal. The former, given by v g = (0.6075 V 0 I 2 ) is several orders of magnitude smaller than the speed of light whereas the latter is positive and much larger than unity. Thus, considering the intra-band transitions only, we have found possibility of only one collective mode exhibiting linear dispersion which could be triggered into an excited state with very low levels of energy input-less than 1 electron-volt. It corresponds to charge plasmons and not spin plasmons as its origin does not lie in the gapless Dirac spectrum. The linear behaviour of the dispersion implies that signals can be transmitted undistorted along the surface. The finding has significant importance in extremely low loss communications.
As regards the inter-band plasmons corresponding to a pair of spin-split bands, we shall have ζζ = −1. The band-overlap of wave functions F α,α,ζ,ζ (K,q) = 1 2 1 + ζζ cos θ α,α,K,q may be written as 1 2 ( 1 − (aq)FF K ). Furthermore, the quantity ε r is the relative permittivity of the surrounding medium and it is positive (negative) for a meta-material (negative dielectric constant medium). We assume the surrounding medium to be a meta-material. Thus, we may write V 0 = −| V 0 |. Equation (43) in this case may be written as ω ≈ 3 − 2 1 3 9 πV 0 (−I 1 + (aq)I 2 ). As in the case of Surfaces 2020, 3 498 frequency dispersion in groups of gravity waves on the surface of deep water, we encounter the exotic possibility of the group velocity being positive, while the phase velocity is negative. Our analysis demonstrates that the same plasmonic system can support both type of solutions, viz. (v g , v p ) > 0 and (v g > 0, v p < 0), depending on parameters.

Surface Spectrum with Rashba Coupling
In the previous sub-section we have obtained the gapless Dirac spectrum Equation (27) under the severe restrictions of Equations (28a) and (28b). If this is true, we expect to have the well-known q 1/2 spin-plasmons [32,33]. Since the restrictions are ironclad, we look for other means to obtain gapless spectrum. None other than the spin-orbit coupling is more suitable for the purpose. We, therefore, propose the consideration of an additional (Rashba) term [2b 2 λ R k y + ik x + H.C] in the surface state Hamiltonian. It is imperative to assume that there must be deposition of particles with considerably high Rashba spin-orbit (RSO) interaction strength λ R , such as Au(111), on the surface of material. The effect of this intrinsic coupling is to be included in the band-structure only. To start with, we consider Equation (23) for ℵ mn slab and add Rashba term to it. We obtain ℵ mn modi f ted as where Upon noting that λ = −6t f 1 + 6b 2 t f 1 and ξ ≈ −3t d1 , the eigenvalues ε of this matrix is given by the quartic Moving over to Equation (47), we notice that the first term (z 0 /2) acts as an in-plane Zeeman term g b µ B B (where B is the magnetic field, and g b is the Lande g-factor, and µ B is the Bohr magneton). The pseudo-Zeeman term of the spectrum in Equation (47) comes into being due the presence of the term 2b 2 λ R 2 (ε − Γ 1 ) 2 K 2 in Equation (46). Without this term, the spectrum reduces to a bi-quadratic (with no Zeeman term) rather than a quartic. Thus, the role of Rashba SOI as the polarization-usherer could be easily understood.
We have plotted the surface state single-particle excitation spectrum Equation (47) as a function of the dimensionless wave vector (aK) in Figure 4 for the considerably high Rashba coupling ( The influences of the coupling on the spectrum is quite conspicuous if one compares this figure with Figure 3. The 'gaplessness at wave vector K = 0' is the striking feature of the spectrum in Figure 4a. In contrast, for the insulating bulk in Figure 4b, this feature is also manifested at K = ±1. The wider spin-splitting under the influence of the Rashba interaction are basically responsible for prominent gapless Dirac spectrum including the cones. The TKI surface, therefore, comprises of 'helical liquids' [39] in the presence of the Rashba impurities. The access to the gapless Dirac spectrum leads to spin-plasmons with the usual wave vector dependence q 1/2 . As we have noted above, the Rashba coupling does not impair the Kondo screening and does not affect the QCP for the bulk. Surfaces2020, 3 FOR PEER REVIEW 16 Moving over to Equation (47), we notice that the first term √(z0/2) acts as an in-plane Zeeman term gbμBB (where B is the magnetic field, and gb is the Lande g-factor, and μB is the Bohr magneton).
The pseudo-Zeeman term of the spectrum in Equation (47) comes into being due the presence of the term (2 ) (ε − ) in Equation (46). Without this term, the spectrum reduces to a bi-quadratic (with no Zeeman term) rather than a quartic. Thus, the role of Rashba SOI as the polarization-usherer could be easily understood.
We have plotted the surface state single-particle excitation spectrum Equation (47) as a function of the dimensionless wave vector (aK) in Figure 4 for the considerably high Rashba coupling ( ~ 0.5). The influences of the coupling on the spectrum is quite conspicuous if one compares this figure with Figure 3. The 'gaplessness at wave vector K = 0' is the striking feature of the spectrum in Figure 4a. In contrast, for the insulating bulk in Figure 4b, this feature is also manifested at K = ±1. The wider spin-splitting under the influence of the Rashba interaction are basically responsible for prominent gapless Dirac spectrum including the cones. The TKI surface, therefore, comprises of 'helical liquids' [39] in the presence of the Rashba impurities. The access to the gapless Dirac spectrum leads to spin-plasmons with the usual wave vector dependence q 1/2 . As we have noted above, the Rashba coupling does not impair the Kondo screening and does not affect the QCP for the bulk. (a)

Discussion and Concluding Remarks
In the present communication, we have started with PAM-a model for a generic TKI. We have used the slave boson approximation (SBA) and then taken the limit U → ∞ to treat the strong interaction between f -electrons in this model. For almost filled f -bands this corresponds to projecting out the states with two or more f -particles per site. This constraint of "no double occupancy", together with the application of a mean-field approximation on the bosonic fields, leads to an effective PAM Hamiltonian that is quadratic in creation and annihilation operators. The quadratic Hamiltonian makes the calculation of the relevant Matsubara propagators extremely simple as one needs to simply solve linear, algebraic simultaneous equations of the Fourier transforms of these propagators and take inverse Fourier transform thereafter. Once this is done, every property of the system (e.g., thermodynamic, transport, optical, thermoelectric, and so on), in principle, could be calculated. In a sequel to this paper, we, in fact, plan to determine the phase diagram of a generic TKI using this approach. Thus, SBA seems to carry great relevance for the PAM. We have also initiated the application of the slave-fermion method to PAM in order to ascertain whether the results obtained here are somehow related to the method used. The initial outcomes are interesting. Furthermore, Equation (28a), which demands that "| t f 1 | should be less than | ε f | by one order of magnitude", is hard to get satisfied unless f -electrons are strongly localized and, therefore, a Dirac-type surface spectrum is possible for a system where NN hopping parameter for f -electrons is smaller than corresponding on-site energy. In fact, the direct evidence [40] of SmB 6 hosting metallic spin helical surface states has been reported recently.
The PAM itself and its extensions are still relevant for the theoretical condensed-matter physics. The examples of extension are those including the on-site repulsion between dand f -electrons given by U d, f d † i.ζ d i,ζ f † i.ζ f i,ζ and thecoulomb interaction V ij f † i.ζ f i,ζ f † j.ζ f j,ζ , between on the nearest-neighbor (NN) sites. For U d, f t d1 , the Kondo singlet state will be destroyed, and the valence fluctuations are enhanced [5]. On the other hand, it was found that, in an investigation [6] of nonmagnetic and magnetic ground states of PAM in the presence of V ij , the magnetic ground state of the PAM in the Kondo regime is unstable for a critical value of V ij . In other words, some V ij critical affects the stability of the magnetic ground state in the Kondo regime in the investigation [6]. We have also introduced an extension to the model here by way of involving the magnetic exchange interaction. The interaction is introduced in the most direct route using only the spin degrees of freedom. We have demonstrated that the exchange field can be used to open a gap at the surface state dispersion. The field, in fact, is expected to play a bigger role, such as in the efficient tuning of the bulk band gap, the plasmon frequency, and so on. Looking backward, we observe there are many unsettled issues. For example, the problem of hybridization of plasmons with optical surface phonons that is likely to occur when TKI is surrounded by a material other than air has not been addressed. We need to suggest the ways and means to curb such a loss.
In conclusion, the inclusion of Rashba impurities in TKI surface gives rise to surface helical states where spin and momentum directions are locked to each other. In the proximity of a superconductor or a magnet, several interesting phenomena, such as the appearance of Majorana anyons, topological Faraday/Kerr effect, fractional charge for the quantum spin Hall effect in 2D, etc., may occur. These offer promising prospects for spintronic applications. Though immense progress has been made in this area over the past decade, never-the-less, we believe that our work may cast new light onto the investigations of how electron correlations and magnetic disorder influence the helical liquid. There are many challenges in the processing of these exotic materials to use the metallic/insulating states in functional devices, and they present great opportunities for the materials science research communities.
Funding: This research received no external funding.

Conflicts of Interest:
The author declare no conflict of interest.