Phase Diagram and Conformal String Excitations of Square Ice using Gauge Invariant Tensor Networks

We thank the referee for her/his thoughtful review of the manuscript.

Following the referee's request, we have modified the title by replacing "Gauge Invariant Tensor Networks" with "Gauge Invariant Matrix Product States", which is indeed more specific to the simulations that we performed using MPS confined to the gauge-invariant sector of the spin-ice model.

Based on the referee's initial remark, we have further added a new section in the appendix (Appendix A) in which we make the connection between our gauge-invariant 1D-mapped MPS and the underlying general construction of (abelian and non-abelian) gauge-invariant tensor networks in [5], based on the quantum link model (QLM) formalism. We hope that these additions make it more clear in how far the construction in this manuscript differs from others (especially [6],[19-22], where various types of gauge-invariances have been encoded in tensors with special internal structures). As the referee has pointed out, the QLM construction parameterizes the gauge-invariant sector of the Hilbert space in the local computational basis, but also introduces artificial abelian U(1) local symmetries on the links, which in turn can be encoded with the well-established framework for global abelian symmetries. We find this a very practical and canonical way to reproduce the special internal (block-)structures required for gauge-invariant tensor networks.

Following the referee's detailed remarks, we have carefully reviewed our manuscript, as reported in the following:

(1)
We have clarified our notation as suggested. In the original manuscript, the symbol $s_i$ was used to denote spin variables, the local Hilbert space, as well as individual canonical spin-1/2 basis states on some lattice link $i$. We now use subscripts $\mu$ to refer to links in general, and use the subscripts $\mu_1$, ..., $\mu_4$ to denote the local Hilbert space which the operators acts on. We also introduce $s_\mu$ to label $\sigma^z$ quantum-numbers (and respective eigenstates in bra-ket notaion) in the beginning of subsection 2.3 which are then used throughout the manuscript. In Eq. (4), the bold-faced $\boldsymbol{s}$ (vector) has been replaced by $\boldsymbol{v}$. In context of spin-correlations of Eq. (8), $s$ still appears as a spin coordinate, but only in subscripts.
In the definition of $\omega_y$ at the end of subsection 2.1, we replaced $s$ by $\sigma_z$.
Furthermore, we now denote charges with $\bar{c}$ to distinguish from the $c$ used for central charge and fit constants of the form $c_i$.

(2)
We have rephrased and added a more descriptive wording, the sentence in question reads now:

"[...] on a lattice of corner-sharing cross-linked squares (the two-dimensional analogue of a pyrochlore lattice), in the limit of strong anisotropy $J_{xy}\ll J_z$ [28] (an equivalent mapping can be obtained in Ising models, see Ref. [27])."

The concept is described in more detail in the Refs. [27][28] cited immediately after the paraphrased results.

(3)
The term originally referred to the antiferromagnetic $\sigma^z \sigma^z$ interaction between linked spins, as discussed in Ref. [28] and reviewed in Ref. [27]. We have recast the sentence such that it is now clear it refers to the references.

(4)
The continuum limit of quantum link models can be reached using dimensional reduction, if the theory allows for a Coulomb (deconfined) phase in one dimension more (see the discussion in Phys. Rev. D 60, 094502 (1999), which we have added as Ref. [42] to the reference list). This is indeed applicable to the model we consider.

There is actually one tricky point which deals with Lorentz invariance, that is explicitly broken on the lattice. In particular, the original proof in [43] strictly requires this symmetry to be present, so it is not immediately clear how the conclusions of [43] are applicable to lattice problems at all. This is partly discussed in [27].

We have modified the sentence mentioned by the referee as follows:
"The fact that confinement is the only possible scenario here is related to the continuum-limit behaviour of the theory (recovered via dimensional reduction in quantum link models, see Ref. [42]), where confinement is due to monopole contributions [43]."

(5)
We thank the referee for bringing this ambiguity to our attention. We have resorted to $\epsilon$ to denote singular values.

(6)
We did not perform measurements of the overhead due to quantum-numbers book-keeping. We agree that a comparison with other gauge-invariant tensor network methods would be highly interesting. The referee points out that we must ultimately expect an exponentially scaling number of blocks in the transverse size, which together with the other scalings discussed in Appendix D limits the tractable system size. This is of course expected for any MPS parameterization of such a two-dimensional system due to the necessary growth in bond-dimension.
Generally speaking, a large number of blocks enables our simulations by lowering the memory- and runtime cost of typical blockwise tensor-operation during the TEBD which scale polynomially with the block-dimensions but usually only linear in the number of blocks. A typical block consists of up to several hundred double-precision tensor elements (depending on state, bond-dimension and system size); each block has a single one-dimensional physical index and two bond-indices for which some dimensions can be deferred from Fig. 15.
By using optimized numerical codes for handling abelian symmetries in tensor networks, we expect that overhead added in form of a few integer calculations and indexing operations per block (e.g. additions and inversions from the quantum-number fusion rule; lookups in, and creation of, the internal memory layout holding those blocks) is small in comparison to the resources spent in blockwise tensor-operations. It could however become well noticeable in absolute terms, especially when approaching the RK-point where block-sizes become small.

(7)
We have added Fig. 3, and expanded on some of the details in Appendix C and Fig. 13.

(8)
We have added a paragraph towards the end of subsection 3.1, referring to (now) Fig. 6, including a sentence with more details on how we obtain the respective error bars for the transition points:

"The reported error-bars account for the uncertainty in determining the intersection, which is also limited by simulation errors affecting the value of $O$ (see Appendix E for estimates), which around $\lambda_c$ tend to be dominated by limitations in the bond-dimension."

Since we do not use DMRG but rather a 'Trotter-ized' time-evolution protocol with dynamically scaled time-step, we do not use the truncation errors of the individual time-steps. Instead, as reported in Appendix E, we can fit most of our results against an empirical model Eq. (26) which is sensitive to the bond-dimension and additionally the targeted "precision" (determining the duration of the imaginary time-evolution). Albeit an exploration of the error in terms of correlation length would be interesting, our method appears to be sufficient in the 2D spin-ice setting where finite-size corrections tend to be much more significant than the error due to finite bond-dimension and simulation precision.

(9)
In the paragraph added with our reply (8) in subsection 3.1, we clarify how we extract the phase-transition point from rescaled curves: "In detail, our values reported in Fig. 6 have been found from the intersection of curves [...] rescaled for the different widths $L_y$. To this end, the exponent $\gamma$ was tuned to make all three curves intersect at the closest possible values $\lambda$.".

Given the restrictions in system sizes, and the weak first-order nature of the transition which limits the scaling regime, we do not report extrapolated critical exponents. In light of this, we refer to the parameter $\lambda_c$ at which the phase transition occurs more consistently as 'transition point'.
We decided against using Binder-cumulants, which can be implemented in form of MPO measurements of the higher order correlations but quickly become impractical with the growing (index-)dimensions due to the underlying 2D systems.

(10)
We are dealing with the full entanglement entropy, which is the sum of the two contributions mentioned by the referee. We are actually planning a systematic study of the two separate contributions in the future, as those are immediately accessible with our algorithm (since the reduced density matrices we are dealing with are already in block diagonal form in what is the electric field basis in the Casini-Huerta-Myers approach to the problem).

(11)
This is a point that requires clarification. Our goal here is not to quantify excess entanglement per se, but rather, to understand how entanglement between two regions is affected by the presence of a string connecting them.

This is very hard to do in a fully rigorous manner: the main issue is that inserting charges in the system does necessarily introduce 1) some additional boundary effects, and 2) since the string has a width comparable in size with $L_y$, also some effects due to self-interaction. This is why we resorted to the simple entropy difference $S_{\text{diff}}$: this quantity is clearly defined, convenient to measure, and describes qualitatively (and, in case the points 1 and 2 above can be neglected, quantitatively) the entanglement between two regions given by the string. We interpret this as the string entanglement entropy -- that is, the entropy generated by the presence of a string on the top of the vacuum. This is reminiscent of the string energy discussion in literature.

In case the string state has the same entropy, this would imply that the string itself does not generate any additional entanglement between two subregions. It is possible that, due to extremely strong boundary effects, this actually does not correspond to a correct result. In the case mentioned by the referee, it might be that the entanglement of the string excitation is compensated by some self-interaction effects due to the finite width. We purposely stayed away from regimes where this can happen, by considering only the $\lambda>-0.2$ region. These regimes would require transverse system sizes larger than the one we could access.

(12)
Our statement concerning the impractically large allocations of 100GB and more referred to "standard TEBD", as opposed to the first-order truncated evolution exponentials used in our implementation. We have rephrased our statement and moved it to Appendix C, above Eq. (24), to clarify this point.

(13)
To the best of our knowledge, this relation had only been discussed in 3+1-d gauge theories. In the context of our work, the goal was to show that tensor network simulations offer an alternative angle to this perspective, allowing to tackle it from an entanglement perspective. An understanding of the string in terms of bosonic theory from microscopics is extremely challenging. In the 3+1-d QLM, a discussion has been presented in Ref. [62]. While we would be willing to add this piece of discussion, we prefer not to do so, as the argument presented therein do no immediately transfer down to the 2+1-d case and might generate some confusion. We actually find the 2+1-d problem an intriguing open one in the field, that deserves future work on its own.

About statistics: This shall indeed be possible to check using simulations of a full (quasi-)adiabatic exchange between two strings. Since with tensor networks one has access to the full wave function, quantities related to interference are accessible. However, this presently lies well beyond our computational capabilities.

Following the referee's comment, we have added the following sentence to the conclusions after "....systems despite achieving modest system sizes.":

"From the theoretical viewpoint, our results motivate the search for a microscopic derivation of the bosonic string theory of U(1) quantum link models that could complement the numerical results presented here, and possibly extend those to parameter regimes in the vicinity of the transition point between RVB and Neel phase."

We thank the referee for her/his report and the suggestions, which we considered in the resubmitted manuscript as follows:

(1)
In subsection 2.3 we have added Fig. 3 to give a pictorial representation of the system, the blocking and the some parts of the symmetric MPS construction. We additionally accompany the paragraph mentioned by the referee with some more general notes in Appendix A.
Following the remarks of the referee on subsection 2.4, we have moved to, and expanded on the details and equations in Appendix C, and added Fig. 13.

(2)
In principle we can use open boundary conditions in the short direction of the MPS, e.g. by explicitly fixing vertical spin-configurations at $m=1/2$ and $m=L_y+1/2$ in the local basis, which in turn also fixes $w_y$. However, the computational cost would be comparable, and we opted for the periodic boundary conditions as to moderate finite-size and boundary effects.

We have added a corresponding sentence at the beginning of Sec. 3, at the end of the first paragraph, where it complements the discussion on aspect ratios:
"While our algorithm is amenable to both periodic- and open boundary conditions along the y-direction at a comparable computational cost, we opted for the cylindrical conditions in all simulations as to minimize finite-size and boundary effects."

(3)
We have cited these works in the introduction together with the other LGT-TN works, as well in our closing remarks in the conclusion right after: "[...] more specialized and powerful and TN classes such as projected entangled pair states".