Fine structure of the nonlinear Drude weights in the spin-1/2 XXZ chain

We study nonlinear Drude weights (NLDWs) for the spin-1/2 XXZ chain in the critical regime at zero temperature. The NLDWs are generalizations of the linear Drude weight. Via the nonlinear extension of the Kohn formula, they can be read off from higher-order finite-size corrections to the ground-state energy in the presence of a $U(1)$ magnetic flux. The analysis of the ground-state energy based on the Bethe ansatz reveals that the NLDWs exhibit convergence, power-law, and logarithmic divergence, depending on the anisotropy parameter $\Delta$. We determine the convergent and power-law divergent regions, which depend on the order of the response $n$. Then, we examine the behavior of the NLDWs at the boundary between the two regions and find that they converge for $n=0, 1, 2$ $({\rm mod}~4)$, while they show logarithmic divergence for $n=3$ $({\rm mod}~4)$. Furthermore, we identify particular anisotropies $\Delta=\cos(\pi r/(r+1))$ ($r=1,2, 3,\ldots$) at which the NLDW at any order $n$ converges to a finite value.


I. INTRODUCTION
Transport phenomena have been a subject of central interest in condensed matter physics. In particular, anomalous transport properties of one-dimensional quantum many-body systems have been actively investigated since they are quite different from higher dimensional ones [1][2][3][4][5]. Nevertheless, our theoretical understanding of them was rather limited to the linear response regime [6] or non-interacting systems. Thus, the theoretical study of nonlinear transport in strongly interacting systems is highly challenging. More recently, nonlinear Drude weight (NLDW) characterizing the nonlinear static transport has been introduced [7,8]. This quantity is a straightforward extension of the linear Drude weight first proposed by Kohn [9] as an indicator to distinguish between a conductor and an insulator in quantum manybody systems. Given that the linear one has played an essential role in characterizing linear transport properties [10][11][12][13][14][15], we expect that its nonlinear counterparts will be equally or even more important in understanding transport phenomena.
There are already a number of previous studies focusing on the NLDWs [16][17][18][19]. In Ref. [16], the NLDWs in the spin-1/2 XXZ chain, which is a paradigmatic example of a quantum many-body system, was examined in detail. It was found that they diverge in certain anisotropy parameter regimes in the thermodynamic limit. In addition, the origin of these divergences was identified as nonanalytic finite-size corrections to the ground-state energy [16]. However, this property was discussed except when the anisotropy parameter takes special values, and thus there are still some cases that have not been thoroughly investigated. Therefore, further research on the NLDWs in this fundamental model is needed to achieve a complete understanding.
In this paper, we investigate the fine structure of the NLDWs at zero temperature for the spin-1/2 XXZ chain in the whole critical regime. The advantage of this model is its solvability by the Bethe ansatz [20,21]. Since the NLDWs can be read off from the higher-order finite-size corrections to the ground-state energy in the presence of a U (1) flux, it is essential to analyze these corrections in detail. This is achieved by using the Bethe ansatz, in conjunction with a mathematical method called the Wiener-Hopf method [16,20,[22][23][24][25]. Furthermore, since the Bethe ansatz enables us to treat very large systems numerically, we can confirm the asymptotic behaviors of the NLDWs in the large system-size limit. From the perspective of the Wiener-Hopf method, we reveal that the finite-size scaling of the ground state energy is quite distinct depending on the value of the anisotropy parameter. The two main findings of this study are as follows. The first one is the behaviors of the NLDWs at their boundaries between the convergent and divergent regions. The detailed analysis suggests that the nth order one there converges for n = 0, 1, 2 (mod 4), while it shows logarithmic divergence for n = 3 (mod 4) in the large system-size limit. By using the exact solutions, we calculate the first several orders of the NLDWs numerically and confirm their behaviors around the boundaries. The other one is the existence of particular anisotropies where all the NLDWs converge. Since higher order ones have wider divergent regions, some of the special anisotropies are surrounded by the divergent region. We confirm this discontinuous behavior in the critical regime by calculating one of the higher order NLDWs numerically.
Our paper is organized as follows: In Sec. II, we review the Bethe ansatz for the XXZ chain with the U (1) flux and introduce the nonlinear Kohn formula to calculate the NLDWs. In Sec. III, the main results of our study are summarized. In Sec. IV, we review the origin of the divergences of the NLDWs and carefully determine the convergent and divergent regions. In Sec. V, by consider-ing logarithmic corrections to the ground state energy, we analytically identify the behaviors of the NLDWs at their boundaries between the convergent and divergent regions and also confirm them numerically. In Sec. VI, we analytically and numerically reveal that there exist some exceptional points where all the NLDWs converge. Finally, the discussion and conclusion of our paper are presented in Sec. VII. In Appendices, we provide the derivation of the finite-size scaling of the ground-state energy based on the Wiener-Hopf method. Furthermore, numerical confirmation of the scaling for several anisotropies is also given there.

II. NONLINEAR DRUDE WEIGHTS IN THE XXZ CHAIN
We consider the spin-1/2 XXZ chain with the U (1) flux Φ defined by the Hamiltonian: whereŜ α l (α = x, y, z) are spin-1/2 operators,Ŝ ± l =Ŝ x l ± iŜ y l , J > 0 is the coupling constant, ∆ is the anisotropy parameter, and N is the number of sites. We impose periodic boundary conditions on the chain and assume that −1 < ∆ < 1 and N is even throughout this paper. See Fig. 1 for a schematic picture of our model. Here it is enough to consider only −π < Φ ≤ π, as H(Φ) and H(Φ+ 2π) have the same spectrum. We note in passing that the case with Φ = 0 corresponds to the spin-1/2 XXZ chain with the Dzyaloshinskii-Moriya (DM) interaction with a uniform DM vector along the z axis [26].
Since the total magnetizationŜ z tot = N l=1Ŝ z l is conserved in this model, we can obtain the lowest energy state in each sector individually by the Bethe ansatz [27]. In the sector with M down spins, the Bethe roots {v j (Φ)} are determined by the following Bethe equation for j =  = π(n − 1)/(n + 3) (n ∈ N). Note that we here denote a set of points where the NLDWs contain logarithmic corrections as S l . Furthermore, we can see that there exist exceptional points Se ≡ S   N (Θ) in the critical regime. The power-law divergent and logarithmic divergent behaviors in the large-N limit are shown in bold characters. Here we denote the boundary between the convergent and divergent regions as ∆ (n) B ≡ cos (π(n − 1)/(n + 3)). Note that O(x) and o(x) are Landau symbols indicating O(x)/x → (const.) and o(x)/x → 0 (N → ∞), respectively. order of response n (mod 4) 1, 2, . . . , M : where p n (v) ≡ 2 tan −1 tanh γ 2 v tan nγ 2 and γ ≡ arccos ∆. In terms of the Bethe roots, the energy density is given by If Φ = 0, it is known that the ground state lies in the sector of M = N/2 [28]. Thus, for sufficiently small Φ the ground state energy density of H(Φ) is e gs (Φ) = e(Φ; M = N/2). Nonlinear Drude weight (NLDW) is a straightforward extension of the linear Drude weight [9] and can be calculated by using the nonlinear Kohn formula [7,8]. At zero temperature, the nth order one D where −π < Θ ≤ π. Note that the finite Θ corresponds to the DM interaction as mentioned above. In the Θ = 0 case, only the odd orders are nonvanishing. This is because the ground-state energy density e gs (Φ) is an even function of Φ, which can be seen

III. OVERVIEW OF THE RESULTS
Here we summarize the main results of our paper. They are shown in Fig. 2 and Table I. As shown in Ref. [16], the nth order Drude weight D (n) N (Θ) has both the convergent and divergent regions. The boundary between them is given by In the following, we denote the convergent region as S and the divergent one as S ≡ {cos (πr/(r + 1)) | r = 1, 2, . . . , (n − 2)/4 } is the set of exceptional points at which the nth order one converges. Note that x is the floor function. The definitions of frequently used symbols are summarized in Table II. The results for the odd order NLDWs are shown in the first and second lines of Table I. In the convergent region S (n) c , i.e., −1 < ∆ < ∆ (n) B , they converge to finite values in the thermodynamic limit. At the boundary point ∆ = ∆ (n) B , they show two distinct behaviors depending on the order of the response. When n = 1 (mod 4), the NLDWs at their boundaries converge to finite values in the thermodynamic limit. On the other hand, when n = 3 (mod 4), the large-N asymptotic behavior of the NLDWs at their boundaries is the logarithmic divergence. In the divergent region S (n) d , the large-N asymptotic behavior of the NLDWs is the power-law divergence of the form D e , the odd order NLDWs converge to finite values unlike the divergent behaviors around the points. This behavior is the same as that in the convergent region S   e , the even order NLDWs vanish unlike the divergent behaviors around the points. This behavior is the same as that in the convergent region S (n) c . The above results are summarized visually in Fig. 2. It clearly shows that higher order NLDWs have the wider divergent regions S (n) d . One might think that the infiniteorder NLDW diverges everywhere in the critical regime. However, this is not the case because, at the exceptional points in S (∞) e = {cos (πr/(r + 1)) | r ∈ N}, all the NLDWs show the convergence as we have discussed above. The situation is illustrated in the bottom panel of Fig. 2.
In the following sections, we derive these results analytically by using the Wiener-Hopf method. Although there are some subtle points in this approach, we confirm our results by directly solving the Bethe ansatz equations numerically.

IV. THE CONVERGENT AND DIVERGENT REGIONS
In this section, we focus on the behaviors of the NLDWs in the convergent and divergent regions. These regions include points where e gs (Φ) contains logarithmic corrections [Eq. (9)]. Here we denote the set of these points by S l ≡ {cos (π(2p − 1)/(2p − 1 + 2q)) | p, q ∈ N}. For later convenience, we also introduce the set of exceptional points S e ≡ S (∞) e = {cos (πr/(r + 1)) | r ∈ N}, at which all the NLDWs converge. For ∆ ∈ (−1, 1) \ S e ∪ S l [29], the detailed analysis of their behaviors was given in our previous paper [16]. From the low-energy effective field theory of the XXZ chain or the Wiener-Hopf method, the finite-size scaling of e gs (Φ) can be cast into the form where ∆ ∈ (−1, 1)\ S e ∪ S l , and A k,l and B k,l,m are coefficients depending on γ (see Appendix A). Note that the smallest exponent of 1/N in the second sum of Eq. (6), namely 2 + 4γ/(π − γ), is always noninteger. In other words, there exist nonanalytic finite-size corrections to the ground-state energy. The straightforward differentiation of Eq. (6) with respect to Φ enables us to identify the large-N asymptotic behaviors of the NLDWs. They read where [30]. From these results, we can see that, in the thermodynamic limit, the odd order NLDWs converge to finite values and the even order ones vanish when n < 1 + 4γ/(π − γ), i.e., −1 < ∆ < ∆   ∩ S e , all the coefficients B k,l,m in Eq. (6) vanish identically as we will see in Sec. VI. This leads to the fact that the odd order NLDWs still converge to finite values and the even order ones still vanish in the thermodynamic limit. In the other case ∆ ∈ S (n) c ∪ S (n) d ∩ S l , the finite-size scaling of e gs (Φ) contains logarithmic corrections. However, even in these cases, when −1 < ∆ < ∆ (n) B , the odd order NLDWs still converge to finite values and the even order ones still vanish in the thermodynamic limit as we will see in Appendix B. Also, when ∆ π−γ in the large-N limit. As a result, we can conclude that every NLDW D         Table I) [31].

V. BOUNDARY BETWEEN THE CONVERGENT AND DIVERGENT REGIONS
As we have discussed in the previous section, the boundary point between the convergent and divergent regions of D = π(n − 1)/(n + 3). This suggests that, when n = 1 (mod 4), i.e., n = 4k + 1 (k ∈ N), the boundary ∆ (4k+1) B = cos (πk/(k + 1)) is included in the set of exceptional points S e . Since the special properties of the NLDWs at these points are discussed in Sec. VI, here we focus on the remaining cases: n = 0, 2, 3 (mod 4), i.e., n = 2k, 4k − 1 (k ∈ N).

The boundaries of D
correspond to the cases when p = k, q = 2 and p = k, q = 1, respectively. In these cases, the detailed analysis in Appendix A shows that the finite-size scaling of e gs (Φ) obeys where ∆ ∈ S l , and C k,l,s and D k,l,m,s are coefficients depending on γ (see Appendix A). From Eq. (9), the NLDWs at their boundaries can be calculated as where Note that the subscript "B" is introduced to indicate the value at the boundary. The above results mean that D (Θ) shows the logarithmic divergence in the large-N limit (see Table I). The analytical form of Y k (Θ) can be obtained from Eq. (4.1) in Ref. [32] as These behaviors can be confirmed numerically. The numerical results for D  (2) B = π/5, we believe that this difference is due to numerical errors in the finite-differentiation and the extrapolation process. We can also see that, in the large-N limit, the behavior of D  N (0) at the boundary fall almost on a straight line in the large-N region. Note that by calculating the quantities related to the derivative of the NLDWs with respect to N , we can avoid observing directly the logarithmic divergence of the NLDWs themselves, which is very difficult to identify numerically. As a result, we can confirm that D N (0) shows the logarithmic divergence at the boundary as we have expected from Eq. (11). This is because Eq. (11) yields and thus, for the case k = 1, we have N ∂ N D Here we have approximated the derivative with respect to N by finite differences. From Fig. 3 N,B (0)/J at 1/N = 0 can be estimated as −0.1414... by the linear extrapolation. On the other hand, the analytical expression for 4!C 1,2,1 can be obtained as −81 √ 3J/(32π 3 ) by setting γ = π/3 + and expanding Eq. (4.1) of Ref. [32] in around = 0. Thus we can indeed confirm that N ∂ N D N (0) above and below the boundary becomes diverging and vanishing in the large-N limit, respectively. The same holds for general k and can be understood from the following relation which follows from Eq. (7). Since the power of the second term in Eq. (14) is the same as one appearing in Eq. (7) for D (4k−1) N (Θ), the above quantity in the convergent region S (4k−1) c vanishes in the thermodynamic limit. On the other hand, the large-N asymptotic behavior of the above quantity is the power-law divergence of the form

VI. EXCEPTIONAL POINTS
Now we focus on the behavior of the NLDWs at the exceptional points S e = {cos (πr/(r + 1)) | r ∈ N}. These points have a special property that all the coefficients B k,l,m in Eq. (6) vanish identically. This can be derived analytically by using the Wiener-Hopf method (see Appendix A). As a result, the finite-size scaling of e gs (Φ) can be written as for ∆ ∈ S e . For example, we can obtain the exact form of e gs (Φ) at the free-fermion point (∆ = 0) [7,33,34] as The large-N expansion of Eq. (16) consists of the terms Φ α /N β with α, β integers and α ≤ β. Thus, nonanalytic finite-size corrections do not appear in the expansion. This is obviously consistent with Eq. (15). Interestingly, in the thermodynamic limit, all the NLDWs converge at any ∆ in S e (see Fig. 4 (a)). This can be seen by noting that which means that the odd order NLDWs remain finite, while the even order ones vanish in the thermodynamic limit (see Table I). These convergent behaviors of all the NLDWs are consistent with the prediction based on numerical studies of small systems [35]. We can confirm these behaviors by numerically solving the Bethe ansatz equations for large system sizes. The exceptional points where all the NLDWs converge are shown in Fig. 4(a). There are infinitely many such points and they accumulate at the ferromagnetic point ∆ = −1. The numerical result for D (11) N (0) is shown in Fig. 4(b). As we can see in Fig. 4(b), the two points ∆ = 0, −0.5 included in S N (0) at these points converge unlike the divergent behaviors around the points. Since the diver-gent behaviors there should be caused by noninteger power terms of N in Eq. (7), we calculated r (11) (1/N ) ≡ D (11) N (0)/(12! B 1,6,1 N 10−4γ/(π−γ) ) numerically. The result for this quantity is shown in Fig. 4(c). This figure clearly shows that each data is on a straight line to the value near 1 in the large N region as we have expected.
Some remarks are in order. First, the spin-1/2 XXZ chain with periodic boundary conditions has a special symmetry related to the sl 2 loop algebra [36,37] at the exceptional points S e [38]. We speculate that this symmetry is responsible for the convergence of all the NLDWs in the thermodynamic limit. Second, the coefficient of the umklapp scattering term (the cosine term) in the low energy effective Hamiltonian of the XXZ chain vanishes at these points (see Eq. (2.23) in [32]). Considering that this term leads to the nonanalytic finitesize corrections and gives the leading contribution to the power-law divergence, we can see that its vanishing is consistent with the convergence of all the NLDWs. Finally, the ground state energy e gs (Φ) has the peculiar adiabatic period at these points. For this case the adiabatic period of e gs (Φ) is of the order of the system size N , while for the other cases the period is 4π [39]. Based on this property, numerical calculations for small system sizes have recently revealed that the current density exhibits nontrivial oscillations, so-called Bloch oscillations, at the points in S e even under an infinitesimal external field [35].

VII. DISCUSSION AND CONCLUSION
In this paper, we examined the fine structure of the NLDWs at zero temperature for the spin-1/2 XXZ chain in the critical regime (see Fig. 2). In order to calculate the NLDWs, we investigated the finite-size corrections to the ground-state energy of the chain with U (1) flux and revealed that its finite-size scaling was quite distinct depending on the anisotropy parameter ∆. Based on the expansions Eqs. (6), (9) and (15), we studied the largesize asymptotic behavior of the NLDWs both analytically and numerically. The analysis determined the convergent and divergent regions of the NLDWs, the boundary of which depends on the order of the response n. We studied the behaviors of the NLDWs at the boundaries in detail and found that they converge for n = 0, 1, 2 (mod 4), while they show the logarithmic divergence for n = 3 (mod 4) in the large system-size limit (see Table I). In addition, we numerically confirmed not only the convergence but also the logarithmic divergence at the boundaries of the first several orders of the NLDWs (see Fig. 3). Furthermore, we revealed that there exist special values of ∆ where all the NLDWs converge in the thermodynamic limit. Since higher order ones have wider divergent regions, some of the special ∆ are surrounded by the divergent region. We confirmed this discontinuous behavior in the critical regime by calculating one of the higher order NLDWs numerically (see Fig. 4).
In order to obtain the finite-size scaling of the groundstate energy, we employed the Wiener-Hopf method for the finite-size system, which is based on the Euler-Maclaurin formula [40,41]. Traditionally, when calculating the leading finite-size corrections to the ground state energy, higher order terms included in the expansion by this formula are often ignored [23,42]. In general, there is no guarantee that these terms are negligible to calculate the corrections in other problems [43,44]. Thus, in our study, we took all these higher order terms into account and obtained the higher order corrections to e gs (Φ) as well as the leading ones. Here we should note that, although this enables us to overcome the above problem, we cannot determine the coefficients of these corrections in closed form within this approach. Also, we assume that e gs (Φ) can be Taylor-expanded around Φ = 0 based on the symmetry of the model and comparison with the analytical results in the thermodynamic limit [16]. Therefore, although we have confirmed our results numerically for several ∆, a more rigorous derivation of the results using another method is desirable and would be an interesting future direction.
Finally, we discuss the implications of our results to the transport phenomena. One might think that the divergent behaviors of NLDWs imply the divergence of a total current density. However, this seems unlikely because contributions to the current density from different orders can cancel each other out. In fact, a similar situation is observed in a single-band tight-binding chain with a defect [18]. Although the NLDWs of this system generally diverge with system size, real-time numerical simulation suggests that the adiabatic current density is suppressed compared to the defect-free case, in which the NLDW remains finite at any order. Although these finite-size corrections based on the same method had been partly discussed in Ref. [23,42], here we expose the mathematical details and illustrate the derivation process for readers' convenience. This detailed analysis also enables us to reveal that there are some cases with logarithmic finite-size corrections. As a result, we derive the general expression of the finite-size corrections including logarithmic ones. Furthermore, after calculating these finite-size corrections, we introduce the U (1) flux into them and obtain the finite-size scaling of e gs (Φ).

A-1. Setup
First, we review the Bethe ansatz and derive some important relations. It is known that the ground-state energy of the above model can be obtained by this ansatz. The Bethe roots {v j } are determined by the following Bethe equations for j = 1, 2, . . . , N/2: where (A4) Note that there exists a unique set of real solutions {v j } satisfying −∞ ≤ v 1 < v 2 < . . . < v N/2 ≤ ∞ and v j = −v N/2−j+1 . Differentiating Eq. (A3) with respect to v, we get where a n (v) ≡ 1 2π Then {v j } gives the ground-state energy density as where A = 2J sin γ/γ. Now we introduce a new useful function S N (v) as This transforms Eq. (A5) into the following form: Here we define a Fourier transformation of a function f (x) asf By using Fourier transformation on both sides of Eq. (A11), we get where the Fourier transform of a n (v) is Then by using Fourier transformation on both sides of Eq. (A13), we obtain where ρ ∞ (v) and R(v) are defined as follows: Note that ρ ∞ (v) is the exact representation of ρ N (v) in the thermodynamic limit. Similarly, Eq. (A9) leads to where we introduced Note that the third line follows from Eq. (A16) and the last from the following relation: Since we can see that only the second term in Eq. (A21) is responsible for the finite-size corrections to the groundstate energy, we only have to evaluate the effect of S N (v) to achieve the goal. Next, we introduce a useful formula to treat S N (v) included in the integral. The derivation of the formula is based on the Euler-Maclaurin formula [40,41]: where f (x) is a continuous function, x is the floor function, and B k (x) is the kth Bernoulli polynomial satisfying By using the recurrence relation (A26) and integral by parts, we naively obtain where B k = B k (0) is the kth Bernoulli number. Although we have B 2l+1 = 0 (l ∈ N), we keep these terms explicit in the following discussion. The above relation and the fact that I j+1 − I j = 1 give us the following relation: It is obvious that the above relation enables us to evaluate the finite-size corrections in Eqs. (A16) and (A21).
Finally, we introduce important relations employed in the Wiener-Hopf method briefly. In the following discussion, we denote v N/2 (= −v 1 ) as Λ. By using Eq. (A33), we get These are the complete representations of the finite-size corrections using N and Λ. Thus, in order to obtain the corrections using only N , we have to derive the relation between N and Λ. (Actually, we can roughly identify e −(π/2)Λ with 1/N as we will see in the following discussion.) Now we introduce new functions where Θ(v) is a Heaviside step function. Then by substituting v + Λ to the argument of Eq. (A34), we have where we introduced coefficients P k and Q k depending on 1/N and ρ (n) N (Λ) = g (n) (0) for n ≥ 0, and superscripts denote numbers of derivatives. Note that Eq. (A34) suggests that all the terms included in P k or Q k can be expressed as follows: where l, m, N n ∈ Z ≥0 , and each power satisfies n N n = m − l and n nN n = l − k. Since we have ρ Here we investigate behaviors of ρ ∞ (v + Λ) and R(v + 2Λ) for v > 0, which appear in Eq.(A38). Since Eqs. (A17) and (A19) give we can see that poles ofρ ∞ (ω) orR(ω) in the lower-half plane contribute to ρ ∞ (v+Λ) and R(v+2Λ), respectively. The position of the poles can be read off from the explicit expressions forρ ∞ (ω) andR(ω) as follows: where p, q ∈ N. SinceR(ω) have poles dependent on the parameter γ, in order to obtain the finite-size corrections, we must consider whether all the poles ofR(ω) are distinct or not. Thus, we perform the following classification shown in Table III. Actually, this classification is essential to ensure convergence of coefficients A k,l , B k,l,m , C k,l,s and D k,l,m,s appearing in Eqs. (A91),(A95) and (A103).
Here we consider the case of γ = π(2p − 1)/(2p − 1 + 2q) or πr/(r+1) (p, q, r ∈ N). In this case, all the poles ofρ ∞ (ω) R (ω) are distinct simple poles. Thus, we have for v > 0. Here we denoted a residue of a function f (x) at x = x 0 as Res(f, x 0 ). It is obvious that poles closer to the real axis contribute to the smaller power of e −(π/2)Λ . Therefore Eq. (A38) implies that g(v) can also be expanded as g(v) = g [1] (v) + g [2] (v) + · · · , (A48) where superscripts denote increasing powers of e −(π/2)Λ or 1/N . Then by substituting Eq. (A48) into Eq. (A38) and extracting the same order terms, we obtain, for examples, [1] , (A49) g [2] (v) − ρ ∞ (v + Λ) + (u)du [2] − R(v + 2Λ) 2N [2] (A50) where superscripts [· · · ] [n] again denote increasing powers of e −(π/2)Λ or 1/N . By using Fourier transformation and integral by parts, we get g [1] + (ω) +g (A51) g [2] + (ω) +g [2] − (ω) − ρ ∞ (ω)e −iωΛ [2] =R(ω)g [2] + (ω) + R (ω)g [1] . (A52) From the above relations, we can obtain the orders of g [1] + (ω) andg [2] + (ω) by splitting the whole into two parts: the part analytic in the upper half-plane and the other part analytic in the lower half-plane. Now we recall that a Fourier transformf (ω) can be split as follows: wheref + (ω) andf − (ω) are defined as and are analytic in the upper and lower half-plane, respectively (actually,g [n] ± (ω) are examples). We also introduce the following convenient factorization [16,20,[22][23][24]: where G + (ω) and G − (ω) are written as and are analytic and non-zero in the upper and lower half-plane, respectively. Here we calculateg [1] + (ω) as an example. By using the above methods for splitting, we can transform Eq. (A51) as g [1] We see that the left-and right-hand side of Eq. (A59) are analytic in the upper and lower half-plane, respectively. Since both of them are analytic on the real axis, the right-hand side of Eq. (A59) is the analytic continuation of the left-hand side, and thus there should be the entirely analytic form P [1] (ω) [23,25]. Although the form of P [1] (ω) is determined so thatg [1] + (ω) → 0 (|ω| → ∞), we do not need the explicit form for our purposes. As a result, we obtaiñ g [1] where c 1 , c 2 and c 3 (∼ O(1)) are certain coefficients. The term with e − π 2 Λ originally derives from the term which is contributed by the simple pole ofρ ∞ (ω) closest to the real axis: ω = −iπ/2. The higher order ones can be calculated in the same way. However, since in that case the poles ofR(ω) contribute tog [n] for example,g + (ω) can contain special power terms like (e − π 2 Λ ) 4γ/(π−γ) . Therefore we havẽ g + (ω) =g [1] + (ω) +g [2] + (ω) +g [3] where A k,l,m and B(∼ O(1)) are certain coefficients. Now we can obtain the relation between Λ and N . Recalling Eqs. (A2) and (A3) we get Then by substituting Eq. (A71) into its right-hand side successively, we obtain where C k,l and D k,l (∼ O(1)) are certain coefficients. Note that all the terms included in D k,l are expressed as follows: where s, m, N n ∈ Z ≥0 , and each power satisfies n N n = m − s. Here since we have where we introduced coefficients E k,l,m and F as the integral values of A k,l,m and B, respectively. Thus by substituting Eqs. (A72) and (A77) into the right hand side of Eq. (A77) successively and using the Maclaurin expansion with respect to e − π 2 Λ ( 1), we obtain ρ (n) where H n,k,l is a certain coefficient. Since Eqs. (A66), (A72) and (A78) give us we finally get the following relation between N and Λ from Eq. (A70): where I k,l (ω) is a certain coefficient. Then the sequential substitution of its right-hand side into e − π 2 Λ yields where J k,l is a certain coefficient. Now we can express e gs,N by using only N . In order to obtain this expression, we use the following relations: where K m,l is a certain coefficient. Therefore Eqs. (A35) and (A82) yield where L k,m is a certain coefficient. This is the finite-size corrections to the ground-state energy density of H(0). By introducing the effect of U (1) flux Φ into the coefficients in Eq. (A88), we get e gs,N (Φ) = e gs,∞ + k≥1,m≥0 Note that e gs,∞ is independent on Φ because this is the value in the thermodynamic limit. Then since the inversion symmetry of the model guarantees that e gs,∞ is an even function of Φ, we naturally expect that the difference e gs,N (Φ) − e gs,N (0) obeys the following finite-size scaling: where we introduced certain coefficients A k,l and B k,l,m . The above expression is intentionally split into two parts based on whether terms are including contribution from the poles ofR(ω) dependent on the parameter γ, namely ω = −iqπγ/(π − γ) (q ∈ N), or not. Although we cannot find any constraint on the summations within this analysis, comparison of the results for the NLDWs calculated from Eq. (A90) in the thermodynamic limit and the analytical ones [16] yields where a new constraint on the summation appears in the first term.

A-3. Case (ii): ∆ ∈ Se
Here we consider the case of γ = πr/(r + 1) (r ∈ N). Since all the values of (A43) are still distinct in this case, the same analysis as we have seen in case (i) is applicable. However, there is a significant difference between case (i) and (ii), i.e., the asymptotic behavior of R(v+2Λ). When γ = πr/(r + 1), the residues ofR(ω) π−γ 2 cos qπγ π−γ · π γ − 1 cos (qπ) (A92) = i sin (r − 1)qπ 2 cos (rqπ) · π γ − 1 cos (qπ) (A93) = 0, which means they are no longer poles ofR(ω). The above fact suggests the second term in Eq. (A46) vanishes, and thus we have for v > 0. Since this results in vanishing of all the terms related to the poles ω = −iqπγ/(π − γ) (q ∈ N), all the coefficients B k,l,m vanish in Eq. (A91). Thus we naturally obtain the following finite-size scaling: This can be understood from the perspective of the c = 1 conformal field theory perturbed by irrelevant operators. For examples, the above discussion is consistent with the fact that the coefficient of the umklapp term (the cosine term) vanishes at these points (see Eq. (2.23) in [32]).
Here we consider the case of γ = π(2p − 1)/(2p − 1 + 2q) (p, q ∈ N). Unlike the other cases we have seen so far, some of (A43) take the same values, which can be written as for l ∈ N. Since this means thatR(ω) have double poles at the above points, we have to use the following asymptotic expansion for v > 0 instead of Eq. (A46): where represent the summation over simple poles, namely poles excluding ω = ω l . The most different point from the other cases is the third term in Eq. (A97). For example, we can see which reveals that a term proportional to Λe −π(2p−1)Λ newly appears in R(v + 2Λ). This fact suggests that the previous relations (A80), (A82) and (A88) are modified as follows: where M k,l,m , N k,l,m and O k,m,s are certain coefficients. Therefore by referring to Eq. (A91), we naturally expect that the difference e gs,N (Φ)−e gs,N (0) obeys the following finite-size scaling: where we introduced certain coefficients C k,l,s and D k,l,m,s .
Appendix B: The detailed analysis of the case (iii): Here we examine the NLDWs in the case (iii). From Eq. (A103), they can be calculated as where Y k (Θ) ≡ l>k (2l)!/(2l − 2k − 1)!D 1,l,1,0 Θ 2l−2k−1 and χ[E] takes the value 1 if E is true and 0 otherwise. The greatest benefit of this analysis is that Eqs. (B1) and (B2) enable us to identify the large-N asymptotic behavior of the NLDWs at their boundaries between the convergent and divergent regions (see the main text). Now we consider the behavior of the NLDWs at the points of case (iii) other than the boundary, namely ∆ ∈ S (n) c ∪ S (n) d ∩ S l . Since Eqs. (B1) and (B2) suggest that the effect of the logarithmic correction can appear in the NLDWs, we investigate their effect in detail. Here we express the region with logarithmic corrections to D log can be evaluated as follows: where γ (n) B = arccos ∆ (n) B and x is the floor function. The equal sign in Eq. (B4) holds when k is even. Note that γ (2k−1) log and γ (2k) log correspond to the cases when p = k/2 , q = 1 and p = (k + 1)/2 , q = 1, respectively, as γ = π(2p − 1)/(2p − 1 + 2q) is monotonically increasing for p and decreasing for q.
Based on the above results, we can specify the large-N asymptotic behaviours of D , the absence of the logarithmic corrections reproduces the same behaviour as Eq. (7), and thus D (2k−1) N (Θ) shows the convergence to a finite value in the thermodynamic limit (see Table I). Similarly, although there exist the logarithmic corrections, D   Table I).
Appendix C: Numerical results for the finite-size corrections We show some numerical results for the finite-size corrections with the U (1) flux. Since the expansion (A90) was studied in our previous paper [16], here we focus only on the new result (A103). In order to confirm the finite-size scaling of e gs,N (Φ), it is better to calculate the NLDWs instead of e gs,N (Φ) itself. Below we discuss several dominant terms of the NLDWs in certain cases.
We show the results for three examples: the secondorder, fifth-order and seventh-order ones.
First, we consider the fifth-order and seventh-order ones, which diverge in the thermodynamic limit at this point. From Eq. (B1), they can be calculated as On the other hand, in order to obtain their sub-leading behaviors, it is useful to differentiate them with respect to N . Therefore we introduce the following quantities: Finally, we consider the second-order one, which converge in the thermodynamic limit at this point. From Eq. (B2), this can be calculated as where we obtained C 1,2,1 at γ = π/3 from Eq. (4.1) in Ref. [32]. Then, in order to evaluate the leading behavior, we introduce the following quantity: The numerical results for the above quantity is shown in Fig.5(a-3). This figure clearly shows that, in the large-N region, the data fall on a straight line to the analytical value −81 √ 3/(32π 3 )Θ = −0.01414... indicated by the orange dotted line, which confirms Eq. (C7). C-2. At γ = π 5 (p = 1, q = 2) We consider the fifth-order one. From Eq. (B1), this can be calculated as As in the previous case, we introduce the following quantities: The numerical results for the above quantities are shown in Fig. 5(b). Since the analytical form of D 1,k,1,0 can be obtained as from Eq. (4.1) in Ref. [32], we have 6!D 1,3,1,0 /J = −0.02881... at γ = π/5 which is indicated by the orange dotted line. This figure obviously shows that the data fall on straight lines to finite values in the large-N region, which confirms Eq. (C10).
C-3. At γ = 3π 7 (p = 2, q = 2) We consider the seventh-order one. Here we focus not only on the leading and sub-leading term but also on the sub-sub-leading term. From Eq. (B1), this can be calculated as D As in the previous case, we introduce the following quantities: