Extracting average properties of disordered spin chains with translationally invariant tensor networks

Comment: Summary: This work presents an original approach to compute disorder-average properties of random-spin chains using translation-invariant tensor networks, by exploiting the restoration of translation invariance upon disorder averaging. Although the idea of using ancilla qudits to represent the disorder and perform the average was already present in the literature, this was so far mainly used for dynamics. The present work highlights how to get finite temperature disorder average results, focusing on two ways of obtaining the correlation lengths. The manuscript is clear and concise and tests the approach on a non-trivial test case, where key signatures of the presence of the infinite-randomness fixed point are recovered. Except for the minor comments below, it seems to me that the work easily meets the expectations and criteria for this journal.

[Comments/questions:] I have minor comments and questions. The first two, though minor, should be addressed. The others are written only out of curiosity and interest in the work, or are marginal suggestions for presentation improvements. They may be disregarded if the Authors do not find them helpful or relevant.

Response: We thank the Referee for the positive assessment of our work and for the recommendation. We also thank the Referee for the constructive questions and suggestions.

Comment: 1- In the text, (\Lambda = \mathrm{Tr}_{\sigma} (N(\tau)e^{-(\tau + \Delta \tau)H})). However, Fig. 1 and its caption suggest that (\Lambda) is obtained by tracing out the spin degrees of freedom in (\tilde{\rho}(\tau)). It seems to me that this distinction is very small, but important. Could the Authors check this and either correct the figure accordingly or clarify why this distinction does not matter?

Response: In Fig. 1, the MPO tensors in (a) and (b) indeed have different meanings. We have modified the figure to make this more clear. We have also corrected the caption of Fig. 1(b) accordingly.

Comment: 2- Context/outlook: It seems to me that the specific context could be slightly better spelled out for the reader, with a small consequence for the outlook. (a) The authors cite the work of Paredes et al, but it would be nice to highlight how the present work is different in spirit. Similarly, I think the work https://scipost.org/SciPostPhys.6.3.031 may be given as part of the context. (b) In the outlook, the authors ask the question of whether a variational version of their algorithm could directly target ground-state properties. With proper context, it is clear that the question is about the properties of their algorithm rather than about targeting ground states, which was already discussed in the previous references. (c) Slightly less directly connected, it may be interesting to comment on the existence of numerical SDRG (see e.g. https://arxiv.org/pdf/2501.02643 for a recent example on a related model) for ground states, and the fact that they can access notably dynamical exponents and typical behavior. This would be an opportunity to highlight better what the present method brings. I also think it may make sense to either cite the review by Igloi & Monthus or the finite-temperature work of Young (PRB 1997) on the transverse-field Ising chain.

Response: We thank the Referee for the helpful comments and suggestions. (a) The works by Paredes et al. and SciPost Phys. 6, 031 (2019) (now added as a citation in the manucsript) are indeed closely related to ours. The main difference is that they primarily focus on real-time dynamics, for which the normalization procedure introduced in our work is not required. We have updated the manuscript to make this distinction clearer. (b) The simulation of ground states has indeed been discussed in Paredes et al. However, their method relies on adiabatic evolution, which is different from the imaginary-time evolution we use. Nevertheless, as the Referee suggested, it is indeed useful to highlight this difference and to clarify our further research direction accordingly. We now do this explicitly in the conclusions section of the new version of the manuscript. (c) In the conclusions section we have added an additional sentence highlighting the possibility of combining our method with the numerical SDRG approach as a potential direction for future research. We have also added references to the work of Igloi and Monthus, and the finite-temperature work of Young, which we indeed should have already cited in the first version of the manuscript as these works are highly relevant. We thank the referee for bringing them to our attention.

Comment: 3- (Optional) Entanglement spectrum: It seems natural to perform the Schmidt decomposition on (\rho) rather than (\tilde{\rho}). However, this density matrix is a prior of a different type than that which occurs in translation-invariant problems. It is therefore tempting to ask if : (a) the entanglement spectrum decays fast enough to justify truncation; (b) there is any interpretation of the validity of the approach (e.g., validity of "typicality") / an expectation of where it should fail (e.g., Bose glass in another model? ). Otherwise stated, could it be that the approach works better close to the IRFP than deep in the localized phases ?

Response: (a) In the updated manuscript we have added a figure of the operator entanglement spectrum of \rho (appendix B). For all temperatures used in our simulations we see a clear exponential decay of the operator entanglement spectrum, which justifies the bond dimension truncation. (b) For the random Ising chain we have also studied the case with non-zero $\delta$ to generate the data in Fig. 7. In our simulations, we encountered no obstruction to going deep in the localized phases. Based on our intuition for tensor networks coming from clean systems, we expect that the more local the physics becomes, the easier it should be to capture it with a tensor network. But it could very well be that disorder actually helps in capturing the critical behaviour, i.e. that infinite randomness physics is easier to describe with tensor networks then criticality in clean systems.

Comment: 4 - (Optional) Possible marginal improvements to the presentation (in order of appearance in the text):

a) Although it is well-known in the community, a specific reference for the statement "The community has reached a point at which low-temperature equilibrium properties of most local spin and fermion Hamiltonians on a one-dimensional lattice can be studied with rather limited computational resources, including many gapless systems" may be helpful to newcomers.

Response: We thank the Referee for the suggestion. We have added references to this statement in the updated manuscript.

Comment: b) (\tilde{\rho}) is introduced in Eq. (3), but its connection to (\rho) is only explicitly clarified at the bottom of page 2. Clarifying this connection from the start may help the reader.

Response: We agree that this ordering of the presentation was not optimal and we have updated the manuscript accordingly. In the updated manuscript, the connection between (\tilde{\rho}) and (\rho) is introduced before we discuss the details of TEBD algorithm.

Comment: c) Fig. 3, right panel: The figure is slightly difficult to read/interpret (in particular because of course, it is a good collapse). I have two suggestions : (i) could the authors also indicate the inverse temperature? (ii) Could the lowest temperature be shown with crosses or other symbols that do not hide the larger temperatures?

Response: To address the comments of the other referee, we have generated new data (to study e.g. the dependence of our data on the number of disorder values) and we have made new plots that now replace the data collapse plot in the current version of the manuscript. We have done this because the new data allows us to make contact with additional theory predictions.

Requested changes: 1- Check if a modification is needed with respect to (\Lambda) in Figure 1 (see point 1 in the report).

2- Better clarify what the existing literature has already established (see point 2 in the report).

Response: As explained in more detailed above, we have updated the manuscript to address both these comments.

Comment: The authors develop a tensor-network approach to simulate quantum systems with disorder avoiding the need to run simulations multiple times to sample over various disorder configurations. The authors achieve this by introducing an auxiliary disorder qudits. They benchmark the algorithm on the paradigmatic model of disordered quantum systems – random transverse field Ising model. I think the problem is timely and this manuscript is an important first step toward the solution. However, there are several reservations that prevent me to recommend the manuscript in its current form to be published in SciPost Physics:

Response: We thank the Referee for recognizing the importance and timeliness of our work. However, we believe that the Referee may not have fully appreciated the main contribution of our work. Besides, among the five main reservations that prevent the Referee from recommending the paper, we note that only the first two are actual concerns about the manuscript; the remaining comments are rather questions/suggestions about the paper.

Comment: 1. First of all, this is a methodology paper – it does not report new physics, but the algorithm that potentially can solve problems in future more efficiently. This is fine, but for the methodology paper, I would expect more thorough analysis of complexity and computational costs of the new method compare to the existing ones. I would naively say that the authors save a lot of computing run time by the price of higher memory costs and sacrificing the access to some information (typical scaling, probability distribution etc). This has to be properly discussed, ideally the gain and loss compared to the standard techniques have to be quantified. I would also recommend to move the discussion on the performance of the algorithm from supplementary materials to the main text.

Response: First, we thank the Referee for pointing out the importance of an analysis of the computational complexity. The complexity of our method follows directly from the standard TEBD approach. The only additional cost arises from the auxiliary disorder qubits and from the additional normalization procedure after each time evolution. In response to the Referee’s concern, we have revised the manuscript to present the scaling of the computational cost (both runtime and memory) more clearly in the main text. We also added discussions on possible directions for further reducing the computational cost arising from these steps.

Regarding the difficulty of accessing information such as typical scaling or probability distributions, we have added a new section to the paper (section 4 in the new manuscript), where we explain how to obtain the distribution of correlation lengths from our MPO, and how to extract the typical correlation length from it. On the methodological side, this is now possible because we modified the algorithm: we now truncate the MPO without tracing out the disorder qudits, and in this way the truncation tries to optimally represent the density matrix for each disorder configuration, and not only the averaged density matrix.

Finally, we emphasize that our approach is a new method with promising potential, but it is still far from its optimal form. Many possible improvements remain to be explored in future studies.

Comment: 2. My second concern is the usage of the random Ising model as a benchmark. The authors made a very big assumption when going from the continuous distribution of random coupling/field to a discrete with a very few values. This assumption has to be justified. Ising model is very dangerous choice for benchmarking – everything works on Ising model; it is exactly solvable. I am not at all convinced that discrete and small number of random parameters would work in general. And this is the main bottleneck of the method --- it would be extremely expensive to go to a larger values of the discrete parameter or to make them different on every site (simple with open boundary conditions, but seems non-trivial for transitionally invariant infinite tensor network).

Response: Regarding the concern about changing the random coupling/field distribution from continuous to discrete values, we note that in Fisher’s original work (Phys. Rev. B 51, 6411 (1995), Sec. IV.B.1) it was shown that a very simple toy model where $h<1$ is constant and $J=1$ with probability $p$ and $J=0$ with probability $1-p$ is already able to reproduce all the thermodynamic properties of the weakly disordered phase of the model with continuous disorder (i.e. the rare-region-dominated quantum Griffiths phase on the disordered side). This shows that discrete disorder is able to reproduce the characteristic physics of the random transverse field Ising chain.

To show that it is possible to increase the number of disorder values in our simulations, we have added new data obtained with 4, 9 and 16 disorder values. We now also mention that this only leads to a linear increase in memory cost. We notice no significant increase in runtime when using more disorder values (so memory is the main bottleneck). In particular, we find that increasing the number of disorder values does not increase the bond dimension required for the MPO representation of the density matrix. For this reason, we do not agree with the Referee’s statement that this becomes "extremely expensive". We also show that the dependence of our data on the number of disorder values agrees with theory. Furthermore, we have increased the values of $\ln\beta$ in the new data in the manuscript (the maximal correlation lengths we get are now $\sim 50$) to further strengthen the claim that our method can indeed achieve competitive results with moderate resources (corresponding to bond dimensions of order 100).

Regarding the concern about using the Ising model as a benchmark, we respectfully note that the model studied here is not the typical Ising model in the sense of the statement that "everything works on the Ising model". In the context of "everything works on the Ising model", one usually refers to the clean 1D quantum Ising model (or, equivalently, the clean 2D classical Ising model). At its critical point, the quantum Ising model is described by a $c = 1/2$ CFT, and due to its small central charge, it is relatively easy for many numerical methods (colloquially speaking, the Ising model is the least critical of all CFTs). In contrast, the disordered Ising model considered in this work is fundamentally different. The disorder acts as a relevant perturbation and intrinsically changes the physics of the system. Its long-distance behavior is governed by a different RG fixed point, namely the infinite-randomness fixed point. Besides DMRG and (free fermion) ED we are not aware of many other numerical methods that have been applied to the random transverse field Ising model. So it is not clear at all that "everything works for the random transverse field Ising model".

It is of course very well possible that only infinite randomness physics can easily be captured by our method, because e.g. it leads to small operator entanglement for the density matrix (this is a possibility -- we don't know if it's true), or because it can be captured with discrete disorder. Whether these properties hold more generally for other disordered models is an interesting open question for future research. This being said, we are not aware of any algorithm that has tried to directly obtain average quantities via a single-shot simulation. So we still believe that the fact that our method works so well for this model is a highly non-trivial result. And even if most disordered models exhibit different universal physics with discrete disorder, the discrete case still presents a general class of well-defined disordered models and it remains an interesting question to understand their physical properties.

Although we focus on translationally invariant infinite tensor networks, the spirit of our work can be readily extended to finite systems with open boundary conditions. In that setting, it is also straightforward to allow the disorder distribution to be different on each site. On the other hand, we believe that such an extension goes against the main motivation of our work, which is to develop an algorithm that explicitly exploits statistical translational invariance.

Comment: 3. The authors analyzed only the correlation function and comment that extraction of the entanglement entropy is complicated by construction of the algorithm. But what about the energy gap and/or the gap between the Schmidt values?

Response: The distribution of energy gaps should only make sense for finite systems, as these are finite-size gaps. We work directly in the thermodynamic limit, so we are not really sure what energy gaps we should try to extract.

In the appendix of the new manuscript (Fig 4) we plot the Schmidt values corresponding to bipartition of the MPO in a left and right infinite half. However, we note that these Schmidt values determine the operator entanglement entropy, and not the physical entropy. The physical entropy of a subregion is obtained by tracing out the complement of that region, and then diagonalizing the resulting reduced density matrix $\rho_R$ to obtain $-\text{tr}(\rho_R\ln\rho_R)$. This is in general very costly, even for an MPO representation of the density matrix (this is related to the fact that MPOs can have volume law scaling for the physical entropy). But for the efficiency of our algorithm, only the operator entanglement matters, and we now check that this remains small and an MPO with managable bond dimension can be used.

Comment: 4. The authors claim that they lose access to the typical values and the overall statistical distribution. However, I wonder whether this information is in fact contained in the final network that the authors optimize. I also wonder if there is a way to post-process the final optimized MPS to extract the distribution, and how computationally costly this could be.

Response: We thank the referee for this question and for stimulating us to further think about typical correlations. As mentioned above, we have added a new section to the manuscript where we explain how to obtain typical correlations and the distribution of correlation lengths from our MPO. Obtaining the distribution of correlation lengths from the MPO is extremely cheap: it does not involve any further optimization or bond dimension trunction, one just has to multiply $D\times D$ matrices.

Comment: 5. The authors use known critical exponents of the infinite-randomness critical point to produce a data collapse. However, what would be the strategy and performance of the method if these exponents are not known or non-universal (as in the case of the random XXZ chain)?

Response: For systems where the exponents are unknown, the data collapse procedure can indeed be more challenging. However, this challenge arises in any numerical approach, and our method does not introduce additional limitations in this respect. Also, in the new version of the manuscript we have replaced the data collapse with plots containing new data, such as e.g. the average spin-spin correlation for different disorder distributions. We have made this choice because this allows us to show that our results agree with additional theory predictions.

Comment: In addition, I have a few minor comments that the authors might want to consider: 1. In the introduction, the authors write that: “randomness can cause entanglement to be inhomogeneously distributed… which significantly increases the number of sweeps needed for convergence.” In my experience, DMRG convergence for disordered systems is in fact much faster than for the clean one. Could the authors provide references reporting slow convergence, or did the authors have in mind the works of Ref. 7–9?

Response: The statement is based on works in the literature, such as e.g. [11] (Brenden Roberts and Olexei I. Motrunich, "Infinite randomness with continuously varying critical exponents in the random XYZ spin chain", Phys. Rev. B 104, 214208 (2021), section III. A). In that work, the authors report that DMRG is susceptible to spurious convergence to excited states in systems with strong disorder. We also explicitly mention that this is most likely to occur in random-singlet-type models -- we do not claim that this happens in all disordered models. But this statement is made in the introduction, and is therefore only meant to provide a broad context.

Comment: 2. At the end of the first column on page 1, the authors refer to a “relatively small bond dimension.” It would be useful to provide at least the order of magnitude here.

Response: We adjusted the manuscript and now provide a numerical value.

Comment: 3. The authors claim that they use a modified version of TEBD to ensure the correct normalization. More details here (on the normalization conditions and how TEBD helps with it) would be useful.

Response: The modified TEBD is "modified" in the sense that we introduce an additional normalization step after each time-evolution step to ensure that the density matrix remains properly normalized within each disorder sector. The relevant normalization condition is given in Eq. (3) and further discussed in the paragraph following that equation. For the normalization procedure, we employ a variational procedure (the VOMPS algorithm detailed in Appendix A). The TEBD algorithm itself does not help with the normalization procedure, and it is the subsequent variational step that enforces the correct normalization.

Comment: Figure 3 appears on top of page 3, but it is first discussed only on page 4.

Response: We have adjusted the location of the figure accordingly.

Comment: 10.a The authors say they use a “uniformly distributed” random field h_n “between [0.73, 1.3].” Is this uniformity in the linear or in the logarithmic scale?

Response: We have adjusted the main text to make it more clear that it is uniform on the linear scale.

Comment: 10.b Also, how do the authors “fine-tune” the system to the critical point? The value of delta in Eq.6 seems very large: in numerical simulations it is typical to keep this value below 10^{-3}-10^{-4}.

Response: We have performed new simulations with new disorder distributions. The data in the new manuscript is obtained exactly at $\delta = 0$. In the random Ising model, in order to have $\delta=0$, one only needs to make sure that the disorder configuration of couplings and fields are the same. In the appendix we show the correlation length divided by $(\ln\beta)^2$ as a function of $\delta$, for different temperatures. We see a clear crossing of the different temperature curves exactly at $\delta =0$. This shows that our numerics indeed finds the correct location of the critical point.

Comment: 1. The comment on “important contributions from rare regions” requires a reference.

We added a citation to Fisher's 1992 paper, which mentions the importance of rare regions already in the introduction.