Long-range spin transport on the surface of topological Dirac semimetal

We theoretically propose the long-range spin transport mediated by the gapless surface states of topological Dirac semimetal (TDSM). Low-dissipation spin current is a building block of next-generation spintronics devices. While conduction electrons in metals and spin waves in ferromagnetic insulators (FMIs) are the major carriers of spin current, their propagation length is inevitably limited due to the Joule heating or the Gilbert damping. In order to suppress dissipation and realize long-range spin transport, we here make use of the spin-helical surface states of TDSMs, such as $\mathrm{Cd_3 As_2}$ and $\mathrm{Na_3 Bi}$, which are robust against disorder. Based on a junction of two FMIs connected by a TDSM, we demonstrate that the magnetization dynamics in one FMI induces a spin current on the TDSM surface flowing to the other FMI. By both the analytical transport theory on the surface and the numerical simulation of real-time evolution in the bulk, we find that the induced spin current takes a universal semi-quantized value that is insensitive to the microscopic coupling structure between the FMI and the TDSM. We show that this surface spin current is robust against disorder over a long range, which indicates that the TDSM surface serves as a promising system for realizing spintronics devices.


I. INTRODUCTION
Transmission of signals over a long distance is essential in designing integrated information devices. While charge current in normal metals is inevitably subject to the Joule heating, spin current is now intensely studied to realize a long-range transmission with less dissipation, in the context of spintronics [1][2][3]. Spin current, namely the flow of spin angular momentum, is carried by various types of quasiparticle excitations in materials. In metals, spin current is carried by conduction electrons [4]. Electron spin current can be generated by current injection from magnetic metals [5], spin pumping from magnetic materials with magnetization dynamics [6][7][8], the spin Hall effect [9][10][11][12][13], etc.
Spin waves (or magnons) in magnetic materials, namely the dynamics of the ferromagnetic or antiferromagnetic order parameters, are also elementary excitations that carry spin and heat currents [14][15][16][17][18]. Magnon spin current can be generated by the magnetic resonance under a microwave [19,20], by the spin Seebeck effect under a temperature gradient [21,22], etc.
While those spin-current carriers are available in various materials commonly used in experiments, their propagation length is inevitably limited due to dissipation. Conduction electrons are subject to scattering by phonons and disorder, which results in the Joule heating. Magnon spin current in insulators is considered to be advantageous in that it is free from the Joule heating [23][24][25][26], whereas the Gilbert damping of spins leads to the dissipation of spin and energy. Due to such dissipation effects, transmission of spin current over a long range is a challenging problem. Recent theoretical and experimental studies showed that phonons in nonmagnetic insulators may realize long-range transport of spin current [27][28][29][30][31][32]: since circularly-polarized transverse phonons carry angular momentum, they can mediate spin current between magnets via magnetoelastic coupling. However, for the efficient interconversion of magnons and phonons, one needs fine tuning of the magnon frequency. Such limitations in long-range spin transport restrict the de- There are four spin-polarized channels on the surface of the TDSM connecting FM1 and FM2. The magnetization dynamics induces population imbalance among the four edge channels, which yields a net spin current from FM1 to FM2.
sign of highly integrated spintronics devices.
In the present work, we theoretically propose a longrange transmission of spin mediated by the surface electronic states of topological Dirac semimetals (TDSMs). TDSM is a class of three-dimensional (3D) crystalline materials having pair(s) of Dirac nodes in the electronic band structure in the bulk [33][34][35]. The TDSM phase is experimentally realized in Na 3 Bi [36,37] and Cd 3 As 2 [38][39][40][41][42]. TDSMs have quasi-1D gapless states on the surface, which arise as Fermi arcs connecting the Dirac points projected onto the surface Brillouin zone [44,45]. These surface states are spin helical, i.e. spin-↑ and spin-↓ states propagate along the surface oppositely to each other, and are protected by the Z 2 topology in the bulk [46,47]. These features are analogous with the helical edge states of 2D quantum spin Hall insulators (QSHIs) [48][49][50], and hence the surface states of TDSMs are robust against disorder as long as the system preserves time-reversal symmetry [51,52]. From these features, we can expect that the helical surface states of TDSMs are suitable for long-range spin transmission.
Indeed, in the context of 2D QSHIs, it was theoretically proposed in numerous literatures that the helical edge states are capable of interconversion between spin angular momentum and electric current, based on the topological field theory, the numerical simulations, the scattering theory, and the Floquet theory [53][54][55][56][57][58][59][60]. The electric current and the spin torque arising from the interconversion take quantized values irrespective of the microscopic coupling structure between the edge electrons and the spins in magnets, which is traced back to the band topology of the 1D edge states. The similar spincharge interconversion is expected also on the helical surface states of TDSMs, as long as the Fermi level is tuned in the vicinity of the Dirac points so that the bulk transport may be negligible [61]. However, the spin-charge conversion discussed in those works occurs locally at the interface of a magnet and a TDSM (or a QSHI), and a theory for nonlocal transmission of spins with the helical edge states over a long distance, which is essential for device application, is not well established.
Based on the above background, we here consider the transmission of spin angular momentum between two ferromagnetic insulators (FMIs) connected by a TDSM, as schematically shown in Fig. 1 (a). We assume a magnetization dynamics in one of the FMIs (FM1), and discuss how the spin current transmitted through the TDSM exerts a spin torque on the other FMI (FM2), by constructing analytical and numerical schemes to evaluate the nonlocal spin transmission between the two FMIs separated at a distance. As a result, we find that the transmitted spin current takes a semi-quantized value, which is determined only by the configuration of Dirac points in momentum space and the frequency of magnetization dynamics. This semi-quantization of spin current can be understood analytically as the electron transport on the helical surface states, which comes from the imbalance of electron population among the edge channels driven by the magnetization dynamics [see Fig. 1(b)]. Moreover, from the numerical simulation of the real-time dynamics of electrons in the whole 3D system, we directly confirm that this semi-quantized spin transmission is robust against moderate disorder in the bulk. These results imply that the TDSM surface may serve as a promising system for highly-integrated spintronics devices with a long-range spin transmission.
This article is organized as follows. In Section II, we review the generic characteristics of TDSMs, and give a detailed explanation about our model setup with a TDSM and FMIs shown in Fig. 1. In Section III, we give analytical expressions of the flow of electrons and spin on the surface, based on the 1D scattering theory. (The detailed calculation processes are shown in Appendices.) In Section IV, we show the results of our numerical simulations within the whole 3D system, and discuss their consistency with the analytical expressions and their robustness against disorder. Finally, in Section V, we give some experimental implications from our calculations and conclude our discussion. Throughout this article, we take the natural unit = 1.

II. MODEL SETUP WITH TDSM
In order to demonstrate spin transmission through a TDSM, we construct a model that we shall use both for the analysis and for the numerical simulation, as shown in Fig. 1. We first review the generic characteristics of TDSMs, and give a detailed explanation about our model setup, with a junction of a TDSM and two ferromagnetic insulators (FMIs), on the basis of those characteristics.
TDSMs are characterized by a pair of Dirac points in momentum space, which are protected by rotational symmetry around a crystal axis [33,46,47]. If we take this axis as z-axis, the Dirac points are located on k z -axis, which we denote as k ± D = (0, 0, ±k D ). Due to the rotational symmetry, the spin component σ z serves as a good quantum number around k z -axis, which means that the spin-↑ and spin-↓ states are degenerate around the Dirac points. The minimal model for such a band structure at low energy is composed of four degrees of freedom, with twofold spins and twofold orbitals [34,36,38], Here the Pauli matrices σ x,y,z act on the spin subspace and τ x,y,z on the orbital subspace. In the effective model of Cd 3 As 2 [38], for instance, the basis functions with τ z = + correspond to the 5s-orbitals of Cd with the total angular momentum J z = ±1/2, while those with τ z = + correspond to the 4p-orbitals of As with J z = ±3/2, and σ z = ± represents the sign of J z . The parameter m 0 characterizes the band inversion, which leads to the Dirac points of k D = m 0 /m 1 if m 0 , m 1 > 0.
Due to the band inversion from spin-orbit coupling, the system shows the intrinsic spin Hall effect. By fixing k z in the band-inverted region −k D < k z < k D and considering the 2D slice, as shown in Fig. 2, the Hamiltonian reduces to that for the 2D quantum spin Hall insulator (QSHI), with the quantized spin Hall conductivity σ s(2D) xy = e/2π. Therefore, by multiplying the number of the 2D slices ν z = 2k D /2π per unit length in z-direction, the spin Hall conductivity of the TDSM in 3D takes the "semi-quantized" value if the Fermi level is in the vicinity of the Dirac points [34,35]. Another consequence of the band inversion is the emergence of surface states. On the surfaces parallel to the rotational axis (z-axis), which we call the side surfaces, there emerge spin-helical states, with the spin-↑ and spin-↓ modes propagating along the surface oppositely to each other [44]. These surface states can be regarded as the collection of 1D helical edge states of the 2D QSHI at fixed k z . They form a pair of Fermi arcs connecting the Dirac points projected onto the surface Brillouin zone, which are robust against disorder as long as time-reversal symmetry is preserved [52]. The contribution of these surface states to the electron transport was observed experimentally, as the quantum oscillation under a magnetic field [40,51,[62][63][64][65][66][67].
In order to make use of the helical surface states for spin transmission, here we consider the model setup as shown in Fig apart by the distance L x , and each of them is attached to the TDSM by the length L y . We assume a situation where the magnetization of FM1 is steadily precessing around z-axis, which is maintained by providing angular momentum and energy externally by microwaves, etc. Under such a setup, we estimate the spin torque exerted on FM2, which corresponds to the spin current transmitted from FM1 to FM2 via the TDSM, both analytically and numerically.

III. TRANSPORT ANALYSIS ON THE SURFACE
In this section, we treat the spin transmission through the TDSM analytically, by focusing on the spin transport mediated by the helical surface states on the surface. If the Fermi level is in the vicinity of the Dirac points, the bulk transport becomes negligible and the surface transport becomes dominant. In order to evaluate the spin transmission between two FMIs phenomenologically, we first formulate the transmission of charge and spin at a single interface with a FMI. By using this single-FMI picture as a building block, we formulate the spin transmission between the two FMIs in our model setup shown in Fig. 1.
As mentioned in the previous section, we regard the helical surface states of the TDSM as the collection of 1D helical edge states, which reside at every k z in the band-inverted region (−k D < k z < k D ). As long as the translational symmetry in z-direction is satisfied, k z serves as a good quantum number, and the contribution from 1D helical edge states at each k z can be treated separately. Therefore, in this section, we first consider the spin transmission by a single pair of 1D helical edge states, and then multiply the number of 2D slices ν z = 2k D /2π per unit length in z-direction, to evaluate the overall contribution from the 2D surface states.

A. Charge and spin pumping by a single FMI
In a manner similar to the theoretical treatment of spin pumping and injection by a ferromagnet [7,8], we formulate the role a FMI coupled to the helical edge as a timedependent scatterer. We consider the scattering process in the hypothetical 1D space along the edge, where we denote its 1D coordinate as x (see Fig. 3).
As in the conventional Landauer-Büttiker formalism in mesoscopic systems [68][69][70][71][72], the charge and spin currents can be derived by comparing the numbers of the incoming and outgoing electrons for the scatterer. Since the nonmagnetic edge states are spin-helical, there are two incoming channels and two outgoing channels: as the incoming channels, spin-↑ electrons come from the left (x < 0) and spin-↓ electrons come from the right (x > 0). As the outgoing channels, spin-↑ electrons go out to the right (x > 0) and spin-↓ electrons go out to The pumping process is mapped to the scattering problem in the quasi-1D space (x ) along the edge, by regarding the interface region with the precessing magnetization n(t) as the time-dependent scatterer. By solving the scattering problem as described in Section III B, we see that the outgoing channel with spin-↑ becomes more populated than that with spin-↓.
the left (x < 0). We denote the annihilation (creation) operators of an electron with energy in these channels as a ↑/↓ ( ) for the outgoing channels, respectively.
The electron distributions in these channels with energy are given by taking the quantum average · of the operators defined above [69][70][71], which we use throughout this section as the main tool to evaluate the transmission of spin current. With these distribution functions, the numbers of electrons coming into / going out of the scatterer with spin-↑ / ↓ per unit time are given by By using these notations, the charge current flowing from the left to the right is given by The two formalisms are equivalent due to the charge conservation at the scatterer. On the other hand, spin can be transferred between the electrons and the FMI, and thus the net spin current flowing out of the scatterer can be nonzero. Noting that each electron carries spin ±1/2, the spin current pumped by the FMI, namely the net spin angular momentum flowing into and out of the FMI per unit time, is given by From these relations, we can immediately see a simple relation between the charge and spin currents, where the right-hand side is determined only by the numbers of incoming particles and is independent of the scattering process. In particular, if the numbers of spin-↑ and spin-↓ electrons entering the scattering region are equal, its right-hand side vanishes and reduces to the simple relation, I s = −I/e. In order to evaluate the charge current I and the spin current I s separately, we need relations between the incoming and outgoing distributions that are determined by the scatterer. If the magnetization n is periodically precessing as n(t) = (sin θ cos(Ωt + φ), sin θ sin(Ωt + φ), cos θ), (8) where Ω is the precession frequency and θ is the azimuthal angle, the energy of electron is not conserved in the scattering process. Such a time-dependent scattering problem can be solved by taking the "rotating frame" of spin: by the time-dependent unitary transformation on the edge electrons, which rotates their spin by the angle Ω per unit time around z-axis, the magnetization direction is fixed to n 0 ≡ n(t = 0) in this rotating frame [57]. Since this transformation U (t) shifts the energies of spin-↑ / ↓ electrons by ±Ω/2, respectively, the operators in the rotating frame, which we denote byã ↑/↓ andb ↑/↓ , are related to those in the rest frame as The operators for the incoming and outgoing channels are related by the S-matrix. By using the S-matrix in the rotating framẽ which can be obtained by solving the time-independent scattering problem with the fixed magnetization (see Appendix A for detail), the operatorsã ↑/↓ andb ↑/↓ in the rotating frame are related as [69][70][71] b Note that the components in the S-matrix satisfy the reversibility relations and the unitarity condition R( ) + T ( ) = 1 (15) due to the time-independence of the scatterer in the rotating frame, where R( ) is the reflection rate and T ( ) is the transmission rate.
With the S-matrix defined above, we are ready to evaluate the electron distributions in the outgoing channels. By substituting Eqs. (10) and (12) to Eq. (3), we obtain the relations for the distribution functions in the rest frame as which we shall use as the fundamental relations throughout this section to evaluate the spin transmission. By integrating over the energy and using the unitarity relation R( ) + T ( ) = 1, we obtain the relations between the numbers of incoming and outgoing electrons,

B. Scattering rates and quantized pumping
The scattering rates R( ) and T ( ) are defined in the rotating frame of spin, by a FMI with its magnetization fixed to the direction n 0 . The Hamiltonian for the edge electrons coupled to this magnetization in the rotating frame is given as where we define the 1D momentum along the edge as k . We can immediately see from this matrix form that the in-plane component of the magnetization opens an exchange gap J ⊥ ≡ |J sin θ| around zero energy in the edge spectrum, which influences the scattering of the edge electrons by the FMI as we shall see in the following discussions. By evaluating the S-matrix in the rotating frame, whose detailed derivation process is shown in Appendix A, the scattering rates are given as with K = − J 2 ⊥ /v edge the wave number inside the interface region coupled with the FMI. (Note that these forms are valid for inside the exchange gap, | | < J ⊥ , as well, where K becomes pure imaginary and the wave function inside the interface region exponentially decays by x .) The numerical behavior of R( ) is shown in Fig. 4, by varying the electron energy and the length of the interface region L.
The most important feature in the reflection rate R( ) is that it reaches unity for | | < J ⊥ , which means that the electron inside the exchange gap is totally reflected, if the length L of the magnetic region is long enough. This is because the electron wave function in the magnetic region, at energies inside the exchange gap, decays exponentially. The decay length of the wave function at = 0 is given as and hence the tunneling through the magnetic region is fully suppressed if L l 0 . On the other hand, if is out of the exchange gap, the reflection rate R( ) oscillates as a function of due to the formation of resonance states inside the interface region.
By using the scattering rates obtained above, we can evaluate the charge and spin currents pumped by the FMI. If we assume that the distributions of the incoming channels f  18) and (19) yield the balance of incoming and outgoing electron numbers, Schematic picture of the model setup shown in Fig. 1 sliced at a fixed kz. On this 2D slice, FM1 and FM2 are connected by four edge channels, two with spin-↑ and two with spin-↓. We consider the electron distributions in these channels, to understand the flow of spin I s 1 and I s 2 .
up to the first order in Ω for slow magnetization dynamics. This relation means that the number rate of outgoing electrons with spin-↑ is raised and that with spin-↓ is lowered by R( F )Ω/2π due to the magnetization dynamics, as schematically shown in Fig. 3.
In particular, if the Fermi level F is inside the exchange gap (| F | < J ⊥ ), the electrons at the Fermi level are fully reflected, i.e. R( F ) ≈ 1. The right-hand side of Eq. (24) thus reduces to Ω/2π, which corresponds to one electron per a precession period T p = 2π/Ω. As a consequence, the electric current pumped through the magnetic region becomes quantized as which is consistent with the previous literatures on the helical edge states of QSHI [53][54][55][56][57][58][59][60]. The spin injection rate (per unit time) from the FMI into the edge electrons is also quantized asĪ which satisfies the relation in Eq. (7). The overall contribution from the 2D surface states of TDSM is given by multiplying the factor ν z = 2k D /2π to those quantized values.

C. Spin transfer between two FMIs
We now consider the model setup constructed in Section II. By slicing the system at fixed k z in the bandinverted region (−k D < k z < k D ), FM1 and FM2 are connected by four channels, namely two pairs of counterpropagating modes with spin-↑ and ↓, as shown in 2↑ ( ) after a time T x = L x /v edge , and in similar manners for the other channels. Electron propagation on these channels leads to spin transmission between FM1 and FM2. We do not take into account the contribution from the bulk electrons here, which is a valid approximation if the Fermi level is set in the vicinity of the Dirac points so that the density of states in the bulk is small enough. Now we consider the spin transmission, with the magnetization in FM1 precessing around z-axis by the frequency Ω and that in FM2 fixed in x-direction. Here, the transmission and reflection coefficients at FM1 are the same as those obtained in the previous subsection (with L → L y ), and those at FM2 are given by setting Ω = 0. Therefore, in a manner similar to Eqs. (16) and (17), the incoming and outgoing distribution functions are related as follows: For simplicity of discussion, we assume here that electron transmission and reflection at each magnetic region occur instantaneously, which is satisfied if L y v edge T p .
Based on the above relations, we evaluate the flow of charge and spin between FM1 and FM2, driven by the magnetization dynamics in FM1. As the initial condition, we start with the system in equilibrium, where all the edge channels are in the equilibrium distribution f 0 ( ) filled up to the Fermi level F the Fermi level. We then switch on the magnetization dynamics in FM1 at time t = 0 adiabatically, so that the switch-on process may not cause any significant disturbance in the electron distributions, and estimate the transient behavior of the edge electrons after the switch-on by considering the following steps (a)-(d). (The schematic pictures corresponding to these steps are shown in Fig. 6.) (a) Before the magnetization dynamics is switched on (t < 0), all the edge channels are in the equilibrium distribution f 0 ( ).  are given as In particular, if the Fermi level F is inside the exchange gap and the magnetization dynamics is adiabatic (Ω J ⊥ ), we can apply the same discussion as in Eq. (24) in the previous subsection. While the number rates of the incoming electrons per unit time I gets lowered by Ω/2π, due to the magnetization dynamics in FM1 (see Fig. 6(b)). Therefore, both the electric current passing through the magnetic region and the spin current pumped from FM1 reach the quantized values, per a single 2D slice at k z .
(c) For time T x t 2T x , the electrons going out from FM1 at the step (b) reach FM2, and are reflected or transmitted by FM2. As a consequence, the outgoing distributions for FM1 given by Eqs. (31) and (32) is transferred from the edge electrons to FM2. By substituting Eqs. (29)-(32), this spin current reads which is the general form applicable to arbitrary precession frequency Ω and equilibrium Fermi energy F . In particular, if the precession frequency Ω in FM1 and the Fermi level F are inside the exchange gap J ⊥ (i.e. Ω, F J ⊥ ), we can again derive the quantized pumping. If the interface region is sufficiently long (i.e. L y l 0 ), the incoming electrons in the vicinity of F are fully reflected at FM2. The incoming electron with spin-↑, whose number rate per unit time is raised by Ω/2π, flips its spin on the reflection process at FM2 and injects spin 1 to FM2 for each electron, whereas that with spin-↓ is lowered by Ω/2π and injects spin −1 to FM2 (see Fig. 6(c)). Therefore, the spin current I s 2 transmitted to FM2 reaches the universal valuē which means that spin 2 is injected into FM2 during a precession period T p per a single 2D slice at k z . For the 3D TDSM, the injected spin current (per unit length in z-direction) takes the semi-quantized valuē This (semi-)quantization of spin current is one of the main results in our analysis, which is determined only by the number of helical channels on the surface and is independent of the microscopic structures in the bulk. Note that this relation is satisfied only if the electrons at the Fermi surface are fully reflected by FM2. If Ω or F is out of the exchange gap, or once magnetization dynamics is driven in FM2, some electrons are transmitted through the magnetic region of FM2, and I s 2 becomes not exactly twice of I s 1 . (d) For time 2T x t, the electrons reflected at FM2 at the step (c) reaches FM1. At this step, all the edge channels are no longer in equilibrium distribution. Moreover, once the magnetization in FM2 acquires dynamics, the incoming electrons at FM2 are no longer fully reflected. Therefore, the spin currents I s 1 and I s 2 deviate from the (semi-)quantized valuesĪ s 1 andĪ s 2 at this stage. In order to make use of the (semi-)quantization of spin current, the magnetization dynamics in FM1 should be in a pulse shorter than the time scale T x .
The important feature seen from the analysis above is that the value of the spin current transmitted by the surface states of TDSM (during the steps (b) and (c)) is universal, which is determined only by the location of the Dirac points k D and the precession frequency Ω, as given by Eq. (37). It only requires the existence of a sizable exchange gap in the helical surface states induced by the FMIs, and is insensitive to the microscopic structure and value of the exchange coupling. Moreover, the transmitted spin current is independent of the system size L x,y , since only a single pair of helical edge states contribute to the spin transmission for each 2D slice at k z . This analytical estimation is valid as long as the surface states are robustly present against disorder, which shall be checked by the numerical simulation in the next section.

IV. NUMERICAL SIMULATION ON LATTICE
In this section, we present our numerical simulation of the spin transmission process via a TDSM, which is performed with the 3D lattice model of a TDSM. For the numerical simulation, we use the model constructed in Section II, with two FMIs (FM1/2) connected by a TDSM. By following the real-time evolutions of the wave function of all the electrons in the TDSM and of the magnetization in FM2, we evaluate the flow of spin driven by the magnetization dynamics in FM1. As a result, we find that the transmitted spin current reaches the semiquantized value at the early stage after switching on the magnetization dynamics. This semi-quantization behavior of the spin current agrees with the surface transport picture employed in the previous section, which implies that the spin transmission in the TDSM is dominated by the surface states. Moreover, we observe that this semiquantized spin transport is robust under disorder even at a long range.

A. Model
For the numerical simulation, we use a lattice model of a TDSM [36,38]. On a hypothetical cubic lattice with the lattice spacing a, the tight-binding Hamiltonian reproduces the low-energy effective model in Eq. (1) around k = 0, with the correspondence of parameters v = au, m 0 = r 0 , m 1 = a 2 r 1 /2. This lattice Hamiltonian gives a pair of Dirac points located at k ± D = (0, 0, ±k D ), with k D = a −1 arccos(1 − r 0 /r 1 ). Throughout our calculation, we fix the parameters r 0 = r 1 = u, which gives k D = π/2a.

In real space, the Hamiltonian becomes
in the operator formalism, where c ( †) r is the fourcomponent annihilation (creation) operator at the lattice cite r. The vectors a i=0,x,y,z are defined as a 0 = 0 and a x,y,z = ae x,y,z , with the Cartesian unit vectors e x,y,z , and the matrices D i=0,x,y,z are defined as In order to simulate the setup shown in Fig. 1(a), we here take open-boundary conditions in x-and y-directions, with the size represented by L x and L y . For z-direction, we take a periodic boundary condition, with the size L z . The number of sites N x,y,z in each direction is related to the system size by L x,y,z = aN x,y,z . The magnetizations in the FMIs are defined as macrospins, with their directions denoted by the unit vectors We investigate their dynamics by solving the Landau-Lifshitz-Gilbert (LLG) equation, as we discuss in detail in the next subsection, and hence we do not implement their dynamical properties in the Hamiltonian of the TDSM. We require that the magnetizations n 1,2 are coupled to the electron spins in the TDSM at the boundaries x = 0 and x = L x , respectively. The coupling is described by the Hamiltonian with the phenomenological coupling constant J exc . The matrix Σ characterizes how the exchange coupling depends on the atomic orbital (s or p) that each electron in the TDSM belongs to [73]. Here we define it as Σ = (1 + τ z )σ, so that the structure of the exchange coupling shall be invariant under a C 4 rotation around z-axis. By incorporating this coupling term, H = H TDSM + H exc is the full Hamiltonian for the electrons in the TDSM.

B. Simulation method
Based on the lattice model defined above, we perform a numerical simulation of the dynamics of the electrons and the magnetization. The aim of this simulation is to evaluate the influence of the magnetization dynamics in FM1 n 1 (t) on the magnetization dynamics in FM2 n 2 (t), which characterizes the spin current transmitted via the TDSM. In order to evaluate them, we simultaneously solve the time-dependent Schrödinger equation for the many-body wave function |Ψ(t) for all the electrons in the TDSM [74,75], and the LLG equatioṅ for the magnetization in FM2 n 2 (t), with γ the gyromagnetic ratio and α the Gilbert damping constant. We introduce the dynamics of n 1 (t) as the input, and do not evaluate the feedback on n 1 (t) from the electron dynamics. The Hamiltonian H(t) for the electrons depends on n 2 (t), and the effective magnetic field B eff (t) for the FM2 depends on |Ψ(t) via the exchange coupling. In particular, if we define the number and magnitude of spins in FM2 as N s and S, the effective magnetic field B eff (t) for each spin is given by where O (t) denotes the expectation value of the operator O evaluated with the many-body wave function |Ψ(t) . Equations (43) and (44) are thus correlated, from which we can evaluate the spin current transmission via the TDSM. As the initial condition for t < 0, we set n 1 (t < 0) = n 2 (t < 0) = e x , and take |Ψ(t < 0) as the Slater determinant of the occupied states in equilibrium, where all the eigenstates in the TDSM below the Fermi energy F = 0 are occupied. At time t = 0, we switch on the in-plane magnetization dynamics in FM1 n 1 (t) = (cos Ωt, sin Ωt, 0), (46) with the precession periodicity T p = 2π/Ω, and solve the time-dependent equations (43) and (44) simultaneously. In order to evaluate the effect of the transmitted spin current exclusively, we neglect the Gilbert damping α and solve Eq. (44) solely with B eff from the exchange coupling. We suppose that the spins in FM2 are residing on the lattice sites at the interface x = L x , which yields N s = N y N z , and fix S = 1 for the simplicity of calculation. Throughout our simulations, we fix N y = 28 and N z = 16.

C. Spin current vs spin torque
Before showing our simulation results, we discuss how the spin torque on FM2 calculated from the simulation is related to the spin current flowing into FM2, to compare the simulation results with the analytical estimations given in the previous section.
The dampinglike torque tilts the magnetization toward zaxis, which originates from the spin angular momentum injected into the magnet. We need to check if the surface-mediated spin current estimated in the previous section is the main contribution to the dampinglike torque τ DL (t). From the discussion in the previous section, the net spin current flowing into FM2 reaches the semi-quantized valuē in our lattice system, if the spin transport is dominated by the helical surface states. On the normalized magnetization n 2 (t) in FM2, this spin current may exert a spin-transfer torque [14] τ where M (tot) = γN s S denotes the net magnetic moment in FM2. Therefore, in order to check if the surface transport picture is valid, we compare the dampinglike torque τ DL (t) with this surface-mediated spin-transfer torquē τ stt (t), by evaluating their ratio ρ(t) ≡ τ DL (t)/τ stt (t). By using the particular settings k D = π/2a, N s = N y N z , and S = 1 employed in our simulation, this ratio can be derived from the time evolution of n 2 (t), which we shall plot in the following figures. If this ratio ρ(t) reaches unity, we can understand that the spin transmission in the TDSM is dominated by its surface states.

D. Results and discussions
We now show our results of the numerical simulations. First, we demonstrate a typical time-evolution behavior of n 2 (t) in Fig. 7(a). We here take the parameters J exc = 0.5u, L x = 16a, and T p = 20u −1 . As mentioned  above, n 2 is fixed to x-direction as the initial condition (t < 0). After the magnetization dynamics n 1 (t) in FM1 is switched on at t = 0, the magnetization n 2 (t) in FM2 also deviates from its initial direction, which implies that spin angular momentum is transmitted from FM1 to FM2 via the TDSM. We can see that the outof-plane component n 2z evolves first at the early stage of the magnetization dynamics, which can be considered as the effect of the dampinglike torque from the transmitted spin current. In order to see the nature of the transmitted spin current in more detail, we plot in Fig. 7(b) the time evolution of ρ(t), namely the ratio of the dampinglike torque τ DL (t) from this simulation in comparison with the spintransfer torqueτ stt (t) from the surface transport picture, with several values of T p . We can see that, for any value of T p in these calculations, ρ(t) reaches unity at the time t ∼ 20u −1 , which implies that the spin current is dominated by the helical surface states of the TDSM, as predicted in Section III C.
The time evolution of the dampinglike torque τ (tot) DL (t) can be associated with the transient steps (b)-(d) described in Section III C (or Fig. 6) as follows. Its zerovalue plateau for t u −1 can be regarded as step (b), with the spin signal from FM1 propagating toward FM2. The semi-quantized plateau for 20u −1 t 40u −1 corresponds to step (c), where the signal from FM1 is reflected by FM2 and is injecting spin angular momentum into FM2. For 40u −1 t, the injected spin current deviates from the semi-quantized value. This behavior can be associated with step (d), where the signal reflected by FM2 returns to FM1 and gradually enters FM2 again, enhancing the spin injection into FM2.
We next investigate how the transmitted spin current is affected by the exchange coupling parameter J exc at the interfaces of the TDSM and the FMIs. Figure 8(a) shows the time evolution of the ratio ρ(t) between τ DL (t) and τ stt (t) for several values of J exc , with T p = 20u −1 . While ρ(t) for J exc 0.4u shows a plateau close to unity at the early stage of the magnetization dynamics, the plateau for J exc = 0.2u is lower than unity. The suppression of the plateau for small J exc can be clearly seen by plotting the time-averaged value of ρ(t), for t ini = 25u −1 and t fin = 40u −1 , which is shown in Fig. 8(b). The dependence on J exc can again be understood from the surface transport picture: the semiquantized spin current is achieved if the magnetization dynamics is adiabatic, which requires the exchange gap 2J exc on the surface spectrum to be much larger than the precession frequency Ω = 2π/T p ≈ 0.3u −1 of the magnetization n 1 (t) (for T p = 20u −1 ).
In order to check the robustness of spin transmission against disorder, we introduce the local random potential where V r takes a random value V r ∈ [−W/2, W/2] for each lattice site r, with W characterizing the strength of the disorder. With 20 profiles of the random disorder potential V r , we simulate the time evoltion of n 2 (t), and take average of n 2 (t) over the 20 profiles to evaluate the disorder-averaged behavior. The time evolution of the ratio ρ(t) = τ DL (t)/τ stt (t) and its time-averaged value ρ ave in the window t = 25u −1 -40u −1 are shown in Fig. 9(a) and (b), with T p = 20u −1 , J = 0.5u, and  L x = 16a. We calculate both the standard errors of the dampinglike torque with respect to the disorder and the disorder-averaged oscillations of the dampinglike torque, and plot the larger one as the error bar for each data point. We can see that the plateau ρ(t) ≈ 1 is achieved for a weak disorder W 2u, due to the robustness of the helical surface states under disorder with time-reversal symmetry. The plateau value is gradually suppressed under a strong disorder, once its magnitude exceeds the bandwidth ∼ 2u of the Dirac bands in the bulk. Finally, in order to check the robustness of the surface spin transport over a long range, we vary the system size L x and observe its effect on the torques. In Fig. 10(a), we show the time-evolution of the ratio ρ(t) = τ DL (t)/τ stt (t) for several values of L x , with the disorder strength W = u. We can see that ρ(t) rises to the plateau ≈ 1 at t rise ≈ L x /v edge (here v edge ≈ ua), as shown in Fig. 10(b), namely the time when an electron propagating from FM1 arrives at FM2. Moreover, the semi-quantized plateau is not significantly violated by the disorder, even for a long distance L x = 24a. From these results, we can conclude that the helical surface states of the TDSM can realize a long-range spin transport that is robust against a moderate disorder below the bulk bandwidth.  junctions of two FMIs and a TDSM as a model setup, as shown in Fig. 1, we have investigated the spin transfer between the two FMIs driven by the magnetization dynamics in one FMI (FM1). We have evaluated the spin transfer both analytically by evaluating the electron numbers in the helical surface channels based on the 1D scattering theory, and numerically by simulating the real-time evolution of all the electrons in the TDSM on a lattice model. As a result, we have found that the spin transfer between the two FMIs at charge neutrality is dominated by the helical surface states, and that such a surface spin transport is almost insensitive to the disorder keeping time-reversal symmetry.
In particular, at the early stage of the spin transmission after turning on the magnetization dynamics, we have found that the transmitted spin current reaches the semi-quantized valuej s 2 , which is a universal value determined by the precession frequency Ω of the magnetization dynamics in FM1 and the number of helical channels ν z on the surface, corresponding to the distance of the Dirac points 2k D in momentum space. This semiquantized spin current is achieved if the magnetization opens a large exchange gap on the helical surface states in comparison with the frequency Ω. This condition is in common with the quantized charge pumping driven by magnetization dynamics on the helical edge states of QSHI [53][54][55][56][57][58][59][60][61]. Indeed, as discussed in Section III, the (semi-)quantized charge current and spin current are described in the unified framework based on the edge (surface) transport picture, and they are related independently of the coupling to the FMIs, as in Eq. (6). Since our analysis and simulation show that the trasmitted spin current will deviate from the semi-quantized value after a long time of magnetization dynamics in FM1, we expect that the semi-quantized spin current can be measured if the magnetization dynamics is in a short time, e.g. driven by a microwave pulse.
From our findings in this article, we expect that the helical surface states of the TDSM are advantageous for a long-range spin transport, in comparison with conduction electrons in normal metals or magnons (spin waves) in magnetic insulators. Plus, in comparison with 1D edge states of 2D topological insulators (QSHIs) and Chern insulators, the helical surface states of 3D TDSMs are advantageous in that they consist of many 1D channels and are capable of transferring a large spin current. The recent transport measurement of a heterostructure of a TDSM Cd 3 As 2 and a FMI indicates the effect of exchange splitting on the surface states [41], and hence we may expect that the long-range spin transport can be possibly measured with such heterostructures of Cd 3 As 2 . We here evaluate the scattering rates R( ) and T ( ) introduced in Section III A, which are defined in the rotating frame of spin. Assuming that the FMI is of length L and located at x = 0, the Hamiltonian for the edge electrons in the rest frame is given as where π L (x ) is the rectangular function taking a value 1 for −L/2 < x < L/2 and 0 otherwise, and p = −i∂/∂x is the momentum operator along x . We take the precession of magnetization around z-axis, where the magnetization direction n(t) is given as n(t) = (sin θ cos(Ωt + φ), sin θ sin(Ωt + φ), cos θ). (A2) We introduce J z = J cos θ and J ⊥ = J sin θ for later discussions. J ⊥ gives the exchange gap, if the magnetization is stationary. The magnetization n(t) is fixed in the rotating frame of spin. By the time-dependent unitary transformation which rotates spin by the angular velocity Ω around zaxis, the Hamiltonian becomes time-independent, as schematically shown in Fig. 11. Therefore, we can apply the conventional scattering theory with a timeindependent scatterer in this rotating frame, to treat the charge and spin pumping by the precessing magnetization. We here fix the energy˜ in this frame, evaluate the eigenstate in each region, and connect the obtained eigenstates to derive the scattering solution.
(i) In the nonmagnetic region, the solutions are simply the eigenstates of spin. The plane-wave solutions for right-moving spin-↑ electrons and left-moving spin-↓ electrons are given as respectively.
(ii) The solutions in the magnetic region are given as where the momentum k ± is defined by and φ ± is the eigenvector of the matrix for the eigenvalue˜ . In particular, one can write φ ± as with without normalization. If˜ is inside the exchange gap (|˜ | < J), k ± becomes imaginary and the solutions become exponentially-growing and decaying functions.
With the above solutions, we can construct the overall wave function as The boundary conditions at x = ± L 2 lead to the relations From the definitions v edge k ↑/↓ = ±˜ + Ω 2 and v edge k ± = ±v edge K − J z + Ω 2 , this relation further reduces as Here we need to recast the above relation into the form of S-matrix, The relation can be rewritten as which yields The relation can be rewritten as The obtained reflection and transmission rates in the rotated frame satisfy the detailed balance relations |r ↑↓ (˜ )| 2 = |r ↓↑ (˜ )| 2 (A35) which are R(˜ ) and T (˜ ) defined in the main text. From the above calculations, these rates are explicitly obtained as where We can check that they satisfy the unitarity condition R(˜ ) + T (˜ ) = 1. (A40)