Quantum Spin Pump in S=1/2 antiferromagnetic chains -Holonomy of phase operators in sine-Gordon theory-

In this paper, we propose the quantum spin pumping in quantum spin systems where an applied electric field ($E$) and magnetic field ($H$) cause a finite spin gap to its critical ground state. When these systems are subject to alternating electromangetic fields; $(E,H)=(\sin\frac{2\pi t}{T},\cos\frac{2\pi t}{T})$ and travel along the {\it{loop}} $\Gamma_{\rm{loop}}$ which encloses their critical ground state in this $E$-$H$ phase diagram, the locking potential in the sine-Gordon model slides and changes its minimum. As a result, the phase operator acquires $2\pi$ holonomy during one cycle along $\Gamma_{\rm{loop}}$, which means that the quantized spin current has been transported through the bulk systems during this adiabatic process. The relevance to real systems such as Cu-benzoate and ${\rm{Yb}}_4{\rm{As}}_3$ is also discussed.


Introduction
Transport phenomena in magnets such as the colossal magnetoresistance 1 and anomalous Hall effect [2][3][4][5][6][7][8] are long-standing subjects in solid state physics. Recently, exotic characters of the anomalous Hall effect are closed up in ferromagnets, [9][10][11][12] antiferromagnets 13 and spin glass 14,15 with non-coplanar spin configurations and collinear ferromagnets with strong spin-orbit couplings. [16][17][18][19][20] Their common origin is gradually recognized as the topological character of magnetic Bloch wavefunctions in its ordered phase, sharing the same physics with the two dimensional (2D) quantum Hall effect. Thereby the spontaneous Hall conductivity is related to the fictitious magnetic field defined in the crystal momentum space and its quantized value is directly related to the first Chern number associated with the fiber bundle whose base space is spanned by crystal momenta.
Meanwhile, in early 80's, Thouless and Niu proposed the quantized adiabatic particle transport (QAPT) in gapped fermi systems, where an adiabatic sliding motion of the periodic electrostatic potential in one dimensional (1D) system pumps up an integer number of electrons per cycle. 21,22 In fact, this quantized particle transport is another physical manifestation of the topological character associated with the Bloch wavefunctions. Specifically this particle (polarization) current is related to the fictitious magnetic field defined in the generalized crystal momentum space which is in turn spanned by the crystal momentum and deformed parameters. Because of this analogy to the 2D quantum Hall currents, this 1D quantized particle current is known to be topologically protected against other perturbations such as disorders and electron-electron correlations. 22 Due to this peculiar topological protection, the QAPT became recently highlighted in the realm of the spintronics, where electron pumpings have been experimentally realized in meso and/or nano-scale systems. [23][24][25][26] However the quantized * E-mail address: shindou@appi.t.u-tokyo.ac.jp.
value observed in these experiments should be attributed to the Coulomb blockade 23,25,26 and is different form the aforementioned Thouless quantization, which is assured only in the thermodynamic limit. 23 In this paper, based on the original idea of the QAPT, we propose the quantized spin transport (quantum spin pump) in the macro-scale magnets. There we discuss systematically the topological stability of its quantization. Specifically, we clarify in this paper the relation between its quantization and holonomy of the phase operator of sine-Gordon theory in 1D quantum field theory. By using this relation, its topological stability can be quantitatively judged by seeing whether the expectation value of this phase operator acquires 2π holonomy during the cyclic process or not.
As for the experimental relevant systems to our theory of quantized spin transports, we have the S = 1 2 quantum spin chain such as Cu-benzoate and Yb 4 As 3 (charge ordered phase), where its unit cell contains two crystallographically inequivalent sites. Accurately speaking, both the translational symmetry by one site (T : S j → S j+1 ) and the bond-centered inversion symmetry which exchanges the nearest neighboring sites (I bond : S i−j ↔ S i+j+1 are crystallographically broken. As a result of this peculiar crystal symmetry, the g-tensors of its two sublattices are in general different in these systems; (1) where [g u ] and [g a ] represent the uniform and staggered component of the g-tensors respectively. 27 Late 90's, these quantum spin systems, whose ground states are nearly critical because of their strong one dimensionalities, are experimentally revealed to show the spin gap behaviors under the external magnetic field. Furthermore, following theoretical works showed that the field-induced gap behaviors (∆ ∼ H 2/3 ) are originated from this staggered component of the g-tensor, where the effective staggered magnetic field in the eq. (1) endows its critical ground state with a finite spin gap. 29,30 On the other hand, the exchange interaction J in these S = 1/2 quantum spin chains, J i S i · S i+1 , does not have an alternating component, since these systems are invariant under the site-centered inversion symmetry which exchanges the nearest neighbor bonds (I site : S i−j ↔ S i+j ). However, when we break this symmetry by applying an electric field E along an appropriate direction, 31 the exchange interaction in general does acquire a staggered component; Futhermore a site-centered inversion operation with the sign of E reversed requires that ∆ must be an odd function of E. This dimerizing field ∆ is also expected to induce a finite mass to its critical ground state, causing the spin-Peierls state. Based on these observations, we propose a method of generating the quantized spin current in this type of spin systems (T, I bond : broken, I site : unbroken) by using electromagnetic fields. Specifically, we study a following S = 1/2 Heisenberg model with time-dependent staggered field and bond alternation: where staggered Zeeman field h st (t) and bond alternations ∆(t) are supposed to be controlled by the applied electromagnetic fields. In § 2, we will argue by using bosonization technique that the quantized number of zcomponent of spins are transported from one end of the system to the other during one cycle along the loop which encloses the critical ground state at (∆, h st ) = (0, 0) in the ∆-h st plane. In order to uphold this bosonization argument, we also demonstrate the numerical calculations in § 3, where we evolve the ground state wavefunction according as the above time-dependent Hamiltonian. Thereby we attibute the quantized spin transport to the Landau-Zener tunnelings 40 which indeed happen between several energy levels during this cycle process.
In addition to eq.(3), we also have a uniform component of Zeeman fields in real systems, H(t) · S j , whose magnitude is usually larger than that of the staggered component h st (t). This uniform field, especially its x or y component, seems to easily flip the spins accumulated at both boundaries and might spoil the physical consequence of the spin transport. Then we also discuss in § 3 the effect of this uniform field and the remedy against it, where it turns out that the sweeping velocity should be appropriately chosen so as to avoid the relaxation process via this uniform Zeeman field.

Bosonization
For clarity, let us first illustrate the physics of our quantum spin pumping in a simple limiting case; where we take the direction of the staggered Zeeman field as the z-direction. 32 In this simplification, we neglected the exchange coupling between the z-component of spins.
In this limit, we can get a quadratic form of Hamiltonian in terms of the Jordan-Wigner (JW) fermionŜ z Then H XY forms a cosine band in its momentum space, ǫ(k) = −J cos kα (k,α are crystal momentum and lattice constant). In the absence of ∆ and h st , the fermi points locate at k = ± π 2α since fermions in the ground state fill up all the k points but those with positive energy ǫ(k).
Nonzero ∆ and/or h st introduce a finite gap at these two Fermi points. The dimer state induced by a finite ∆ can be understood as a Peierls insulator, where the JW fermions occupy the bonding orbitals between the two neighboring sites and form a valence band in its k-space. On the contrary, the antiferromagnetic state induced by the effective staggered magnetic field h st along the zdirection can be interpreted as an ionic insulator, where the JW fermions stay on every other site. The conduction band and the valence band touch at k = ± π 2α , only when (∆, h st ) is taken at the origin in this ∆-h st plane. Because of the periodicity of π α along the k-axis, the double degeneracy point at (k, ∆, h st ) = ( π 2α , 0, 0) is identical to that of (k, ∆, h st ) = (− π 2α , 0, 0). Elementary analyses show that this double degeneracy point becomes the source (sink) of the vector field B +1 (B −1 ) defined in the generalized crystal momentum space (k-∆-h st space); where K = (k, ∆, h st ) and ∇ K = (∂ k , ∂ ∆ , ∂ hst ). Here |n(K) represents the periodic part of the Bloch function for the valence band (n = +1) and the conduction band (n = −1). The vector field A n,µ changes by 1 2π ∇ Kµ φ under the U(1) gauge transformation of these Bloch wavefunctions; |n(K) → exp[iφ(K)]|n(K) , while its rotation, i.e., B n remains invariant. Because of this gauge invariance, we often call the latter vector field B n as a fictitious magnetic field or flux. Correspondingly, the double degeneracy point located at (± π 2α , 0, 0) will be referred to as the fictitious magnetic charge (magnetic monopole in the generalized momentum space), whose magnetic unit can be shown to be 1: Here S 1 represents the arbitrary closed surface which encloses the double degeneracy point at (k, ∆, h st ) = (+ π 2α , 0, 0) ≡ (− π 2α , 0, 0). Based on these observations, let us consider the adiabatic process where the two parameters (∆, h st ) are changed along a loop Γ loop enclosing the origin ((∆, h st ) = (0, 0)); Then, according to the original idea of the QAPT, 21,22,35 the total number of JW fermions (I) which are transported from one side of this system to the other (in the positive direction) during this adiabatic process is equal to the total flux for the valence band (n = 1), B +1 , which penetrates the 2D closed sphere spanned by Γ loop and the Brillouin zone : The physical meaning of this quantized fermion transport is nothing but the quantized spin transport, since the JW fermion density is related to the z-component of the spin density. In other words, the total S z around one end of this system decreases by 1 while that of the other end increases by 1 during this adiabatic cycle. This quantization is topologically protected against the other perturbations as long as the gap along the loop remains finite, 22,35 in other words, as far as the double degeneracy points do not get out of (enter into) the 2D closed surface Γ loop × [− π 2α , π 2α ]. Then, we naturally expect that this quantized spin transport is stable against the weak exchange interactions between the z-components of spins; In the following, we will prove that this is indeed the case, by introducing the exchange interactions between S z . As a first step, we will treat H XY term as a nonperturbed term and review which perturbations are relevant among H Z , H dim and H st by using the bosonization technique. Namely, we first introduce the slowly varying fields, R(x), L(x): Then we rewrite the spin operator in terms of these fermion fields R(x j ) and L(x j ): where : ... : stands for the normal order. 36 Here we retain the lowest order terms with respect to the lattice constant α both for the uniform and staggered components. According to the bosonization recipe, 37 we can rewrite these fermion operators by using the phase operatorθ + (x) and its canonical conjugate fieldΠ(x). Namely, the spin operators in eqs. (12) and (13) read; Then we substitute these equations into eqs.(4)− (6) and (11) and obtain the following expressions for H XY , H Z , H dim and H st in the continuum limit, where we neglect the rapid varying components such as j (−1) j cosθ + (x j ) and etc.: Here the origin of the third term of eq.(17) were discussed elsewhere, 38 which we will omit by redefining ∆ in eq. (18) in the followings. Then our final expression for Hamiltonian reads, where the velocity v is given by v = J 2 (1 + 2|γ| πα ). We also introduced the quantum parameter η as which measures the strength of the quantum fluctuation. Namely, when |γ| in eq.(11) grows from 0, this parameter decreases monotonically from 2. Even though eqs. (16)−(21) were derived in the weak coupling limit (|γ| ≪ 1), the final form of eq. (20) is known to be valid until |γ| reaches 1 (isotropic Heisenberg model), where the SU(2) rotational symmetry of correlation functions ( S x j S x i ∼ S z j S z i ) requires that η = 1 at |γ| = 1. 36 By using this quantum parameter η, the renormalization group (RG) eigenvalue of sin(θ + + ϕ) is represented by 2 − η 2 while that of cos 2θ + is given by 2 − 2η. This indicates that, as long as |γ| ≤ 1 (1 ≤ η ≤ 2) , the exchange coupling between the z-component of spins ( |γ|J 2π 2 α 3 cos 2θ + ) is always irrelevant in the sense of renormalization group analyses, while the dimerizing field ∆ and staggered field h st ( R πα 2 sin(θ + + ϕ)) are equally relevant and lock the phase operatorθ + on π 2 − ϕ + 2πn. As the system is always locked by sin(θ + + ϕ) for |γ| ≤ 1, we naturally expect that the quantized spin transport in the case of |γ| = 0 could be generalized into the case of finite |γ|, at least up to |γ| = 1. This expectation is easily verified when we notice the physical meaning of the phase operatorθ + 2π . Since the spatial derivative of the phase operator corresponds to the zcomponent of spin density, this phase operator is nothing but minus of the spatial polarization of the z-component of spins, i.e., −P S z ≡ − 1 N N j=1 jŜ z j . The equivalence between these two quantities is discussed in detail in the appendix. Then, through the adiabatic process eqs. (8) and (9), θ + decreases monotonically and acquires −2π holonomy after one cycle (see Fig. 1). In other words, P S z increases by 1 per one cycle, i.e., This relation always holds as far as the system is locked by the sliding potential sin(θ + + ϕ(t)), which is true for |γ| ≤ 1. Then we will argue the physical consequence of eq. (22). Generally speaking, when the bulk system has a finite spin gap, the effects of boundaries range over its magnetic correlation length from the both ends. Then we naturally divide the contribution to P S z into following two parts; where the j-summation in P in S z are restricted within the interior of the systems (Ω in ), while spins within the edge region (Ω edge ) contribute to P edge S z . The "interior" of the system is defined as the region where the spin state is same as that of the periodic boundary condition (p.b.c.). Namely, the presence of the boundaries does not influence on the spin state of the interior. On the other hand, the spin state within the "edge" region is strongly affected by the boundaries. Then, by definition, the spins within the interior would turn back to the same spin configuration as that of the initial state after an adiabatic evolution along Γ loop , δP in S z = 0. Therefore the difference of the spatial polarization of Ŝ z j should be attributed to the change of spin configurations in the edge region, In eq.(23), we next approximate j for sites around the right end by N/2 and those of left by −N/2; where we take the origin of site index j at the center of the system. This approximation is allowed, since the edge region ranges over the magnetic correlation length from both ends, which is at most several sites in gapped spin systems. Thereby the semi-equal sign in eq.(24) could be safely replaced by the equal sign in the thermodynamic limit. When we bear in mind that the system has been staying on the eigenspace ofŜ z tot = 0 during this process; eq. (24) indicates that the total S z around the right end (j ≃ N/2) increases by 1 while that of the left end (j ≃ −N/2) decreases by 1 along this cyclic evolution: Namely, the spin current quantized to 1 flows from the left end of the system to the right end during this one cycle.
In summary, the quantized spin transport is always assured as long as |γ| ≤ 1. This story does not alter even with the finite sliding speed c = 2π T , as far as it is slower than the height of this locking potential, in other words, the spin gap along the loop.

Numerical Analyses
In order to confirm the above argument, we performed the numerical calculation in the case of |γ| = 1. Namely, we generated numerically the temporal evolution of the ground state wavefunction under the following timedependent Hamiltonian: where T represents the time order operator, T is the cycle period during which the system has swept the loop (Γ loop ) one time and |g I is the ground state wavefunction at (h st , ∆) = (R, 0). In Fig. 2, we show the expectation value ofŜ z j taken over the final state, |φ(t = T ) , both with the p.b.c. and with the open boundary condition. The parameter R is fixed to be 0.3 and J is taken as 1.5 (0.6) in Fig. 2a(b). The cycle period T is taken as 40 for the former and 80 for the latter, both of which are sufficiently long compared with the inverse of the spin gap observed along the loop.
For the p.b.c., as far as the sweeping velocity is low enough, the final state gives completely same configurations of Ŝ z j as that of the initial antiferromagnetic ground state. However in the case of the o.b.c., the spins around both boundaries are clearly modified as seen in Fig. 2. When we read the total S z around the left end ( 3 j=1 Ŝ z j ) and that of the right end ( 14 j=12 Ŝ z j ) from Fig. 2a(b), the former increases by 0.67 (0.9) while the latter decreases by 0.67 (0.9) after this cyclic process. This result is qualitatively consistent with our preceding arguments.
In order to understand the difference between the result of the o.b.c. and that of the p.b.c., we also calculated the instantaneous eigenenergy as a function of ϕ ≡ 2πt T in Fig. 3, where we take the offset as the ground state energy. For the p.b.c. (inset of Fig. 3a) , there is always a finite energy gap from the ground state for all ϕ. On the contrary, the gap for the o.b.c. reduces strongly at ϕ = 3π 2 (Figs. 3a and 3b). Furthermore, this reduced energy gap ∆ 1 becomes smaller and smaller when we take the system size N larger (inset of Fig. 3b, N=6,8,10,12),  (25)) and dimers are formed between them, i.e., where S i = S i+1 represents the singlet bond. Then, in the case of the o.b.c., the spin at j = 1 and that of j = N cannot form a spin singlet due to the absence of S 1 · S N term. Namely, in the limit of J = R, the following four states would be degenerated, The last two states ((28),(29)) do not belong to the eigenspace of S z tot = 0 and we ignore them henceforth. 39 When we introduce finite J − R > 0, the N 2 -th order perturbation in terms ofŜ + l is expected to be a magnetic correlation length. In fact, we can fit the size dependence of ∆ 1 by this exponential function as in the inset of Fig. 3b, from which the correlation length at ϕ = 3π 2 is estimated around 3 sites. Because of this size dependence, we expect that ∆ 1 would Online-Journal Subcommittee vanish when we took the system size sufficiently large compared with this correlation length.
Because of this crossing character of the energy spectrum at ϕ = 3π 2 , the system which has resided on the ground state (n = 0) transits, at ϕ = 3π 2 , into the first excited state with n = 1 which is represented by the bold solid line in Fig. 3b. We can see this transition in Fig. 4a, where the projected weights of |φ(t) onto each instantaneous eigenstate are given as a function of ϕ = 2πt T . Futhermore, its transition probability, P n=0→n=1 , is well fitted by the Landau-Zener formula 40, 41 as in Fig. 4b, where we plotted P n=0→n=1 for various cycle period T and various system size. From this fitting, c in the Landau-Zener function is estimated around 2.6± 0.2. 42 As shown in the Fig. 4a, after the first transition took place at ϕ = 3π 2 , the second transition from the first excited state to the second excited state with n = 2 (represented by the bold broken line) happens at ϕ ≃ 5.9. However its transition probability P n=1→n=2 is around 70% and not so perfect compared with P n=0→n=1 . This is mainly because the gap ∆ 2 around ϕ ≃ 5.9 between the first excited state and the second excited state is a substantial amount compared with ∆ 1 (see Fig. 3b).  However, as in the inset of Fig. 3b, its size dependence is also scaled by where the magnetic correlation length l is estimated around 4.6 sites. Therefore, when we take the system size large enough, we can naturally expect that the gap ∆ 2 would reduce exponentially, which enhances the transition probability P n=1→n=2 drastically. In order to verify this expectation, we perform another numerical calculations. As it is very difficult to calculate numerically with larger system size (N > 14), we, instead, decrease the ratio J/R in order to reduce the magnetic correlation length l, which is also expected to reduce these gaps according to eq.(30). The result for J = 0.6, R = 0.3 and N = 14 is summarized in Fig.  5. Figure5a indicates the energy spectrum for the lower four excited states withŜ z tot = 0 and Fig. 5b shows the projected weight, | n|φ(t) | 2 for n = 0, 1, 2, 4. Figure5c represents the expectation value ofŜ z j taken over the |n = 4 ϕ=2π≡0 as a function of site index j. Then, as is expected, the gaps we observed in Fig. 3 (∆ 1 , ∆ 2 ) decrease drastically. In addition to them, there appears a new gap ∆ 3 between the excited state with n = 2 and that with n = 4 at ϕ = 6.15. This gap ∆ 3 is also very tiny and scales as exp(−N/2.25).
As a result of these crossing characters of energy spectrum at ϕ = 3π 2 , 5.97 and 6.15, three Landau-Zener type tunnelings take place as shown in Fig. 5b. Namely, the wavefunction |φ(t) changes its weight from the ground state with n = 0 to the first excited state with n = 1 (represented by the bold solid line in Fig. 5a) at ϕ = 3π 2 , from the first excited state to the excited state with n = 2 (bold broken line) at ϕ = 5.97 and lastly from the excited state with n = 2 to the fourth excited state with n = 4 (bold dotted line) at ϕ = 6.15. All these transition probabilities are almost 100% in accordance with our previous expectations. 43 After all, through these transitions, the wavefunction is raised from the ground state onto the fourth excited state (n = 4) during this last quarter of the cycle. When we see the expectation value ofŜ z j taken over this fourth excited state at ϕ = 0 ≡ 2π (see Fig. 5c), it is almost identical to the solid line of the lower panel in Fig. 2, where the totalŜ z on the left end ( 3 j=1 Ŝ z j ) decreases by 0.9 while that of the right end ( 14 j=12 Ŝ z j ) increases by 0.9 in comparison with that of the initial ground state wavefunction.
To summarize, the numerical observations found in Fig. 2 can be ascribed to the difference between the energy spectrum with the p.b.c. and that with the o.b.c. In the latter case, a certain excited state decreases its energy level until ϕ = 3π 2 and then picks up the system which has been staying on the ground state. Then, through the Landau-Zener tunnelings, several excited states carry the system hand in hand perfectly upon a particular excited state at ϕ = 0, whoseP S z increases by almost 1 in comparison with that of the initial state. These observations remain unchanged if the sweeping velocity is always higher than This lower limit for the sweeping velocity is expected to vanish in the thermodynamic limit, since the ∆ 1,2,3 reduces exponentially to zero as the system size N becomes larger. Therefore the conclusion in eq. (22) is consistent with the numerical results of this section.
Finally let us mention about the relevance of our story to real systems such as Cu-benzoate and Yb 4 As 3 . There the effective staggered field along the z-direction accompanies the uniform effective Zeeman fields along the x(or y)-direction, H x(y) u . Furthermore, since these systems break the bond-centered inversion symmetry, there is a Dzyaloshinskii-Moriya (DM) interaction. Its DM vector must alternate bond by bond in the absence of the electric field. In addition to them, the uniform component of the DM vector would be also induced by the applied electric field in general. Even if these terms are gradually introduced to our model calculations, δθ + remains invariant as far as the gap along the loop remains open, in other words, the system along the loop remains locked by the sliding potential sin(θ + + ϕ(t)). This is a manifestation of the fact that the quantized δθ + (in the thermodynamic limit) is a topologically protected quantity. However, detail studies on the effect of these terms might be remaining interesting problems.
When the system travels around this critical ground state M times, total z-component of spin around one end increases by M , while that on the other end decreases by a same amount. Such an inhomogeneous magnetic structure can be detectable around sample boundaries and/or magnetic domain boundaries by using some optical probes. However, we must mention that the excited state with P S z = m (In the case of m = 1, this corresponds to |n = 4 ϕ=0 in Fig. 5) might fall into the ground state |n = 0 ϕ=0 by way of the x-component of the uniform magnetic field H x u at around t = mT , when we sweep the system slower than 2µ B H x u . In order to avoid this relaxation process, the sweeping velocity should be taken faster than 2µ However, in order to keep the system on a particular minimum of the locking potential cos(θ + + 2πt T ), the velocity must be also slower than the spin gap observed along the loop, Therefore we need to have a finite window between the lower limit of the sweeping velocity and the upper limit. In the case of Cu-benzoate, the spin gap of 0.1meV is induced by the applied magnetic field of 1T (2µ B × 1T ∼ 0.1meV). Thereby, the lower limit of this sweeping velocity is unhappily comparable to the higher limit in this system. However, spin gaps around 0.1meV could be also induced by the magnetic field smaller than 1T in those quantum spin chains with relatively larger J, where substantial amount of the window between the upper limit and the lower limit would exist. Therefore, it is not so hard to realize our story experimentally on the quasi-1D quantum spin systems with relatively large intrachain interaction and with two crystallographically inequivalent sites in its unit cell.
We also want to mention about the magnitude of the electric field required in order to induce the dimerization gap enough to be observable and also enough to be comparable to the spin gap induced by the magnetic field, which is around 0.1meV. Since it goes beyond the scope of this paper to quantify microscopically the change of the exchange interaction induced by external electric fields, we will pick up some reference data from other materials. Let's see magneto-electric materials, where the applied electric field often changes its effective exchange interactions J due to its peculiar crystal structure and its change ∆ could be quantitatively estimated via the resulting magnetization. 45 In the case of famous magnetoelectric material Cr 2 O 3 , the intra-sublattice exchange interaction for the + sublattice and that for the − sublattice should be equal to each other because of the inversion symmetry (J + = J − = J). Then an electric field applied along its principle axis breaks this symmetry and induce a finite difference ∆ = J + − J − , where ∆/J was estimated around 10 −5 under E ∼ 1kV/cm. 45 Namely, with its linear coefficient 10 −5 [cm/kV]. When we apply this coefficient to our Heisenberg model, eq.(2) with J ∼ 0.1eV, we need an electric field of order of 1kV/cm in order to induce the dimerizing gap J 1/3 ∆ 2/3 of order of 0.1meV, which is still one order of magnitude smaller than the typical value of the Zener's breakdown field (≥ 10 4 V/cm at T ≤ 100K) of 1D Mott insulators. 46

Conclusion and discussion
In this paper, we propose the quantum spin pumping in S = 1/2 quantum spin chains with two crystallographically inequivalent sublattices in its unit cell (T, I bond : broken and I site : unbroken) and with relatively large intrachain exchange interactions. Due to its peculiar crystal structure, an applied electric field (E) and magnetic field (H) endow its critical ground state with the finite spin gap via the dimerizing field (∆ ∼ E) and the staggered magnetic field (h st ∼ H), respectively, where the phase operatorθ + in sine-Gordon theory is locked on a particular valley of the relevant potential h 2 st + ∆ 2 sin(θ + + ϕ) : θ + = π 2 − ϕ + 2πn. When the system is deformed slowly along the loop which encloses the critical ground state in the E-H plane, the locking potential h 2 st + ∆ 2 sin(θ + + ϕ) in the sine-Gordon model slides gradually (see Fig. 1). After one cycle along this loop, the expectation value of the phase operatorθ + for the ground state acquires 2π holonomy. This means that a spatial polarization of z-component of spins P S z increases by 1 after this cycle. 56 In other words, the zcomponent of the spin current quantized to be 1 flows through the bulk system during this cycle.
These arguments are supported by the numerical analyses, where we performed the exact diagonalization of the finite size system and developed the ground state wavefunction along this loop in a thermally insulating way.

Acknowledgments
The author acknowledges A. V. Balatskii, Masaaki Nakamura, M. Tsuchiizu, K. Uchinokura, T. Hikihara and A. Furusaki, S. Miyashita, S. Murakami and N. Nagaosa for their fruitful discussions and for their critical readings of this manuscript. The author is very grateful to K. Uchinokura for his critical discussions which are essential for this work.
Appendix: Equivalence between the phase operator and the spatial polarization operator In this appendix, we will show that the continuum limit ofP S z corresponds to the phase operator − 1 N α θ +(y) 2π dy, where related work was done by Nakamura and Voit. 47 Historically the polarization operator is known to be an ill-defined operator on the Hilbert space with the periodic boundary condition (p.b.c.) and is consequently formidable to treat in solid state physics. However, recently King-Smith and Vanderbilt 51, 52 and Resta [48][49][50]53 showed that the derivative of the polarization with respect to some external parameter λ can be expressed by the current operator which is in turn welldefined in the Hilbert space with the p.b.c.. Based on this observation, they quantitatively estimated electronic contributions to macroscopic polarizations in dielectrics and semiconductors such as GaAs, KNbO 3 and III-V nitride. Following their strategies, we will show the equivalence between the derivative of − 1 N α dy θ +(y) 2π with respect to λ and that of P S z .
According to the standard perturbation theory, 49 the derivative of the polarization with respect to some mechanical parameter λ is given by where we assume the system has a finite spin gap from the ground state. The time derivative of P S z in the numerator reduces to the summation of the current opera-tor over all bonds; Thereby the contribution of every single bond is always O(1/N ), which is not the case withP S z = 1 N N j=1 jŜ z j . As a result, we are allowed to estimate with the periodic boundary condition the r.h.s. of eq. (A·2) in the thermodynamic limit, which we should have figured out with the open boundary. This is because their difference is attributed to the spin currents on the several bonds 54 around the boundaries, which is at most O(1/N ) in our spin gapped system. Then we will consider eq. (A·2) with the p.b.c. and rewrite it in terms of the bosonization language. That is to say, we express the current operator in eq. (A·2) in terms of the phase operatorθ + (x) and its conjugate fieldΠ(x), Accordingly, the time derivative of P S z in the continuum limit can be expressed in the following form, From now on, we will show that this expression is identical to the time derivative of the phase operator: The first term in the r.h.s. of eq.(A·4) comes from the commutator betweenθ + and H XY given in eq.(16), On the other hand, the commutator between the phase operator and H dim given in eq.(18) vanishes and does not produce the last two terms of eq. (A·4). This is because eq. (18) is not an accurate expression for H dim in the continuum limit and needs some additional terms, whose commutators with the phase operator correctly yield the last two terms of eq. (A·4). In order to prove that this is indeed the case, we will carefully reexamine the bosonization ofŜ + jŜ − j+1 +Ŝ + j+1Ŝ − j 55 : = (r.h.s. of eq. (15)) Then, by taking its staggered components, we obtain the following expression for H dim , instead of eq. (18), Here we want to mention that, irrespective of this modification, cosθ + in H dim always locks the phase operator in combination with sinθ + in h st as far as |γ| ≤ 1 (see the text). This is because the RG eigenvalues of the additional terms such as (∂ xθ+ ) 2 cosθ + , (Π) 2 cosθ + and etc. are all negative and thereby irrelevant in the sense of the renormalization group study. However these additional terms in H dim cannot be discarded, since the commutator between these terms and the phase operator produces the last two term of eq.(A·4): After all, by comparing eqs.(A·6) and (A·9) with eq.(A·4), we can safely replace −i[P S z , H] by 2π , H] in the continuum limit and rewrite eq. (A·1) into a following simple form by using the bosonization language: From this paper, we change its notations as follows; [φ R (x), φ L (y)] = iπ 38) S.Eggert and I.Affleck, Phys. Rev. B 46 (1992) 10866 39) The energy level which reaches around 0.9 at ϕ = 2π in Fig.3a corresponds to doubly degenerate eigenstates which belong to the eigenspaces ofŜ z tot = ±1 respectively. These doubly degenerate states at ϕ = 3π 2 corresponds to the two state given in eqs. (28) and (29). On the contrary, the eigenstate labeled as n = 1 (represented by the bold solid line in Fig. 3b) and the ground state belong to the eigenspace ofŜ z tot = 0. Since our time-dependent Hamiltonian conserves the total spin density S z tot , |φ(t) does not have projected weights onto those eigenstates with S z tot = 0. 40) C. Zener: Proc. R. Soc. London, Ser. A137 (1932) 696 41) S. Miyashita: J. Phys. Soc. Jpn. 64 (1994) 3207 42) We identify P n=0→n=1 as 1 − | φ(t = T )|n = 0 ϕ=2π | 2 . 43) Here the wavefunction does not have the projected weight on the third excited state with n = 3 (thin line in Fig. 5a) during the cycle. 44) The uniform component of the Zeeman field does not cause the spin gap to the critical ground state of the isotropic Heisenberg model as far as it is small. In contrast the staggered component is always relevant even if it is infinitesimally small. However when the uniform component of the magnetic field Hu is very strong, the system will be governed by the fixed point where the ground state is ferromagnetically polarized with a finite spin gap. This fixed point is different from the fixed point at which the system is locked by the staggered component of the magnetic field, i.e., hst sinθ + . As far as the system under the applied magnetic field is governed by the latter fixed point, δθ + is quantized to be 1.
within its magnetic correlation length from both ends. 55) In the final line of eq.(A·7), we drop those terms which can be summarized into a form of the total derivative. 56) P. Sharma and C. Chamon: Phys. Rev. Lett. 87 (2001) 096401; ibid. Phys. Rev. B 68 (2003) 035321. Similar arguments as we have employed in § 2 are also found in this work, where onedimensional Luttinger liquid (LL) between two separated contacts reduces to two decoupled LL wires due to the back scattering potential in between. Namely, the spin part of their back scattering potential is also described by sin(θ + (0) + t), where θ + (0) denotes the phase operator at the midpoint of two contacts. This backscattering potential turns out to be relevant in the presence of repulsive electron-electron interactions and tunneling probability of spins between two contacts becomes vanishing. As a consequence of this vanishing tunneling probability of spins (TPS), number of spins transported between these two contacts is also quantized in their work. Even though our work treats a pumping in macro-scale systems and is different from their work which is about meso (nano)-scale systems, there is an interesting correspondence between these two works. Namely, in both stories, quantizations of spin transports result from the vanishing TPS between two "contacts". Two contacts in their configurations correspond to the two sample boundaries in this paper and TPS between these two "contacts", namely tunneling probability between the state of eq. (26) and that of eq.(27) decreases exponetially in thermodynamic limit, which is also a necessity of the quantized spin transfers in this paper (see § 3).