Multiple relaxation times in perturbed XXZ chain

Understanding the properties of nearly integrable systems seems to be crucial for establishing the relation between theoretical studies of integrable models and experimentally realizable systems. The essential consequence arising from the integrability-breaking perturbations concerns the relaxation of correlation functions which occurs at a finite time scale. Then, there is an obvious question of whether a single perturbation introduces a single relaxation time or multiple relaxation times. Up to our best knowledge, we provide the first direct numerical evidence for the latter scenario. We also formulate a conjecture which explains why typical observables show exponential decays with a single relaxation time. We believe that this message is important for the ongoing research on integrable and nearly integrable systems.

Motivated by the Referee’s criticism as well as by the first remark of the Referee 1 we have significantly changed the abstract and conclusions. Refer to our response to Report 1 for more details.

**The referee writes:**

>It would be nice if the authors could clarify what is difficult or not when trying to identify the relevant changes needed to explain a given result. For example, an algorithm for determining such relevant charges is mentioned but never used.

In the revised manuscript (below Eq. 6) we describe the basic idea of the numerical algorithm and explain why it is applicable only to systems that are much smaller than the systems studied in the present work.

>Page 2, the last paragraph is not very clear. Another study of the non integrable XXZ chain focusing on anisotropies Delta=0.5 and Delta=1 is mentioned, however, all numerical data in the previous sections is for Delta=0.5.

We have newly written this paragraph. We stress that the case Delta=0.5 has been used throughout the manuscript. However, in the last section concerning the decay of operators which do not commute with S^z_{tot} we additionally study the case with Delta=1. In contrast to other commensurate values of Delta, for Delta=1 one may easily find a conserved quantity that does not commute with S^z_{tot}. This simplicity was the motivation for discussing the Delta=1 case in the last section of our manuscript.

>What is the origin of the oscillations which start appearing at very late time in several plots? Are those numerical artifacts?

The time dependence of the correlation functions consists of a smooth part (that is relevant for the thermodynamic limit) and tiny finite-size fluctuations which are unavoidable in an arbitrary finite system. The smooth part decays in time and eventually becomes of the same order of magnitude as the finite-size fluctuations. Since we plot the numerical data on the log-scale, the fluctuations seem to be amplified for a large time.

>Page 11, in the first paragraph, it is mentioned that the observables were carefully selected to have small overlap onto the modes with slow decay. How was this done exactly ? By trying several parameters values and observing the late time results, or from first principles?

In the revised manuscript (page 11) we explain the main idea, which allowed us to tune the overlap. Namely, we made use of the result from Ref. [3] that the projection of spin current on the energy current is proportional to S^z_{tot} thus the overlap may be tuned via changing the magnetization sector.

We thank the Referee for pointing out this unclear point of our presentation. We have formulated our conjecture for the multi-scale relaxation scenario at the end of the second paragraph of conclusions that starts with “Typically, the opposite holds true.” Still, it is just a conjecture that requires additional studies, and, for this reason, it was not sufficiently stressed in the previous version. We have modified this part of our manuscript. We have also substantially changed the abstract.

**The referee writes:**
>My main criticism about the approach is the following. The authors study numerically chains with L=25 sites up to times of order 10^3. One would expect that to address relaxation one should consider times t\lesssim L, avoiding revival effects. I would appreciate if the authors could address and clarify this point.

In systems with periodic boundary conditions (which we use) the size L can directly determine characteristic time scales in correlation functions only in certain anomalous cases, e.g., in integrable systems with ballistic transport. On the other hand, for perturbed integrable models, the transport is (normal) dissipative and the relevant length scale is the effective mean-free path decreasing with increasing the perturbation strength, g. Provided that the mean-free path is smaller than the system size, correlation functions should become essentially L independent (which we monitor by comparing with results for smaller systems, which we, however, do not display in the manuscript). Therefore the accessible system sizes set the lower bounds mostly on the accessible values of g. Since we are not able to directly evaluate the mean-free path, we use a parameter range for which correlation functions do not reveal any significant L dependence but also exhibit clear quadratic dependence of g. In order to address this problem we have rewritten the paragraph on page 3 that starts with: “As a first step, we establish the range of the perturbation strength,…”
The decay of operators which don’t commute with S^z_{tot} does not show g^2 dependence because the decay originates from a different physical mechanism, i,e. from lifting the degeneracy. In this case, we explicitly discuss that the numerical results do not show any significant L-dependence.