Key Observable for Linear Thermalization

For studies on thermalization of an isolated quantum many-body system, the fundamental issue is to determine whether a given system thermalizes or not. However, most studies tested only a small number of observables, and it was unclear whether other observables thermalize. Here, we study whether `linear thermalization' occurs for all additive observables: We consider a quantum many-body system prepared in an equilibrium state and its unitary time evolution induced by a small change $\Delta f$ of a physical parameter $f$ of the Hamiltonian, and examine whether \emph{all} additive observables relax to the equilibrium values in a manner fully consistent with thermodynamics up to the linear order in $\Delta f$. We find that the additive observable conjugate to $f$ is key for linear thermalization in that its linear thermalization guarantees, under physically reasonable conditions, linear thermalization of all additive observables. Such a linear thermalization occurs in the timescale of $\mathcal{O}(|\Delta f|^0)$, and lasts at least for a period of $o(1/\sqrt{|\Delta f|})$. We also consider linear thermalization against the change of other parameters, and find that linear thermalization of the key observable against $\Delta f$ guarantees its linear thermalization against small changes of any other parameters. Furthermore, we discuss the generalized susceptibilities for cross responses and their consistency between quantum mechanics and thermodynamics. We demonstrate our main result by performing numerical calculations for spin models. The present paper offers an efficient way of judging linear thermalization because it guarantees that examination of the single key observable is sufficient.


I. INTRODUCTION
Thermalization of isolated quantum many-body systems has long been studied as the question of how the thermal equilibrium state arises from the quantum unitary dynamics .
Many studies have also been devoted to the eigenstate thermalization hypothesis (ETH) [1,4,5], which states that every energy eigenstate represents an equilibrium state. This hypothesis leads to thermalization, and is expected to hold in quantum chaotic systems [11-14, 16, 17, 23]. However, since there also exist many systems that do not satisfy the ETH, how universally it holds and when it fails are still the subjects of active research [24][25][26][27][28][29][30][31][32][33]. Furthermore, while the ETH is a sufficient condition for thermalization, whether it is also a necessary condition seems to be controversial because the answer depends on the choices of the initial state and observables that are employed to check thermalization [22,24,25].
In such nonthermalizing systems, the entanglement entropy of energy eigenstate often behaves anomalously [108,118,119,125,126], and hence is sometimes employed to discriminate between thermal and nonthermal eigenstates [108,[127][128][129][130]. However, even if the entanglement entropy agrees with thermodynamic entropy, it does not necessarily imply that the eigenstate represents the equilibrium state because there exist many quantum states with the same entanglement entropy. For this reason, many studies examined the expectation values of observables to discriminate between thermal and nonthermal states [18,42,47,[131][132][133][134].
When testing thermalization using observables, most of previous studies examined only a small number of observables [11-14, 16, 18, 22, 42, 47, 131-134]. However, thermalization of such observables does not necessarily guarantee thermalization of other observables. For example, the authors showed in Fig. 2(b) and (c) of Ref. [135] that, in the XXZ and the XY spin chains, the magnetization of any nonzero wavenumber (such as the staggered magnetization) thermalizes while the uniform (zero wavenumber) magnetization does not. A simpler example is the case where thermalization trivially occurs by symmetry. For instance, when both the Hamiltonian and the initial state are symmetric under spin rotation by π around the z axis, the x component of magnetization is kept 0 under the unitary time evolution (because of the symmetry) and hence thermalizes trivially. Nevertheless, if this system is integrable many other observables such as interactions between nearest neighbor spins do not thermalize, as in the case of Model III of this paper. Another example, which seems nontrivial and most interesting, is the possibility that all one-body observables thermalize whereas two-and more-body observables do not. All these examples show that thermalization of some observables does not necessarily imply thermalization of all observables.
These issues have also been long-standing questions even in the linear nonequilibrium regime. For example, Kubo stated in his famous paper on the 'Kubo formula' [136] that if the system has an 'ergodic property' then his formula would give the isothermal response. However, his discussion turned out wrong [137][138][139]: Later studies proved that the Kubo formula gives the adiabatic response under certain conditions [135,[137][138][139], and it can also give the isothermal response under another condition by taking the limit of vanishing wavenumber [135]. However, the conditions given by these studies were for individual observables, which do not guarantee the conditions for all observables. In other words, the above-mentioned fundamental issue of thermalization has left unsolved even in the linear nonequilibrium regime.
In this paper, we study the unitary time evolution induced by a small change (quench) ∆f of a physical parameter f , and examine whether the quantum state relaxes to an equilibrium state that is fully consistent with thermodynamics up to O(∆f ). We call the relaxation in this sense linear thermalization, which is obviously necessary for thermalization against an arbitrary magnitude of ∆f . We place a particular emphasis on the full consistency. That is, when comparing the quantum state with a thermal equilibrium state we examine all additive observables because thermodynamics assumes that all additive observables take macroscopically definite values in an equilibrium state [140][141][142], although the state is specified by only a small number of variables (such as temperature). Furthermore, we consider the case where the initial state is an equilibrium state because thermodynamics basically treats transitions between equilibrium states.
We find that the additive observableB which is conjugate to f is the key observable in linear thermalization. We show rigorously that, if its expectation value relaxes to the value predicted by thermodynamics, so do the expectation values of all additive observables. Even when B is a simple one-body observable, its relaxation guarantees relaxation of all other additive observables including two-and more-body ones. In addition to this theorem, we prove two propositions which state that, under reasonable conditions, the time fluctuations of the expectation values and the variances of all additive observables are sufficiently small. These results mean that linear thermalization of the single observableB guarantees linear thermalization of all additive observables. This linear thermalization occurs in timescale of O(|∆f | 0 ). We prove that it lasts at least for a period of o(1/ |∆f |). Furthermore, we show that linear thermalization ofB against the quench of ∆f implies linear thermalization ofB against the quench of any other parameters. Moreover, for the generalized susceptibilities (crossed susceptibilities) our theorem gives a necessary and sufficient condition for the consistency between quantum mechanics and thermodynamics. As demonstrations of the theorem, we present numerical results for three models of spin systems.
These results dramatically reduce the costs of experiments because a single quench experiment on the key observable in a timescale of O(|∆f | 0 ) gives rich information about all additive observables, a longer time scale o(1/ |∆f |), and the quench of other parameters.
The paper is organized as follows. Section II explains the setup. Section III defines linear thermalization by introducing three criteria. The theorem (main result) and two propositions are summarized in Sec. IV. We discuss the timescale of linear thermalization in Sec. V, and prove the third proposition about a longer timescale. Generalized susceptibilities are analyzed in Sec. VI, where two corollaries are presented. Numerical demonstrations are given in Sec. VII. We prove the theorem and propositions in Sec. VIII. In Sec. IX, we show that our results are also applicable to the case of continuous change of f , and explain a relation between our main result and the ETH. Section X summarizes the paper.

II. SETUP
We study an isolated quantum many-body system that is defined on a lattice with N sites and described by a finite dimensional Hilbert space [143]. The system obeys the unitary time evolution generated by the Hamiltonian H(f ), which depends on a physical parameter f [144] such as an external magnetic field. We use units where the reduced Planck constant ℏ and the Boltzmann constant k B are unity. We investigate the system prepared in an equilibrium state and its time evolution induced by a change of f . We specifically consider a quench process, in which f is changed discontinuously. This does not mean any loss of generality because the same results can also be obtained for continuous change of f , as shown in Sec. IX B.
Since we are interested in the consistency with thermodynamics, we consider the case where each equilibrium state of the system is uniquely specified macroscopically by an appropriate set of variables. Here, the number of the variables in the set is finite and independent of N . (This is one of the basic assumptions of thermodynamics [145].) To be specific, we here assume that f , N , and the inverse temperature β is such a set of variables. Then, equilibrium states can be represented by the canonical Gibbs statê where Z(β, f ) is the partition function. For time t < 0, f takes a constant value f 0 , which defines the initial Hamiltonian, and the system is in an equilibrium state at a finite inverse temperature β 0 , represented bŷ At t = 0, f is changed from f 0 to another constant value f 0 +∆f discontinuously, as shown in Fig. 1(a). Due to this quench of f , the state for t > 0 evolves according to the Schrödinger equation described by the postquench HamiltonianĤ Consequently the expectation values of additive observables evolve in time, as shown in Fig. 1(b). Here, we say an observableÂ is additive when it is the sum of local observablesâ r over the whole system, where we say an observableâ r is local when its support consists of sites within a distance of at most O(N 0 ) from the site r.
Note that a local observableâ r is not necessarily a onebody observable. It can be, say, a two-spin observable such as the one given by Eq. (41) below. Note also that the additive observableÂ can be noninvariant under a spatial translation. For example, for a magnetic field of a wavenumber k with magnitude f , its interaction with spins, is an additive observable withâ r = −h(r)σ z r . This example also shows that our theory is applicable to the case where a spatially-varying external field is applied. Since we allow such an r-dependent coefficient inâ r , we exclude strange cases where the operator norm ∥â r ∥ ∞ diverges as |r| → ∞ by imposing where C A is a constant independent of N and r. From this restriction, any additive observableÂ satisfies We also make a natural assumption thatĤ(f ) and its derivative,B are additive observables. Indeed this assumption is satisfied for the above example ofĤ(f ) [see also Eq. (32) below]. We say that the observableB is conjugate to the parameter f , and vice versa. We examine whether the system relaxes to the equilibrium state for small |∆f |. For this purpose, we consider the case where an equilibrium state exists not only for f = f 0 but also for f = f 0 + ∆f . That is, we exclude, for instance, electric conductors in a uniform electric field. Furthermore, to exclude uninteresting divergences, we also limit ourselves to the case where no phase transition occurs in the neighborhood of the initial equilibrium state. This implies that the initial equilibrium state is stable against the changes of parameters conjugate to any additive observables. (This condition is precisely expressed by Eq. (53) in Sec. VIII A.)

III. CRITERIA FOR LINEAR THERMALIZATION
Since we consider the case of small |∆f |, we examine the consistency with thermodynamics up to the linear order in ∆f : We say linear thermalization occurs for an additive observableÂ if the following three criteria are satisfied up to O(∆f ).
Criterion (i) (consistency of expectation value): The long time average of the expectation value ofÂ after the quench is sufficiently close to the equilibrium value. More precisely, for sufficiently long time T . (See Sec. V for detailed discussions on T .) Herê is the Heisenberg operator ofÂ evolved byĤ ∆f , and, for any t-dependent quantity g(t), Furthermore, ⟨•⟩ eq 0 := Tr ρ eq 0 • is the expectation value in the initial stateρ eq 0 , while ⟨•⟩ eq ∆f := Tr ρ eq ∆f • is the expectation value in The latter stateρ eq ∆f represents the final equilibrium state predicted by thermodynamics. Its inverse temperature β ∆f is determined from energy conservation Note that β ∆f ̸ = β 0 even in O(∆f ), as explicitly given by Eq. (A3) in Appendix A 1. Criterion (ii) (equilibration): At almost all t > 0, the expectation value ofÂ is sufficiently close to its long time average. That is, their difference is macroscopically negligible in the sense that for sufficiently long time T . Criterion (iii) (smallness of variance): At almost all t > 0, the variance ofÂ is macroscopically negligible such that for sufficiently long time T . Here Var 0 [•] := ⟨(•) 2 ⟩ eq 0 − (⟨•⟩ eq 0 ) 2 . Thermodynamics assumes that all additive observables take macroscopically definite values in an equilibrium state [140][141][142]. Since we are interested in full consistency with thermodynamics, we examine whether linear thermalization occurs, i.e., whether these criteria are satisfied, for all additive observables.

IV. SUMMARY OF RESULTS
Among three criteria of the previous section, Criterion (i) has been studied most intensively because in many cases it discriminates thermalizing systems from nonthermalizing ones. For instance, it was observed in many integrable systems that relaxation to some steady states occurs, indicating that Criteria (ii) and (iii) are satisfied, whereas they are nonthermal states, which do not satisfy Criterion (i) [41,43,44,[75][76][77][78].
Therefore we place a theorem about Criterion (i) as our main result. We also obtain additional results which show that Criteria (ii) and (iii) are easily satisfied in our setting, under conditions weaker than those of previous studies [7,146,147]. These results are summarized in this section. We will extend them slightly in Sec. V, and thereby discuss the timescale of the linear thermalization.

A. Theorem for Criterion (i)
Our main result is that if the additive observableB conjugate to f , given by Eq. (10), satisfies Criterion (i) of Sec. II up to O(∆f ), then so do all additive observables. To be precise, we obtain Theorem (main result): holds for every additive observableÂ if and only if it holds forÂ =B, This theorem identifiesB, which is conjugate to f , as the key observable for linear thermalization of all additive observables. That is, one can judge whether the long time average of the expectation value of every additive observableÂ is consistent with thermodynamics by examining the single observableB. This striking result has the following significant features.
Firstly, it is in sharp contrast with the existing approaches to Criterion (i) [11][12][13][14][15]148], such as the ETH. In fact, previous studies clarified that the ETH for an observable implies Criterion (i) for that observable. However, it did not guarantee either the ETH or Criterion (i) for other observables. By contrast, our Theorem does guarantee that Criterion (i) is satisfied by all additive observables up to O(∆f ) if it is satisfied byB. (We will demonstrate this point numerically in Sec. VII.) This result may be most surprising in the case whereB is a one-body observable: If Eq. (19) holds for that one-body observable, then it also holds for all other additive observables including two-and more-body (O(N 0 )-body) ones.
Secondly, the ETH forB implies condition (19) but the converse is not necessarily true, as will be discussed in Sec. IX C.
Thirdly, our Theorem has significant meanings about the generalized susceptibilities for cross responses, as will be discussed in Sec. VI.
Fourthly, we can extend the result such thatÂ and B in Eqs. (18) and (19) are not restricted to additive observables, as will be discussed in Sec. IX A.
for every additive observableÂ. [A detailed expression of the right-hand side (r.h.s) will be given in Sec. VIII B.] This proposition is obtained by adapting the argument of Short and Farrelly [146] to our setting. While their result contains the effective dimension of the initial state, our condition (20) instead contains its purity Tr (ρ eq 0 ) 2 . Since the postquench HamiltonianĤ ∆f is different from the initial HamiltonianĤ 0 , evaluation of the effective dimension of the initial stateρ eq 0 in terms ofĤ ∆f is not an easy task. By contrast, its purity can be evaluated easily as shown in Appendix A 2.
Note that our condition is weaker than the "nonresonance condition," D res = 1, of the previous studies [6,149,150]. The nonresonance condition often fails even in nonintegrable systems, e.g., when energy eigenvalues have degeneracies due to symmetries. By contrast, condition (20) is expected to hold not only in such systems but also in wider classes of systems, including interacting integrable systems (such as Model III of Sec. VII), because at any nonzero temperature 1/β 0 as shown in Appendix A 2.

C. Proposition for Criterion (iii)
Under the reasonable condition that the fluctuation of every additive observable in the initial equilibrium state is sufficiently small, the variances of all additive observables after the quench remain small enough such that Criterion (iii) is satisfied up to O(∆f ). To be precise, we find Proposition 2: If the fourth order central moment of every additive observableÂ in the initial state satisfies for every additive observableÂ. [A detailed expression of the r.h.s will be given as Eq. (62) in Sec. VIII C.] Note that condition (23) also bounds the variance in the initial state as which follows from the Cauchy-Schwarz inequality. By inserting this into Eq. (24), we have [151] Var (26) which shows that Criterion (iii) is satisfied up to O(∆f ).
Since we exclude phase transition points, our condition (23) seems plausible. By contrast, previous results on Criterion (iii) [59,147] required a condition that is stronger than the ETH. The above proposition shows that the criterion is satisfied under a much weaker condition in our setting.
To sum up these theorem and propositions, linear thermalization of the single key observableB guarantees linear thermalization of all additive observables, under physically reasonable conditions.

V. TIMESCALE OF LINEAR THERMALIZATION
When studying thermalization theoretically, it is customary to take the limit T → ∞ [59,61], as we did above. However, more detailed information about the timescale is necessary because in experiments thermalization occurs in reasonably short timescales [42,47]. The timescale is particularly nontrivial when "prethermalization" [61,[152][153][154][155][156][157] occurs, i.e., when the system first relaxes to a nonthermal quasi-steady state, and at some time stage after that, it relaxes to a true thermal equilibrium state. In nearly integrable systems, the timescale for the relaxation to a true thermal equilibrium state crucially depends on the magnitude of ∆f and is typically of Θ(1/|∆f | 2 ) [158][159][160][161]. Thus the ∆f dependence of the timescale is very important. In this section, we investigate it for linear thermalization. Our results of Sec. IV, namely Theorem and Propositions 1 and 2, take the two limits T → ∞ and ∆f → 0 in the following order [162]: [function of T and ∆f ]. (27) In this order of limits, T is smaller than any timescale that grows as ∆f → 0 [163]. That is, it extracts the behavior of the system in the time interval [0, T ] such that T = Θ(|∆f | 0 ). This leads to the following observations: Timescale in Theorem and Propositions 1 and 2: (1) If linear thermalization occurs in the limit (27), as in Theorem and Propositions 1 and 2, then it occurs in some t = O(|∆f | 0 ), as illustrated in Fig. 2(a). (2) On the other hand, if linear thermalization does not occur in the limit (27), then it does not occur in any timescale of O(|∆f | 0 ), as in Fig. 2 That is, the timescale of thermalization is independent of the magnitude of ∆f . Behaviors in longer timescales will be discussed in the following two subsections.

B. Extension to a longer timescale
We now extend our results to a longer timescale o(1/ |∆f |), and explain its implications for linear thermalization.
Let T ∆f be a timescale that grows as ∆f → 0 satisfying T ∆f = o(1/ |∆f |). As will be shown in Sec. VIII D, we find that the values of the left-hand sides of Eqs. (18), (21) and (24) That is, we obtain the following proposition [164]. Proposition 3: For any T ∆f = o(1/ |∆f |) and for any additive observableÂ, This proposition shows that Theorem, Propositions 1 and 2 hold even when T is replaced with a longer timescale T ∆f , and leads to the following observations: This result, together with the results of Secs. IV and V A, implies, in particular, that linear thermalization of the single key observableB in the timescale of O(|∆f | 0 ) guarantees linear thermalization of all additive observables not only in the same timescale but also in a longer timescale of o(1/ |∆f |).

C. Stationarity and implications for prethermalization
We here discuss stationarity of the system and its implication for prethermalization.
From Eqs. (29)-(31), we obtain the following observations: Stationarity throughout a certain time region: Suppose that conditions (20) and (23) for Propositions 1 and 2 are fulfilled. Then, throughout a time region from O(|∆f | 0 ) to o(1/ |∆f |), all additive observables take macroscopically definite and stationary values, up to O(∆f ). In other words, the system relaxes to a macroscopic state and stays in the same macroscopic state throughout this time region, up to O(∆f ).
Note that this stationary state can be either thermal or nonthermal, depending on whether condition (19) of our Theorem is satisfied or not. If the state is thermal as illustrated in Fig. 2 (19) is not satisfied] as in Fig. 2 The latter case includes systems which exhibit prethermalization [61,[152][153][154][155][156][157]. The prethermalization often occurs when ∆f switches the system from integrable to nonintegrable. (Shiraishi proved the existence of a system in which such switching is possible by an arbitrary nonvanishing value of ∆f [165].) In a typical case of such prethermalization, the system first relaxes to a nonthermal stationary state in some timescale of O(|∆f | 0 ), and then relaxes to the true thermal equilibrium state at some timescale of Θ(1/|∆f | 2 ) [61], as illustrated in Fig. 2(c). For such a system, our results detect the nonthermal stationary state in the time region from O(|∆f | 0 ) to o(1/ |∆f |) and the absence of linear thermalization in this time region.

VI. GENERALIZED SUSCEPTIBILITIES FOR CROSS RESPONSES
Our Theorem has significant meanings about the generalized susceptibilities for cross responses.
In response to the change of a parameter (such as an external field), not only its conjugate observableB but also other observables often change their values. The magnetoelectric effect and the piezoelectric effect are well-known examples. Such responses of observablesÂ that are not conjugate to the changed parameter are called cross responses, and have been attracting much attention [166][167][168][169][170]. They are characterized by the gen-eralized susceptibilities (crossed susceptibilities), which we denote by χ(A|B).
For example, when a magnetic field of an arbitrary wavenumber k, Eq. (6), is applied to a spin system, Eq. (7) yields the additive observable conjugate to f aŝ which is the total magnetizationM k of wavenumber k.
Then, if we takeÂ to be the total electric polarizationP k of wavenumber k the generalized susceptibility χ(A|B) is the magnetoelectric susceptibility at wavenumber k [171], whereas if we takeÂ =B =M k the corresponding susceptibility χ(B|B) is just the ordinary magnetic susceptibility at wavenumber k.
We here compare two types of generalized susceptibilities. One is χ(A|B) obtained in quantum mechanics, which is defined by the Schrödinger dynamics. The other is χ(A|B) predicted by thermodynamics [172], where the final stateρ eq ∆f is determined by thermodynamics and equilibrium statistical mechanics using the energy conservation, Eq. (15).
From these definitions, Eqs. (18) and (19) are equivalent to respectively. Therefore, our Theorem can be rephrased as follows.
Consistency of cross responses: Equation (35) holds for every additive observableÂ if it holds forB. Furthermore, the following symmetries follow from Eqs. (A2) and (A5) in Appendix A 2: (The latter is just a Maxwell relation of thermodynamics.) From these relations, Eq. (35) can be rewritten as Therefore, we obtain the following result. That is, the quantum-mechanical response of an additive observableB to the parameter f A conjugate to an arbitrary additive observableÂ is consistent with thermodynamics if the quantum-mechanical response ofB to its own conjugate parameter f is consistent with thermodynamics.
Moreover, this corollary has the following implication for linear thermalization. Corollary 2: Under the conditions (20) and (23) for Propositions 1 and 2, linear thermalization ofB against the quench of its conjugate parameter f implies linear thermalization of the same observableB against the quench of any other parameter f A that is conjugate to an arbitrary additive observableÂ.
These corollaries dramatically reduce the costs of experiments and theoretical calculations of linear thermalization and cross responses. For example, suppose that one wants to examine the cross response of the magne-tizationM of a spin system against the quench of an interaction parameter J, whose quench is, however, technically difficult. In such a case, one can perform an alternative experiment in which an external magnetic field that is conjugate toM is quenched. If the quantummechanical response ofM is consistent with thermodynamical one in the latter experiment, then Corollary 1 guarantees their consistency in the former experiment.

VII. EXAMPLES
In this section, we demonstrate our Theorem using onedimensional spin systems. For a nonintegrable model, we first show linear thermalization ofB, which implies, according to our Theorem, linear thermalization of all otherÂ's. We demonstrate it for typicalÂ's. We also present integrable models in which linear thermalization does not occur neither forB nor for typicalÂ's.

A. Models
If a system had many symmetries, one would have to investigate thermalization separately in individual symmetry sectors not to overlook the degeneracy of energy eigenvalues. To avoid such complicated procedures, we construct our model Hamiltonian by adding two extra terms to the one-dimensional XYZ model in a magnetic field [165] aŝ whereσ µ N +1 =σ µ 1 (µ = x, y, z). The last two terms are the extra terms that break all known symmetries, except for the translation symmetry, of the conventional XYZ model in a magnetic field. Indeed, the J yz term breaks the lattice inversion symmetry and the spin π rotation symmetry around z axis. Furthermore, the J yz and the h x terms break the complex conjugate symmetry. We believe this model, when all parameters are nonzero and J xx ̸ = J yy , is nonintegrable in the sense that it has no local conserved quantities other than the Hamiltonian, because so is the conventional model with J yz = h x = 0 [165]. As a support of this belief, we have confirmed in Appendix B 3 that the energy level statistics in each momentum sector is described by the Gaussian unitary ensemble (GUE) of the random matrix theory [88,173,174].
We tabulate the values of the parameters used in the numerical calculations as Model I in Table I. The values are taken in such a way that the ratio of every two of them is irrational in order to avoid a possible accidental symmetry.
For comparison, we also study Model II in which J xx = h z = 0 (see Table I). Although this model is slightly different from known models [175], we find it integrable in the sense that it can be mapped to a noninteracting fermionic system by the Jordan-Wigner transformation. We give its analytic solutions in Appendix C.
These models cover three typical types of systems: the nonintegrable systems, the "noninteracting integrable systems," and the "interacting integrable systems" that are solvable by the Bethe ansatz.
In all these models, we choose J zz as the quench parameter f in order to demonstrate that our additive observables are not restricted to one-body observables. In fact, the additive observable conjugate to f = J zz is the two-spin operator, We writeĤ of Eq. (40) asĤ(f ), andĤ 0 =Ĥ(f 0 ), where the initial value f 0 of f = J zz is given in Table I. Taking the initial state as the canonical Gibbs stateρ eq 0 ofĤ 0 with the inverse temperature β 0 = 0.15, we study the quench process in which f = J zz is changed suddenly.
We calculate χ QM N (A|B) and χ TD N (A|B) for several choices ofÂ including the case ofÂ =B. For this pur- We here present the results for Models I and II, whereas those for Model III are presented in Appendix B 2.

B. Susceptibilities ofB
First we calculate the two susceptibilities ofB, χ QM N (B|B) and χ TD N (B|B), and examine whether condition (36), which is equivalent to condition (19), is satisfied.
The two susceptibilities of Model I are plotted against the system size N in Fig. 3(a). They approach the same value as N is increased. The inset of Fig. 3(a) shows the N dependence of their difference, χ TD N (B|B) − χ QM N (B|B), in a log-log plot, indicating a power-law decay. From these results we conclude that Model I satisfies condition (36) and, equivalently, condition (19).
For comparison, we plot the two susceptibilities of Model II against N in Fig. 3(b). The solid lines depict their thermodynamic limits, which are calculated from the analytic solutions (C20) and (C21) given in Appendix C. These results clearly show that Model II violates condition (36) and hence condition (19).

C. Susceptibilities ofÂ
Next we calculate the susceptibilities of additive ob-servablesÂ that are not conjugate to the quench parameter f = J zz in order to demonstrate our Theorem (in the form rephrased in Sec. VI). That is, we demonstrate that Eq. (35) is satisfied for suchÂ's in Model I while it is violated in Model II.
As typicalÂ's, we choose the following observables for the demonstration. The first one is the sum of single-site observables,M The second one is the sum of two-spin observables, The third one is also the sum of two-spin observables, but the two spins are the next nearest to each other, whereσ z N +2 =σ z 2 . The susceptibilities ofÂ =M x of Model I are plotted against the system size N in Fig. 4(a), and those of A =M xx andM z1z are plotted in Figs. 5(a) and 6(a) of Appendix B, respectively. In each figure, χ QM N (A|B) and χ TD N (A|B) approach the same value with increasing N . The insets of these figures show the difference of the susceptibilities in a log-log plot. They indicate power law decays, as in the case of Fig. 3(a). From these results, we confirm that Eq. (35) is satisfied for all the three additive observables.
For comparison, the susceptibilities ofÂ =M x ,M xx , andM z1z of Model II are plotted in Figs. 4(b), 5(b), and 6(b), respectively. In all these figures, χ QM N (A|B) and χ TD N (A|B) deviate from each other. Therefore we conclude that, in Model II, Eq. (35) is not satisfied for the three additive observables [185]. This is consistent with our Theorem because this integrable model violates condition (36).

VIII. OUTLINES OF THE PROOFS
In this section, we describe outlines of proofs of Theorem and Propositions of Secs. IV and V. A combination of these outlines and the detailed discussions in Appendix A gives the complete proofs. It is known that the generalized susceptibilities can be expressed in terms of the canonical correlation [136,186,187], which is defined for arbitrary operatorsX andŶ by We start from such expressions. We then introduce an projection superoperator. We finally make use of the fact that the canonical correlation defines an inner product, i.e., it satisfies all of the axioms of an inner product (hence is called the Kubo-Mori-Bogoliubov inner product).
In Appendix A 1, we have expressed χ QM N (A|B) and χ TD N (A|B) using the canonical correlations as Eqs. (A2) and (A5), respectively. These expressions yield where δ• :=• − ⟨•⟩ eq 0 is the deviation from the initial equilibrium value, and we have introduced the Heisenberg operator,X and its long time average, Here, we have put the superscript 0 becauseX 0 (t) evolves by the initial HamiltonianĤ 0 . Since [X 0 ,Ĥ 0 ] = 0 holds for any operatorX, we can rewrite the first term of Eq. (46) as Now we introduce the following projection superoperator, It projects an operatorX onto the operator subspace whose elements (operators) commute withĤ 0 and are orthogonal to1 andĤ 0 (under the Kubo-Mori-Bogoliubov inner product). By using this superoperator, Eq. (46) can be written as Since the r.h.s. is an inner product, we apply the Cauchy-Schwarz inequality, and obtain where we have used χ QM N (A|A) ≥ 0, which follows from Eq. (A2), and χ TD N (B|B) ≥ χ QM N (B|B) [137,139], which is obvious from Eq. (51). Since we exclude phase transition points as stated in Sec. II, it is required that for anyÂ (53) from thermodynamics and equilibrium statistical mechanics [188]. Therefore, Eq. Let |ν⟩ be an eigenstate ofĤ 0 with an energy eigenvalue E ν . We take these eigenstates such that they form an orthonormal basis even if the eigenvalues are degenerate. The maximum number of resonances in condition (20) is defined by As shown in Appendix A 2, the left hand side of Eq. (21) can be rewritten as Here χ QM N (A|B; t) is the time-dependent susceptibility, which is related to χ QM N (A|B) via The r.h.s. of Eq. (55) is bounded from above by The proof of this inequality, shown in Appendix A 2, is similar to that by Short and Farrelly [146]. Combining Eqs. (55) and (58) with the condition (20), we prove Eq. (21) for every additive observableÂ. We remark that Eq. (58) can be extended to finite T , as in Ref. [146].
Here the last line follows from Eq. (A26) of Appendix A 3. By using this inequality, we have Since the time average of the left-hand side (l.h.s.) of this inequality can also be bounded from above by the same quantity, we have Here, we have used interchangeability of the time integration and the limit ∆f → 0, which is shown in Eqs. (A43)-(A45) of Appendix A 4. Combining this with the condition (23) and its consequence (25), we obtain Eq. (24). It should be remarked that Eq. (61) holds at an arbitrary time t > 0 without taking time average. That is, the variance remains small at all t > 0, although Criterion (iii) requires it only for almost all t > 0.

D. Outline of the proof of Proposition 3
We use the following inequalities that are proved in Appendix A 4, Hereχ and D 1 , · · · , D 9 are nonnegative constants of O(|∆f | 0 ) that are independent of T . By noting that the righthand sides of these equations vanish in the limit (28) Using Eqs. (33), (34) and (67), we obtain Eq. (29). Next we evaluate By taking the limit (28), we can evaluate the first and second terms from Eqs. (64) and (67), respectively. Then we have Combining this with Eq. (55), we obtain Eq. (30). Finally we evaluate By taking the limit (28), we can evaluate the first and third terms from Eqs. (65) and (67), respectively. Note that the second term vanishes in this limit because of Eq. (64). As a result, we have In the last line, we have used which follows from Eqs. (A22) and (A23) of Appendix A 3. From Eqs. (A43)-(A45) of Appendix A 4, the time integration and the limit ∆f → 0 can be interchanged, and hence Eq. (71) yields Eq. (31).

IX. DISCUSSIONS
A. Extension to nonadditive observables In the above proof, additivity ofÂ andB is not crucial. In fact, we can extend our Theorem as follows. Extension to nonadditive observables: Suppose that the observableB defined by Eq. (10) is not necessarily additive. If Eq. (18) holds forB then it holds for any observ-ableÂ that is not necessarily additive but satisfies This result will be particularly useful when discussing nonlocal properties of the system, such as the entanglement. Nevertheless, we have focused on additive observables because we are interested in consistency with thermodynamics in this paper.

B. Case of continuous change of f
We have studied the case of a quench process, in which f is jumped discontinuously. Our Theorem holds also for processes in which f is changed continuously as Here, λ(t) is a continuously differentiable function such that λ(t) = 0 for t ≤ 0 and λ(t) = 1 for t * ≤ t, where t * is a constant independent of ∆f and T . The validity of our Theorem for this process is shown as follows. χ QM N (A|B) of this process agrees with that of the quench process because χ QM N (A|B) is the zerofrequency component of a linear-response coefficient and hence independent of the time profile of f , as will be shown in Sec. IX E. χ TD N (A|B) is also independent of the time profile of f because it agrees with the adiabatic susceptibility regardless of the time profile of f . In fact, the entropy does not change in O(∆f ) because if the entropy increased by O(∆f ) then it would decrease when the sign of ∆f is inverted, in contradiction to the second law of thermodynamics. Therefore, our Theorem holds independently of details of the process.

but they are inequivalent
In this subsection, we show that the ETH forB implies condition (19) but the converse is not necessarily true.
In Ref. [135], we have shown that if the ETH (referred to as the "strong ETH" there) is satisfied for the uniform magnetization then condition (8) of Ref. [135] holds. This condition is equivalent to Eq. (4) of Ref. [135], which states that the k = 0 components of two types of magnetic susceptibilities, denoted by χ qch N (0) and χ S N (0) there, coincide. We can easily show that the proof is applicable when χ qch N (0) and χ S N (0) are replaced with the generalized susceptibilities χ QM N (B|B) and χ TD N (B|B) introduced in Sec. VI, respectively, and the uniform magnetization withB. Therefore, the ETH forB implies Eq. (36) of the present paper, which is equivalent to condition (19) as explained in Sec. VI.
Note that this result is not so obvious because, as shown in Sec. V, the timescale of linear thermalization in condition (19) is shorter than the timescale where the ETH is usually applied. In other words, when applying the ETH, one usually take T → ∞ before ∆f → 0, which differs from the limit (27) employed in condition (19).
On the other hand, condition (19) forB is weaker than its ETH. In fact, Shiraishi and Mori [24,25] constructed systems in which a certain observable violates the ETH but satisfies Criterion (i) when the initial state is an arbitrary equilibrium state at a nonzero temperature. Hence, when f is chosen as the parameter conjugate to that observable, these systems satisfy condition (19) but violate the ETH forB.

D. How to test our results experimentally
In this section, we discuss a way of testing our results experimentally. For concreteness, we discuss how to measure χ QM N (A|B) and χ TD N (A|B). Other quantities such as fluctuation can also be measured in a similar manner.
To measure χ QM N (A|B) experimentally, prepare the system in an equilibrium state [189]. This can be achieved, for example, by making the system be in thermal contact with a heat bath of inverse temperature β 0 . Obtain ⟨Â⟩ eq 0 by measuringÂ in this state. After that, detach the heat bath, so that the system undergoes the unitary time evolution. Then, ⟨Â ∆f (t)⟩ eq 0 evolves as schematically shown in Fig. 1. By measuring this evolution, one can judge whether the system relaxes to a (quasi)-steady state or not. Note that the (quasi)-steady state can be a nonthermal state because the relaxation occurs much more easily than thermalization, as pointed out at the beginning of Sec. IV. When the relaxation occurs, one obtains ⟨Â ∆f (t)⟩ eq 0 T and the relaxation time. By taking |∆f | small enough and T long enough [190], and extrapolating the data to ∆f → 0 and T → ∞, one obtains χ QM N (A|B) as Eq. (33). One also obtains the ∆f dependence of the relaxation time from the ∆f dependence of the data.
For χ TD N (A|B), its measurement might look difficult when the system does not thermalize under the unitary evolution. Fortunately, one can fully utilize a heat bath to realize equilibrium states because χ TD N (A|B) is a state function and therefore independent [apart from a finitesize-effect term of o(N 0 )] of how the equilibrium state is prepared. Firstly, set f = f 0 and β = β 0 , and obtain ⟨Â⟩ eq 0 by measuringÂ. Then, supply a small amount of energy dU to the total system (composed of the system and the bath), and measure the resultant decrease dβ of the inverse temperature and the change dB of the equilibrium value ofB. Since the thermal capacity of the bath is known, one obtains the specific heat c 0 of the system from dU and dβ. Furthermore, dB/dβ gives Then, using this and Eq. (A4), one can determine the inverse temperature β ∆f , given by Eq. (A3), of the final equilibrium state that is predicted by thermodynamics after the quench of ∆f [191]. Finally, set f = f 0 + ∆f and β = β ∆f , and obtain ⟨Â⟩ eq ∆f by measuringÂ. By substituting the measured values for ⟨Â⟩ eq 0 and ⟨Â⟩ eq ∆f of Eq. (34), one obtains χ TD N (A|B).

E. Relation to linear response theory
The quantum-mechanical susceptibility, χ QM N (A|B), is closely related to that given by the linear response theory. Hence, the consistency of χ QM N (A|B) with χ TD N (A|B), discussed in Sec. VI, is also related to the validity of the linear response theory. We finally discuss these points.
Following the pioneering studies by Callen, Welton, Greene, Takahashi, and Nakano [192][193][194][195][196], Kubo established the linear response theory for quantum systems [136]. He assumed that the system, which is initially in an equilibrium state, is isolated from environments while an external force is applied adiabatically, which means, in our notation, that Here ε is a small positive number and ω is the frequency. Then, the Kubo formula for the response ofÂ/N reads whereÂ 0 (t) is defined by Eq. (47). Note that the validity of this formula is never obvious because it depends on the degree of complexity of the dynamics. The conditions for the validity have been discussed by Kubo himself [136] and by many authors [135,137,139,186]. For technical reasons, these discussions have been made for finite N , although it is sometimes stressed that the thermodynamic limit should be taken before taking other limits such as ε → +0 [197][198][199][200][201].
If we keep N finite following these studies, we can show that [135] Therefore, the validity of the Kubo formula can be judged from the validity of χ QM N (A|B), which was discussed in Sec. VI. That is, the validity for the response of the key observable against ∆f guarantees the validity for those of other observables and for the responses against any other parameters. By contrast, the previous conditions for the validity [135-137, 139, 186]

required investigations of individual observables and responses.
It is worth mentioning that the above validity means the agreement of the response between quantum mechanics and thermodynamics. On the other hand, the classical limit of Eq. (77) yields the "fluctuation-dissipation theorem" (FDT). It states that [response at ω] = β × [correlation spectrum at ω], (79) where the r.h.s. is not a formal one but the observed spectrum. In classical systems the FDT, despite its name, holds even for non-dissipative responses [195]. Many authors discussed its quantum corrections [136,186,192,196,[202][203][204][205][206][207][208][209][210]. Although their results partially disagree with each other, it can be interpreted as due to different assumptions on the ways of measurements. However, it is recently shown rigorously that the observed fluctuation in macroscopic systems is independent of details of the ways of measurements when the measurement is performed as ideally as possible [211,212]. This universal result clarified when the FDT [in the form of Eq. (79)] holds and when it is violated [211][212][213][214].
These studies on the FDT also demonstrate that the nonequilibrium statistical mechanics is never trivial even in the linear response regime.

X. SUMMARY
We have studied how quantum mechanics is consistent with thermodynamics with respect to infinitesimal transitions between equilibrium states. We suppose that an isolated quantum many-body system is prepared in an equilibrium state, and then a parameter f of the Hamiltonian is changed by a small amount ∆f , which induces the unitary time evolution. By inspecting the expectation values and the variances of all additive observables, we have investigated whether linear thermalization occurs, i.e., whether the system relaxes to the equilibrium state that is fully consistent with thermodynamics up to the linear order in ∆f .
By comparing the long time average of the expectation value of an arbitrary additive observableÂ with its equilibrium value predicted by thermodynamics, we have obtained Theorem described in Sec. IV A. It states roughly that the two values coincide for everyÂ if and only if they coincide for a single additive observable. This key observable is identified as the additive observableB that is conjugate to f . We have also pointed out that this condition forB is weaker than the ETH forB.
To reinforce Theorem, we have then proved two propositions. Proposition 1, described in Sec. IV B, shows that the time fluctuation of the expectation value (after changing f ) is macroscopically negligible for everyÂ as long as the number of resonating pairs of energy eigenvalues (before changing f ) is not exponentially large. Proposition 2, described in Sec. IV C, shows that the variances of allÂ's remain macroscopically negligible if fluctuations ofÂ's in the initial equilibrium state have reasonable magnitudes as Eq. (23). Reinforced with these propositions, Theorem ensures that, under the reasonable conditions, the linear thermalization of the key observableB implies linear thermalization of all additive observables, which means full consistency with thermodynamics.
We have also proved Proposition 3 in Sec. V, which extends the above theorem and propositions to a longer timescale. It show that when linear thermalization occurs, it occurs in a timescale that is independent of the magnitude of ∆f , t = O(|∆f | 0 ), and it lasts at least for a period not shorter than o(1/ |∆f |). On the other hand, when linear thermalization does not occur by t = O(|∆f | 0 ), it does not occur at least in a period of o(1/ |∆f |).
Furthermore, we have shown that Theorem has significant meanings about the generalized susceptibilities for cross responses, which have been attracting much attention in condensed matter physics. Theorem guarantees that the quantum-mechanical susceptibility of every additive observableÂ coincides with the thermodynamical one if those of the key observableB coincide. We have also obtained two corollaries, described in Sec. VI, about response ofB to another parameter f A that is conjugate to an arbitrary additive observableÂ. Corollary 1 states that the quantum-mechanical susceptibility ofB to f A and the thermodynamical one coincide if those ofB to its own conjugate parameter f coincide. This is rephrased in terms of linear thermalization as Corollary 2: Linear thermalization ofB against ∆f implies linear thermalization ofB against any other ∆f A .
We have demonstrated Theorem by numerically calculating the generalized susceptibilities in three models. In a nonintegrable model, we have first shown linear thermalization ofB, and then confirmed that of otherÂ's predicted by Theorem. We have also studied two integrable models; one can be mapped to noninteracting fermions, the other is solvable by the Bethe ansatz. We have found that linear thermalization does not occur either forB or for otherÂ's in these models. This is again consistent with Theorem.
Our results will dramatically reduce the costs of experiments and theoretical calculations of linear thermalization and cross responses because testing them for a single key observable against the change of its conjugate parameter in a timescale of O(|∆f | 0 ) gives much infor-mation about those for all additive observables, about those against the changes of any other parameters, and about a longer timescale. From the linear response theory [136,137,139,186], the time-dependent quantum-mechanical susceptibility defined by Eq. (56) can be written as which results in For χ TD N (A|B), we obtain β ∆f from energy conservation, Eq. (15), as Here, since (β 0 , f 0 ) is not at a phase transition point, the specific heat We derive Eqs. (55), (58) and (22).
First, we derive Eq. (55). The temporal fluctuation of ⟨Â ∆f (t)⟩ eq 0 can be divided into two terms, We evaluate these two terms in the limit ∆f → 0.
The first term will be evaluated in Appendix A 4 as Eq. (A44). From Eq. (A43), which will also be given in Appendix A 4, the second term is evaluated as Therefore, Eq. (A6) is evaluated, in the limit ∆f → 0, as (A8) By taking the limit T → ∞, we obtain Eq. (55). Next, we derive Eq. (58). Using Eqs. (45), (A1), (A2) and the eigenstates ofĤ 0 , we have where and we have used Hence we have where, for E ν1 ̸ = E ν2 , we have introduced Using this expression, the time fluctuation of χ QM N (A|B; t) is evaluated as The term |v (ν1,ν2) | 2 can be bounded from above by using which follows from the inequality, sinh x/x ≤ cosh x.
Combining this inequality with |(ρ ν2 + ρ ν1 )/2| 2 ≤ (ρ 2 ν2 + ρ 2 ν1 )/2, we have From these expressions, we obtain Eq. (58). Finally, we show that the purity ofρ eq 0 is exponentially small with respect to N , i.e., Eq. Equations (A20) and (A21) yield Eq. (22). Note that Eq. (58) can be extended to the case where T is finite, by evaluating the term e i(Eν 3 −Eν 4 −Eν 1 +Eν 2 )t T in Eq. (A15) more carefully, as in Short and Farrelly [146]. However, such extension is meaningful only when 1/T is sufficiently small compared to the mean level spacing, which means that T has to be exponentially large with respect to N [215]. Hence we think that the difference between such extension and the original result (58) is physically unimportant, and for simplicity, we omit the precise description of such a extension. We also remark that a result similar to Eq. (58) was derived in Ref. [216].

Relation used in proof of Proposition 2
First we derive Eq. (59). Its left hand side consists of four terms, For the second and the fourth terms, we have which corresponds to the Leibniz rule of calculus. On the other hand, because the first and the third terms of Eq. (A22) give the time-dependent susceptibility ofÂ 2 and Eq. (A1) is also applicable to this susceptibility, we have Combining these with Eq. (A1), we have which results in Eq. (59). Next, inserting Eq. (A17) of Appendix A 2 into Eq. (45), we can show the following inequality between the canonical correlation and the symmetrized correlation: for an arbitrary operatorX. This was proved by Bogoliubov [217,218].

Relations used in proof of Proposition 3
In this Appendix, we show Eqs. (63)-(65) of Sec. VIII D. We also show that the time integration and the limit ∆f → 0 are interchangeable.
We introduce a function This function satisfies where χ QM N (A|B; t) is defined by Eq. (56). According to Taylor's theorem, for each t there is a constant θ t ∈ [0, 1] such that (Here the differentiability of ϕ(f, t) with respect to f follows from the concrete expression (A32) below.) Let us derive an upper bound of the absolute value of the right hand side of Eq. (A30). From the identity we have Since the integrands are bounded from above by using the operator norms ofÂ and the derivatives ofĤ(f ), we have Here, we have introduced two constants (A35) which are of O(|∆f | 0 ) since they decrease monotonically as |∆f | → 0. Using Eqs. (A30) and (A34), we obtain Eq. (63).
In addition, using Eq. (A30) we have Furthermore, using Eqs. (60), (A26) and (A1), we find that |χ QM N (A|B; t)| is bounded from above by a quantity independent of t, Combining this with Eqs. (A34) and (A37), we obtain Eq. (64). Next we introduce a function This function satisfies whereχ N (A 2 ; t) is defined by Eq. (66). According to Taylor's theorem, for each t there is a constant θ t ∈ [0, 1] such that In a way similar to the derivation of Eq. (A34), we have Combining Eqs. (A41) and (A42), we obtain Eq. (65). Finally, we notice that, in the limit ∆f → 0 Eqs. (63)- which show that the time integration and the limit ∆f → 0 can be interchanged.
Appendix B: Additional numerical results

Susceptibilities of Model I and II
In this Appendix, we perform additional calculations of the susceptibilities in Model I and II. As explained in Sec. VII A, we choose the quench parameter as f = J zz , and henceB is given by Eq. (41). We take the initial state as the canonical Gibbs stateρ eq 0 given by Eq. The results are plotted in Figs. 5 and 6. See Sec. VII C for discussions on these results.

Susceptibilities of Model III
In this Appendix, we calculate the susceptibilities of Model III (see Table I) in the same way as described in Sec. VII A or Appendix B 1.
First we investigate whether Model III satisfies condition (36). In Fig. 7, χ QM N (B|B) and χ TD N (B|B) of Model III are plotted against the system size N . This shows Model III violates condition (36). Hence, our Theorem (in the form rephrased in Sec. VI) indicates that some of the additive observablesÂ do not satisfy Eq. (35). To investigate this point, we calculate, as in Sec. VII C, χ QM (A|B) and χ TD (A|B) for three additive observablesÂ =M x ,M xx ,M z1z .
In the main of Fig. 8, χ QM N (A|B) and χ TD N (A|B) are plotted against the system size N forÂ =M xx . It is seen that these susceptibilities deviate from each other, indicating the violation of Eq. (35). The inset of Fig. 8 is the same plot forÂ =M x . In this case, Eq. (35) is satisfied in the trivial form χ QM N (A|B) = χ TD N (A|B) = 0, which follows from the spin rotation symmetry around z axis. This does not contradict our Theorem, which does not state violation of Eq. (35) for allÂ. Figure 9 shows the same plot as in Fig. 8 forÂ =M z1z . Interestingly signs, which clearly indicates the violation of Eq. (35).

Energy level spacings
In this Appendix, we study level spacing statistics of H 0 for Models I, II and III, which are defined in Table I.
Let D k be the dimension of the eigenspace of the lattice translation with an eigenvalue e −ik where k = 2πn k /N (n k = 0, 1, ..., N − 1) is wavenumber. We here label eigenvalues ofĤ 0 in this subspace using an integer j = 0, ..., D k −1 as E k j , and sort them in ascending order, E k j ≤ E k j+1 for all j. We consider the ratio of two consecutive energy level spacings E k j+1 − E k j and E k j+2 − E k j+1 in the subspace which satisfies 0 ≤ r k j ≤ 1. (When E k j+2 = E k j+1 = E k j , we define r k j := 0.) In order to investigate the level statistics in the bulk of energy spectrum, we use the ratios r k j satisfying D k /4 ≤ j < 3D k /4. By using all these ratios, we define the histogram of the ratio, P k (r). [We take the width of the interval of r for the construction of P k (r) as 0.05.] Atas et al. [174] proposed the following functional forms of P k (r) using the random matrix theory. When the energy levels obey the Poisson law, P k (r) is given by which is expected for typical integrable systems. When the energy levels obey the GUE, its P k (r) is well approximated by This is expected for typical nonintegrable systems that have no symmetries other than the lattice translation. Furthermore, when the levels obey the Gaussian orthogonal ensemble (GOE), its P k (r) is well approximated by P GOE (r) := 2 27 8 (r + r 2 ) (1 + r + r 2 ) 5/2 which is expected for typical nonintegrable systems. We have numerically calculated P k (r) of our Models by the exact diagonalization, and compared the results with the above forms. Note that the symmetry properties of the subspace of wavenumber k sometimes depend on whether N is even or odd. In addition, the properties at special wavenumbers k = 0, π are sometimes very different from those at other wavenumbers. Therefore we consider three values of wavenumber k = 0, 2π/N, π and both even and odd N . (Results for k = π with odd N are absent, since k = π is impossible when N is odd.) We plot P k (r) of Model I for N = 21, 22 in Fig. 10. It is well fitted by the GUE form P GUE (r) (solid line). For comparison, the Poisson form P Poi (r) and the GOE form P GOE (r) are depicted by the dashed lines, which clearly deviate from P k (r). These results indicate that Model I has no local conserved quantities other than the Hamiltonian.
In Fig. 11, P k (r) of Model II is plotted for N = 21, 22. It is well fitted by the Poisson form P Poi (r) (solid line), while it clearly deviates from the GUE form P GUE (r) (dotted line) and the GOE form P GOE (r) (dashed line). These are consistent with the integrability of Model II.  In Fig. 12, P k (r) of Model III is plotted for N = 20, 21. The P k (r) for k = 2π/N is well fitted by the Poisson form P Poi (r) (solid line). This is consistent with the integrability of Model III. On the other hand, P k (r) for k = 0, π show δ-function-like behaviors, which are much different from the GUE form P GUE (r) (dotted line) and the GOE form P GOE (r) (dashed line). The δ-function-like behaviors come from the fact that almost all energy eigenvalues have certain degeneracies [219]. These degeneracies indicate that the subspaces of k = 0, π have additional non-Abelian symmetries, which would be related to the integrability of Model III [220].
The operatorR takes 1 for the states with even number of fermions, whereas it takes −1 for the states with odd number of fermions. Hence depending on whether the number of fermions is even or odd, terms containingR in Eq. (C2) become the antiperiodic or periodic boundary condition, respectively. To resolve this boundary condition that depends onR, we consider two types of the Fourier transformation. We introduce two sets of wavenumbers K e := 2πn k N + π N n k = 0, 1, ..., N − 1 , K o := 2πn k N n k = 0, 1, ..., N − 1 , and perform the following Fourier transformatioñ where e iϕ++ := J ++ /|J ++ | and θ k ∈ (−π/2, π/2) is determined from