Bloch oscillations in the spin-1/2 XXZ chain

Under a perfect periodic potential, the electric current density induced by a constant electric field may exhibit nontrivial oscillations, so-called Bloch oscillations, with an amplitude that remains nonzero in the large system size limit. Such oscillations have been well studied for nearly noninteracting particles and observed in experiments. In this work, we revisit Bloch oscillations in strongly interacting systems. By analyzing the spin-1/2 XXZ chain, which can be mapped to a model of spinless electrons, we demonstrate that the current density at special values of the anisotropy parameter $\Delta=-\cos(\pi/p)$ ($p=3,4,5,\cdots$) in the ferromagnetic gapless regime behaves qualitatively the same as in the noninteracting case ($\Delta=0$) even in the weak electric field limit. When $\Delta$ deviates from these values, the amplitude of the oscillation under a weak electric field is suppressed by a factor of the system size. We estimate the strength of the electric field required to observe such a behavior using the Landau--Zener formula.


I. INTRODUCTION
Imagine a system of electrons in a one-dimensional ring [Fig 1(a)]. Suppose that a static and uniform electric field E is turned on at time t = 0. When E is sufficiently weak, the response of the electric current density j(t) is dominated by the contribution from the Drude weight D 1 of the linear response theory [1]. In free space, electrons keep accelerated and the induced electric current density increases linearly as a function of time, i.e., j(t) D 1 Et. In the presence of a disorder-free lattice with the lattice constant a, on the other hand, noninteracting electrons get accelerated initially but then start to move backward [ Fig. 1(b)], forming an oscillatory motion with the frequency aE/(2π) (we set e = = 1 throughout this work). Such an oscillation is known as a Bloch oscillation and has been observed in semiconductor superlattices [2][3][4][5]. The oscillation can be understood as a result of nonlinear responses of the electric current to higher powers of E [6].
The adiabatic transport property is known to be strongly affected by interactions and disorders [7,8]. For interacting systems, recent studies revealed that the Drude weight D N (N ≥ 2) of nonlinear responses tends to diverge in the limit of large system size [6,9,10]. Similar divergence has been observed in noninteracting systems in the presence of a localized impurity [11]. These results imply that the Bloch oscillations in noninteracting and imperfection-free systems may completely disappear when interactions or disorders are added. This is indeed the case for a tight-binding model with a single defect [11].
In this work, we study the current response of disorderfree interacting electrons in the limit of weak electric field. By analyzing the spin-1/2 XXZ chain, which can be mapped to a model of spinless electrons interacting with each other between neighboring sites, we find that interacting systems can exhibit Bloch oscillations * hwatanabe@g.ecc.u-tokyo.ac.jp  (7) and the current density j(t) in (8) for the model (1) with A = A(t) = 2πt/T and L = 24 in the noninteracting limit ∆ = 0. We set T = 1 in this plot but results are almost identical up to T = 10 6 , the largest period we computed. Both j(t) and ε(t) show a single oscillation over the period T where A increases by 2π. Solid curves represent the analytic expressions ε(t) = −L −1 cos A(t)/ sin(π/L) and j(t) = L −1 sin A(t)/ sin(π/L). The inset shows the amplitude jmax as a function of L, suggesting that it is of O(1) in the large L limit.
for fine-tuned values of the interaction ∆ = − cos(π/p) (p = 3, 4, 5, · · · ). When ∆ is not exactly at these special values, the interaction induces level repulsions to the many-body energy levels as a function of the gauge field A. As a result, the behavior of j(t) becomes qualitatively different -it still oscillates as a function of t but with an amplitude inversely proportional to the system size. At the same time, the frequency becomes LE/(4π) and is proportional to the system size. Such behaviors can be seen only when the electric field E is sufficiently weak, and we give an estimate of the required strength of E using the Landau-Zener formula [12][13][14]. In this work, we discuss the spin-1/2 XXZ chain whereŝ α i (α = x, y, z) is the spin-1/2 operator defined on the site i andŝ ± i =ŝ x i ± iŝ y i , L is the system size, J is the antiferromagnetic coupling constant, ∆ is the anisotropy parameter, and A is the vector potential (with the sign convention opposite to the conventional one) associated to the spin rotation symmetry about z axis. The lattice constant a is set to 1. We impose the periodic boundary condition and identifyŝ α L+1 =ŝ α 1 . The conserved current density of this model is defined bŷ The spin chain (1) can be mapped to a model of spinless electronŝ by the Jordan-Wigner transformation [15], where t 0 = J/2 and U = J∆. The spin current density (2) can thus be regarded as the electric current density of spinless electrons. In the following, we set J = 1 and focus on the gapless regime −1 < ∆ < 1.

B. Periodicity of the adiabatic ground state
Let us examine the periodicity of the many-body energy levels ofĤ(A) as a function of A. Obviously we haveĤ(2π) =Ĥ(0), but ∆A = 2π is not the smallest period of this model in the following sense. The gauge field A is related to the flux Φ = AL piercing the hole of the 1D ring [ Fig. 1(a)], and A = 0 and A = 2π/L give the physically equivalent flux. In fact, they are connected by a unitary transformation, called the large gauge transformation [16]. Therefore, the entire spectrum ofĤ(A) has a period ∆A = 2π/L. This is the case regardless of the specific choice of ∆.
In contrast, the "adiabatic ground state" has a longer period that depends on ∆ [17]. Let |Ψ be the ground state at A = 0. Among the stroboscopic eigenstates of H(A), we keep track of the state that is smoothly connected to |Ψ as a function of A. This is what we call the adiabatic ground state |Ψ AGS (A) . See red curves in Fig. 2 for illustration. In the following, we denote by ε AGS (A) the energy density of the adiabatic ground state.
The flux-threading argument for the Lieb-Schultz-Mattis theorem predicts a level crossing between two levels with distinct momentum eigenvalues 0 and π for translation and time-reversal invariant spin-1/2 chains [16]. Hence, the smallest possible values of the period of the adiabatic ground state is ∆A = 4π/L. This is indeed the case for generic ∆ [18,19] [ Fig. 2(a)]. In particular, the period ∆A vanishes in the large L limit and the energy density ε AGS (A) becomes a constant ε AGS (0) independent of A.
Exceptions occur at special values of ∆, i.e., ∆ = − cos(π/p) (p = 2, 3, · · · ), at which the adiabatic ground state has a much longer period [20], as shown in Fig. 2(b) for p = 3. From numerical investigations, we identify the period ∆A = 2π/(p − 1) for the sequence L = 2(p − 1) ( ∈ N) [21]. The system size dependence is examined in more detail by the exact diagonalization in Fig. 3. For example, in the case of p = 3, the period is ∆A = π for L = 4 and ∆A = 2π(L−2)/L for L = 4 +2. The actual period converges to ∆A = π in the large L limit for the latter case as well [see Fig. 3(a)]. As shown in Fig. 3(b) and (c), similar behaviors are also observed for the p = 4 and p = 5 cases. Most importantly, both the period ∆A and the amplitude of the oscillation of ε AGS (A) do not vanish in the large L limit. This difference results in the qualitatively different behaviors of the induced electric current in response to an infinitesimal electric field, as we discuss in the next section.
In Fig 3, the dotted curve represents whereε GS is the energy density andD n is the n-th order Drude weights in the large L limit. The analytic expressions ofε GS andD 1 have long been known [18,22]. The expressions for higher Drude weights,D 3 andD 5 , have been obtained in recent studies [6,9]. For readers' convenience, we include these expressions in Appendix A. We see an excellent agreement for small A, supporting the numerical stability of our results. Assuming that ε AGS (A) converges to a smooth function of A in the large L limit, all the nonlinear Drude weights remain finite for the special values of ∆ [11]. Before closing this section, let us explain how we kept track of the adiabatic ground state |Ψ AGS (A) in our numerical investigation of ε AGS (A). First, we reduce the number of candidate states by using the spin rotation symmetry about z-axis and the translation symmetry, focusing on the S z = 0 sector and the momentum P = 0 sector. Next, we discretize A ∈ [0, 2π] as A (n) ≡ n δA (n = 0, 1, 2, · · · ). The increment δA needs to be chosen sufficiently small, and in this work we set δA = π/(160L). We start with the unique ground state ofĤ(A) for A = A (0) = 0 and A = A (1) = δA. Then, for n > 1, we proceed step by step by choosing as |Ψ AGS (A (n+1) ) the state whose energy density is closest to the linear interpolation among the eigenstates ofĤ(A (n+1) ) in the S z = P = 0 sector. We confirmed the convergence by changing δA.

III. TIME EVOLUTION
A. Bloch oscillations in the weak field limit Now, following Ref. [23], let us introduce a time dependence to the Hamiltonian (1) by setting A = A(t) = 2πt/T for t ≥ 0 and A = 0 for t < 0. This timedependent gauge field induces a static and uniform electric field E = 2π/T for t ≥ 0. We assume that the system is in the ground state |Ψ ofĤ(0) for t < 0. The time evolution of the state for t ≥ 0 is described by the Schrödinger equation The energy density ε(t) and the induced current density j(t) and at time t > 0 are given by When T is sufficiently large for a given L and ∆, the adiabatic theorem [24,25] gives [6] ε(t) = ε AGS (A(t)), Namely, the energy density ε(t) and the current density j(t) of the time-dependent problem are fully characterized by the adiabatic ground-state energy density ε AGS (A) of the time-independent system. This immediately implies the possibility of Bloch oscillations even under an infinitesimal electric field for the special values ∆ = − cos(π/p) (p = 2, 3, · · · ), where the period and the amplitude of the oscillations of ε AGS (A) remain nonzero in the large L limit. We demonstrate this prediction by numerically solving the Schrödinger equation using the fourth-order Runge-Kutta method. We mainly set dt = 0.01 but sometimes use dt = 0.005 if necessary. Our results are summarized in Fig. 4. As anticipated, even for T = 10 6 , the amplitudes of the oscillations of the energy density ε(t) and the current density j(t) are of O(1) for the special values of ∆ = − cos(π/p) in panels (a)-(c). In contrast, they are O(L −1 ) for the generic case in (d). Over the time period T , the gauge field A(t) increases by 2π and the number of peaks is p − 1 for the special values in panels (a)-(c) and is of O(L) for the generic case in (d). These are all consistent with the behavior of ε AGS (A) discussed in the previous section. Sharp drops of the electric current j(t) around t = 2π(2j − 1)/(LE) (j = 1, 2, . . . ) seen in (d) underlie the divergence of nonlinear Drude weights for generic ∆.

B. Landau-Zener formula
When ∆ deviates from the special values, the manybody energy spectrum opens a tiny gap ∆E between the adiabatic ground state and the next level at A = 2π/L. See the orange dashed box in Fig. 2 (a). When T is large enough, the applied electric field is sufficiently weak and the state |Ψ(t) stays in the adiabatic ground state. In contrast, when T gets smaller, non-adiabatic transition starts to occur. This type of tunneling has been studied for the Hubbard ring of spinful electrons in Refs. 26 and 27. Here we estimate the strength of the electric field E required to avoid non-adiabatic transition in the spin-1/2 XXZ chain 1 based on the Landau-Zener formula.
Let us first assume ∆E = 0. Then the energy density of the adiabatic ground state around A = A 0 ≡ 2π/L can be approximated, to the leading order of L −1 , by See Fig. 5(a) for comparison with the exact diagonalization for L = 12 and ∆ = 0.5. The error in the approximation is a finite-size effect. Hence, when A(t) = Et, the energy difference of the two levels that cross each other at A = A 0 goes as α|t−t 0 | with α ≡ 2D 1 A 0 LE = 4πD 1 E and t 0 ≡ A 0 /E. When ∆E > 0, the Landau-Zener formula [12][13][14] gives the probability of the non-adiabatic transition This formula was originally derived for two level systems and ∆E/2 is identified with the amplitude of the offdiagonal components of the two level Hamiltonian. Here, following Ref. 26, we apply it to the focused two levels in the many-body spectrum [28].
To make use of the formula (12), we now study the system-size dependence of the finite-size splitting ∆E. Here we assume the power-law form [29] The leading contribution to ∆E comes from the least irrelevant perturbation due to Umklapp scattering [19,30], which gives a = 4K − 1 = 2π[arccos(−∆)] −1 − 1, where K is the Luttinger parameter [31]. We append the analytic expression for C in Appendix A. For example, we find ∆E 2 √ 3πL −2 for ∆ = 0.5. We confirm these analytic expressions by fitting with numerical data up to L = 32. As shown in Fig 5 (b) and (c), we find a good agreement.
Finally, let us verify the Landau-Zener formula (12). We compute the probability of adiabatic process in two ways: one by 1 − p LZ using the last expression of (12) and the other by the overlap  (14) at t = 2t0 and the solid curve represents 1 − pLZ given by (12).
at t = 2t 0 , by numerically solving the time-dependent Schrödinger equation (6). For the gap ∆E CL −a , we use the analytic expressions of C and a. As shown in Fig. 5(d) for ∆ = 0.5 and L = 12, they agree unexpectedly well. These results suggest that, in order to avoid nonadiabatic transition, the electric field must satisfy The right-hand side becomes small quickly as L increases, suggesting that E must be chosen very small. In other words, when E = 2π/T , the time interval T must be chosen quite long. Therefore, it is actually nearly impossible to achieve the adiabatic limit such as the one in Fig. 4(d) in a thermodynamically large system.

C. Bloch oscillations under a stronger field
The subtle differences caused by ∆ discussed above are irrelevant under a strong electric field. Indeed, we find qualitatively the same behavior regardless of the value of ∆ for the T = 1 case as shown in Fig. 6. The corresponding electric field E = 2π is so strong that non-adiabatic transitions occur everywhere. We also present the behavior of ε(t) and j(t) for T = 1, 10, 10 2 , · · · , 10 6 for comparison in Fig. 7 in Appendix.
In principle, one can analyze the Landau-Zener transition probability for each level repulsion point, but these plots imply that behaviors of the ε(t) and j(t) depend on the details of the system except for the two limiting cases of adiabatic and non-adiabatic transitions.

IV. DISCUSSIONS
In this work, we demonstrated that interacting electrons can in principle exhibit Bloch oscillations even under a weak electric field when the interaction is finetuned. This is the case when the non-linear Drude weights introduced in Ref. 6 remain finite in the large L limit to all orders of response. The expected behavior of the oscillations in the electric current density is similar to the noninteracting limit in the sense that the amplitude of the oscillation remains nonzero in the limit of large system size. Yet, there are also some differences, for example, in the number of peaks and the period of the oscillations for a given adiabatic process, as one can see by comparing Fig. 1 (b) and Fig. 4 (a)-(c). In particular, the frequency of the oscillation becomes a(p − 1)E/(2π) for ∆ = − cos(π/p).
When the interaction is not exactly at the fine-tuned value, the oscillation in the limit of the weak electric field is suppressed in a large but finite system [ Fig. 4  (d)]. However, even in that case, when the applied electric field is not weak enough, non-adiabatic transitions occur and one observes nontrivial oscillations of the electric current density, as shown in Fig. 6. In fact, our estimate based on the Landau-Zener formula in (15) suggests that the adiabatic limit is difficult to be realized in a thermodynamically large system. This, in turn, implies that pathological behaviors of nonlinear Drude weights observed in Refs. [6,9,10] do not come into play for a realistic strength of the applied electric field.

ACKNOWLEDGMENTS
We thank Kazuaki Takasan and Masaki Oshikawa for collaborations in previous related works. We thank Hosho Katsura for partially informing us of their on-going result based on the Bethe ansatz, which supports the finiteness of nonlinear Drude weights we numerically observed [32]. We also thank Masaki Oshikawa and Hosho Katsura for their valuable comments on the earlier version of the draft. Here we summarize the analytic expressions derived in previous studies. We parametrize ∆ as ∆ = cos γ with 0 ≤ γ < π. The ground state energy densityε GS and the linear Drude weightD 1 in the large L limit are given by [18,22] .