Zitterbewegung-mediated RKKY coupling in topological insulator thin films

The dynamics of itinerant electrons in topological insulator (TI) thin films is investigated using a multi-band decomposition approach. We show that the electron trajectory in the 2D film is anisotropic and confined within a characteristic region. Remarkably, the confinement and anisotropy of the electron trajectory are associated with the topological phase transition of the TI system, which can be controlled by tuning the film thickness and/or applying an in-plane magnetic field. Moreover, persistent electron wavepacket oscillation can be achieved in the TI thin film system at the phase transition point, which may assist in the experimental detection of the jitter motion (Zitterbewegung). The implications of the microscopic picture of electron motion in explaining other transport-related effects, e.g., electron-mediated RKKY coupling in the TI thin film system, are also discussed.

In general, the ZB frequency scales with the energy gap and can be reduced in systems with narrow energy gaps, e.g., narrow gap semiconductors [32] and topological insulators [7]. At the same time, the oscillation of a wavepacket usually decays over time, which results from the interference between oscillations of different momentum-dependent frequencies [2,8]. Therefore, for the ZB effect to be observed, it is crucial to prolong or even indefinitely sustain the oscillatory motion. There have been some proposals to achieve persistent ZB motion, for example, by using semiconductor nanowires [2], or time-dependent systems [33,34]. In principle, we can also design a system in which the ZB oscillation frequency is independent of electron momentum. In this way, we can avoid the interference effect and render the ZB motion persistent and robust against damping.
In this work, we show that such persistent ZB motion can be realized in topological insulator (TI) thin * sonhc85@gmail.com † elembaj@nus.edu.sg films [35][36][37][38]. TI thin films differ from the more commonly studied semi-infinite TI slabs in that they have both a top and bottom surface, each of which can host surface states. The surface states on the two surfaces are coupled to each other due to the finite thickness of the film. In such thin films, the energy gap in the surface states can be controlled by applying an in-plane magnetic field [35] or tuning the thickness of the film [36][37][38]. Topological phase transitions can thus be induced by closing the gap. We show that at the transition point, there exists a momentum-independent oscillation frequency, which can give rise to persistent ZB oscillations of electron wavepackets. Furthermore, we find that the motion of electron in the x − y plane is anisotropic with respect to the injection direction and confined to a certain region of the TI film.The anisotropy of the electron motion due to the ZB effect has consequences for transport-related properties of the thin film system. Here, we focus on the inter-layer interaction between two localized magnetic centers by means of Ruderman-Kittel-Kasuya-Yosida (RKKY) mechanism [39][40][41][42]. The RKKY interaction has been extensively investigated in various systems such as superconductors [43][44][45], topological insulators [46][47][48][49][50][51], Weyl and Dirac semimetals [52][53][54][55][56][57], graphene [58][59][60][61], carbon nanotubes [62,63], semiconductor quantum wires [64,65], and tunneling junctions [66]. The RKKY interaction is mediated by the itinerant electrons. Intuitively, one would then expect the enhancement of the RKKY interaction when the magnetic centers lie along a preferred direction of electron motion, and a corresponding suppression of the RKKY interaction when the electrons are prohibited from moving between the two centers. We find that, indeed, the anisotropy of the RKKY coupling is in line with that of the electron motion. We show that maximum RKKY coupling occurs when the separation between the two magnetic centers is perpendicular to the line connecting the Dirac points.
This manuscript is organized as follows. In section II, we present the model Hamiltonian and derive the dynamics of both plane-wave and wavepacket electrons. We discuss the confinement of the electron trajectory and the regime conditions for persistent ZB oscillation. In section III, the RKKY coupling is calculated in both the weak and strong hybridization limits, and its correlation with the electron motion is also discussed. Finally, section IV contains a summary of our main conclusions.

II. ELECTRON DYNAMICS
We first consider a TI thin film subject to an in-plane magnetic field. For simplicity, we assume that the magnetic field is applied along the x-direction, so that the corresponding gauge field is A B = −ŷBz. As the thickness d of the thin film is comparable to the surface state decay length, the two surfaces are hybridized. The effective Hamiltonian of the system is then [35] where H D (k) = v f (z × σ) · k is the Dirac Hamiltonian describing the topological surface state, in which v f is the Fermi velocity, σ the vector of Pauli spin matrices, and z the unit vector perpendicular to the film (see Fig. 1). ∆ is the hybridization parameter describing the coupling between the top and bottom surfaces, and τ the vector of the Pauli matrices in pseudo-spin space that represents the electron occupancy at the top and bottom surfaces. For simplicity, we set = 1, and introduce the characteristic momenta corresponding to the hybridization energy k ∆ = ∆/v f and the wavevector corresponding to the magnetic field k B = e/2cBdŷ. The eigenenergies of the system are then given by in which s, τ = ± represent the real spin and pseudo-spin indexes respectively, and we define k u = k 2 ∆ + k 2 y , k v = k 2 B + k 2 x , Θ = arctan ky k∆ , and Φ = arctan k B kx . The bandstructure of the TI film is depicted in Fig. (1). An energy gap of Otherwise, the bandstructure is gapless, with the formation of two Dirac cones separated by 2q 0 = 2 k 2 B − k 2 ∆ along the direction perpendicular to the magnetic field.
In particular, at the transition value k ∆ = k B , the two Dirac cones merge to form a single cone. The corresponding eigenstates are given by the fourvectors where N sτ are the normalization factors, and To study the ZB in this multi-band system, we derive the time-evolution of position operator, which is described in the Heisenberg picture asr(t) = e iH0tr (0)e −iH0t , and which at t = 0 is formally represented byr(0) = i∇ k . The time-dependent position operator comprises of a non-oscillatory part that describes the translational motion related to the intraband interference, and an oscillatory part that is associated with the ZB motion [2,8,67,68] and related to the interband interference. Our interest lies in the latter, which is given by [67] (6) in which Ω ij = (E i − E j ), with i = (s, τ ) andÂ ij = iQ i ∇ kQj are the frequency and amplitude of the oscillation, respectively. In the above, we have introduced projection operatorsQ sτ = |ψ sτ ψ sτ |, so that the Hamiltonian (1) can be decomposed as H 0 = sτ E sτQsτ . We can further express the projection operators aŝ whereR andT are involution operators satisfyingR 2 = T 2 = 1, [R,T ] = 0. The explicit forms of these operators are given in Appendix (A). Due to the electron-hole symmetry of the eigenenergies given in Eq. (2), there are only four distinct beat frequencies corresponding to the differences between the energies of interfering eigenstates. These frequencies are given by where k ± is given by Eq. (2).

A. Bound trajectory
Having derived the position operator in Eq. (6), we now trace out the electron trajectory in the system. In general, a free electron can travel in a region as large as the area of the system defined by its physical boundaries, e.g., edges or interfaces. However, we show that in the TI thin film system, the electron trajectory is bound within an area determined by the initial state (spin and momentum) of the electron and the energy gap of the system. Consider an electron injected into the top surface of the TI film with initial spin state in the spin up direction and momentum k which is represented by the planewave |ψ 0 (k) = e ik·r |φ 0 . The position of the electron on the x-y plane at time t can be calculated from Eq. (6) and is given explicitly by where Corresponding expressions for other combinations of injected spin orientation and injection surfaces can be obtained from symmetry arguments. Eq. (1) in terms of k ∆ and k B is, explicitly, Eq. (11) is invariant upon a simultaneous τ reflection about the τ x axis and in-plane spatial inversion, i.e. τ z → −τ z , x → −x, y → −y. This implies that the x and y displacements of electrons injected into the top and bottom surfaces have the same magnitudes but opposite signs. Eq. (11) is also invariant upon a simultaneous spin reflection about σ x (σ y,z → σ y,z ) and reflection along the y axis (x → −x,y → y). This implies that spin up and spin down electrons injected into a given surface (top / bottom) have the same x displacements, and y displacements of the same magnitude but opposite signs. The electron motion of an electron injected in the top surface with initial spin in the +z direction on the x − y plane is depicted in Fig. 2(a) and (b) for different ratios of k ∆ /k B . Taking the initial position of the electron to be the origin, it can be shown that x(t) ≥ 0 for k ∆ /k B < 1, i.e., the electron is always confined in the +x-half of the x-y plane. On the other hand, when k ∆ /k B > 1, the trajectory of the injected electron encompasses the origin as shown in Fig. 2(b).
It can be seen that the electron oscillation comprises both transverse and longitudinal modes. This is a manifestation of the four-band system illustrated in Fig. 1(b), where the quantum dynamics involves not just the evolution of the spin, but also the pseudo-spin degree of freedom, which in our case, represents the surface index (top and bottom surfaces). The electron trajectories in Fig. 2 both the transverse (y) and longitudinal (x) directions. Now in the conventional ZB picture, an electron injected along the x-direction would undergo oscillations in the transverse y-direction, due to the electron spin precession and spin-momentum locking. In this simple picture, the longitudinal oscillations do not seem to play a role. To explain the emergence of the longitudinal oscillations, we need to consider the pseudo-spin (τ z ) degree of freedom. This can be ascribed to the precession of the pseudo-spin, which represents the back and forth tunneling between surfaces. From Eq. (1), this pseudo-spin dynamics is coupled to the longitudinal motion. Indeed, as shown in Fig. 2(d), the electron lies in the positive x-half when it is on the top surface, and would move to the negative x-half after tunneling to the bottom surface. Thus, the back and forth tunneling between the surfaces mediated by the hybridization ∆ translates into the oscillation of the electron motion in the longitudinal x−direction. In the thick TI film limit where the top and bottom surfaces are decoupled, i.e., ∆ = 0 in the Hamiltonian Eq. (1), the motion of the electron is simply given by r(t) = 1 k 2 (z ×k) 1 − 2 cos 2 kv f t . Surprisingly, a spin-up electron initially injected along the x-direction will only move in the y-direction, i.e, its trajectory is confined in a line perpendicular to the injection direction. This can be explained by considering the electron velocity given by v = ∂ k H = v f (z × σ). The electron spin precesses as is perpendicular to the momentum.

B. Wavepacket dynamics
In the previous section, we have considered the trajectory of a single electron. We now consider the more practical case of an electron wavepacket, which is a superposition of different momentum states. In general, the beat frequencies ws as given in Eq. (8) are dependent on the momentum. Thus, when evaluating the expectation value of the position operator for a wavepacket, the resulting interference of oscillations with different momentum-dependent frequencies would, in general, lead to a decay of the ZB over time. In order to sustain the ZB motion, we need to realize a scenario where at least one beat frequency is momentum-independent. We will show that such a scenario can be achieved by the appropriate choice of parameters such as hybridization energy and the in-plane magnetic field.
Suppose that the electron is injected in the y-direction, i.e., k x = 0, k y = k, and the hybridization and magnetic field are tuned so that k ∆ = k B . In this case, the beat frequencies of Eq. (8) are now given by (12) in which we recall that k u = k 2 ∆ + k 2 y . We can see that besides the three momentumdependent frequencies, there is one frequency w 2 that is independent of momentum. At large time scales, we would expect the oscillations associated with the other three frequencies to decay away due to interference, while the oscillations associated with the k-independent w 2 frequency would persist. This is one of the main results of this paper.
To quantitatively verify the above intuitive picture of persistent ZB motion, we consider the electron wavepacket given by where |φ 0 is the initial spin state, and a(k) = 1 √ πδk e − (k−k 0 ) 2 2δk 2 is the Gaussian distribution function that represents the spread of the electron state in momentum space, in which k 0 and δk are the initial momentum and line-width, respectively. The expectation value of the position operator Eq. (6) for the above state is given by where the integration is taken over momentum space. As a consequence of the wavepacket spread in k-space, the ZB will generally decay over time. In order to analytically describe the damping process, we will consider the narrow wavepacket limit, i.e., δk/k 0 1, so that the integration of the Gaussian function in Eq. (14) can be approximated by up to O(δk 4 ). In the above, the first term is the initial ZB oscillation with momentum k 0 , and the second term represents the deviation of the ZB around the packet center. Substituting the position operator in Eq. (6), we have in which a = x, y. The first term in the above describes oscillations with constant amplitude that are in-phase with the initial oscillation. The next two terms have time-dependent amplitudes that are linear and quadratic in time, respectively. Rearranging the equation (15), the ZB of a wavepacket can be expressed as where the decay time is defined as with the beat frequencies given by Eq. (8). At short t, the first term in Eq. (17) can be formally written as r Z (t) ≈ ij r ij (k 0 , t)e −t 2 /T 2 ij , which expresses the exponential decay of the ZB (see Fig. 3).
From Eq. (12), the decay times are obtained as where θ(x) is the Heaviside step function, T 2 corresponds to the combinations of i and j where |Ω ij | = w 2 , and T d the other combinations of i and j. As can be seen, when one of the beat frequencies, i.e., w 2 , becomes independent of momentum at resonance where k ∆ = k B , the associated decay time T 2 in Eq. (19) goes to infinity. This implies that the ZB related to this mode will be persistent. In this case, the steady state transverse oscillation is given by In the limit of large hybridization k ∆ k 0 , the persistent oscillation reduces to y(t) ≈ A ZB cos(ω ZB t) with A ZB = 1 2k∆ and ω ZB = 2v f k ∆ = 2∆ being respectively the amplitude and frequency. This persistent oscillation is depicted by the orange line in Fig. 3(b). Surprisingly, both the amplitude and frequency of the persistent mode do not depend on the initial momentum and width of the injected wavepacket and are instead determined by a single parameter, i.e., the hybridization energy. Following Eq. (19) , the ZB has a sharp transition from a transient to persistent mode at k ∆ = k B at which the bulk gap closes (Fig. 1) and the TI film undergoes a topological phase transition [35]. We can hence refer to the persistent oscillation as a topological mode of electron oscillation.

III. ELECTRON-MEDIATED RKKY INTERACTION
In the previous section, we have shown that the electron trajectory is confined and may be highly anisotropic (see e.g., Fig. 2(a)). This has consequences for the transport-related properties of the system, such as the electron-mediated RKKY interaction. The confinement of the electron trajectory implies that the electrons are not able to mediate information, e.g. angular momentum, between magnetic moments separated by a separation distance that exceeds the confinement region. In order to verify this effect, we consider two magnetic centers S i (i = 1, 2) located at R i . The electron-mediated exchange interaction between the magnetic centers is modeled by where J is the exchange coupling. The exchange interaction can be considered as a perturbation to the Hamiltonian in Eq. (1). For simplicity, we assume that R 1 = (0, 0), and R 2 = R(cos φ R , sin φ R ). We show that the RKKY coupling between the two magnetic centers does not depend on just the distance R, but also on the direction φ R between them.
In the framework of the second-order perturbation theory, the effective interaction between two magnetic im-purities is given by [39-42, 46, 64, 66] where + = + i0 + , Tr stands for the trace over the spin degree of freedom, and the expanded spin operator in spin and pseudo-spin spaces is defined asσ = τ 0 ⊗ σ, in which τ 0 is the identity matrix of rank 2. The Green's function in real space is given by the Fourier transformation where G(q, + ) = [ + − H 0 (q)] −1 is the Green's function in momentum space, and A BZ is the area of the first Brillouin zone. Let us first consider the weak hybridization limit, i.e., k ∆ k B . In this limit, the system is gapless and the two Dirac points are separated by q 0 ≈ k B . The analytical expression of the RKKY coupling can be obtained as (see Appendix B for more details) in which the range functions are In the above,û = (R × z), withR = R/R being the unit vector along R. The RKKY coupling in Eq. (25) consists of three terms: the Heisenberg exchange, the spin-frustrated, and the Dzyaloshinsky-Moriya interaction terms. As shown above, the RKKY coupling exhibits not only the usual R −2 distance dependence [46] in a semi-infinite thick TI slab with only a single surface, but also has an additional direction-dependence due to the cos(2q 0 · R) factor that is absent in the semi-infinite thick slab. This directional dependence stems from the contribution of the surface states on both surfaces of the film in mediating the effective exchange coupling, and the fact that the corresponding Dirac cones are separated in momentum space. In the case where the two magnetic impurities are separated along the x-direction, i.e., along the in-plane magnetic field direction, k B · R = 0 and the RKKY coupling reaches its maximum. This can be explained by considering the process of indirect exchange coupling between the two magnetic moments via the itinerant electrons. When an electron is in close proximity to the first magnetic moment, its spin angular momentum is coupled to that of the magnetic moment. If there is finite electron overlap with the second magnetic moment, then its spin angular momentum is also coupled to the second moment. In this way, an effective exchange coupling arises between the two magnetic moments. The strength of the effective coupling depends on the rate and probability of electron overlap between one magnetic moment and the other. In other words, if the second magnetic moment is located at a position with little electron overlap with the first magnetic moment, then the coupling between the moments would be weak. Conversely, if the second magnetic moment is at a position where the electron has a high probability of overlap, the coupling will be enhanced. In our case, when the magnetic field is applied along the x-direction, the electron motion has a tendency of being confined along the same x-direction [see Fig. 2(a)]. This means that a second magnetic moment placed along the x-direction with respect to the first moment will have high probability of being coupled by an intermediary electron, thus inducing stronger RKKY coupling.
Although Fig. 1 shows only the results of a spin up electron injected on the top surface explicitly, the results of the symmetry analysis following Eq. (11) imply that the electron trajectory will still be confined along the x direction for spins of other orientations injected into both the top and the bottom surfaces.
To quantify the correlation between the RKKY coupling and the electron trajectory, we will analyze the preferred direction of the electron motion. As the electron position oscillates over time as described in Eq. (9), we consider its average valuer(k) = lim T →∞ 1 T T 0 r(t), which is explicitly given bȳ In the limit of weak hybridization energy k ∆ k B , the above reduce tor Eq. (30) is the time-averaged position of an electron with momentum k. Now, averaging the above over momentum space up to the Fermi wave-vector, we obtain which indicates that the electron will preferably move in the direction perpendicular to the direction separating the two Dirac cones. Therefore, when R is parallel tor, and thus perpendicular to k B , the RKKY coupling strength will be maximum. This is in line with the prediction based on the electron trajectory, as discussed above.
We note that in the above, the preferred motion direction was obtained based on the position of the plane-wave electron. Here, we show that the preferred direction is the same if we consider the electron wavepacket treatment. As the Gaussian function in the wavepacket picture is time-independent, it would not alter the position value after time-averaging. From Eq. (14), the average position of an electron wavepacket initially centered at k 0 is simply derived asr pk (k 0 ) = k |a(k − k 0 )| 2r (k) ≈ r(k 0 )+δr(k 0 ), where the deviation δr(k 0 ) = δk 2 4 ∇ 2 kr (k 0 ) follows Eq. (15) for a narrow wavepacket, withr(k) given in Eq. (29). In the weak hybridization limit, applying Eq. (30), we find that the deviation δr(k 0 ) = 0, which means that the preferred direction of motion of a wavepacket coincides with that of a plane-wave electron. This result thus suggests that one may use the wavepacket treatment in understanding properties of the RKKY coupling, besides the conventional plane Bloch wave approaches [39][40][41] in future works.
On the other hand, in the strong hybridization limit k ∆ k B , the RKKY coupling is given by (details are shown in Appendix B) where the range functions are now given bỹ In this case, the surface states are gapped (see Fig.  1(d)), and the Dirac cones vanish. In this limit, the RKKY coupling becomes isotropic, i.e., it is independent of the angle between the magnetic centers. This result is consistent with the calculated electron trajectory in the gapped scenario, where its trajectory is almost isotropic in the 2D plane (see Fig. 2(c)). This can be further verified by considering the time-averaged electron position as outlined above, which is given bȳ The above goes to zero upon averaging over momentum space, so that there is no preferred direction of the electron motion in the 2D plane in this case. We remark here that in the insulating phase, the TI film has been shown to have a diamagnetic response to an in-plane magnetic field [35]. As a consequence, the magnetic moments in the TI film may acquire an additional magnetic response and the steady state magnetization may change accordingly. However, the magnetic susceptibility is extremely small, i.e. on the order of 10 −8 [35], which is several orders of magnitude smaller than even the small diamagnetic susceptibility of typical metals. The effect of the induced magnetization can thus be neglected in the bandstructure of the TI film. Since the RKKY coupling is derived from the bandstructure of the TI film, therefore it will not susceptible to this diamagnetic response.
In addition, we note that if the Fermi level lies within the gap in the insulating phase, the RKKY mechanism is no longer be valid as it relies on itinerant electrons. Instead, the indirect exchange coupling is now described by the van Vleck mechanism as discussed in previous works [69][70][71]. In our work, we assume that the Fermi level is finite, i.e., within the conduction band, and ignore the van Vleck coupling for simplicity.

IV. CONCLUSION
In this paper, we investigated the anomalous motion of electrons in topological insulator thin films. First, we showed that due to the hybridization of the surface states with opposite helicities, a spin-polarized electron will undergo oscillatory motion within a confined region. Furthermore, the oscillation is anisotropic with the preferred direction being along the separation of the two Dirac points, a finding which that be ascribed to the anisotropy of the Fermi circle. As a consequence, the direction and distance dependence of RKKY interactions mediated by itinerant electrons between two magnetic impurities in thin TI films have a strong correlation with the electron motion. Interestingly, it was found that the RKKY coupling is maximized when two impurities at a fixed distance are positioned along the separation direction of the two Dirac points. This finding is consistent with the preferred direction of the confined electron motion.