LiHoF4 as a spin-half non-standard quantum Ising system

We would first like to thank the referees for their very valuable comments. These, we believe, helped us improve the quality and presentation of our work. Here, we will address each comment and point to corresponding amendments made in the revised manuscript.

Referee 1:

1. Carefully check the formulas, correcting apparent errors (e.g. the presence of a two-body interaction in H_{3B}) and ambiguities of notation (e.g. the lack of a site index on the crystal field terms and the apparent confusion between H_{full} and H_{micro}).

We added a site index to the crystal field term in Eq. (1) and clarified the caption of Fig. 2, explaining precisely what $H_{micro}$ means. It is $H_{full}$, as is defined in this work, with the addition of a hyperfine term. Other curves also include the hyperfine interaction, but since, in those cases, it is omitted from the Hamiltonian and later accounted for by a temperature-dependent renormalization, we deemed it appropriate to use a different label for when it is included from the start. Regarding the presence of a two-body interaction in $H_{3B}$, we maintain that it is indeed a three-body interaction due to its dependence on the existence of a third spin $i$, as evident in the presence of the index $i$ in, e.g., $V_{ij}^{xz} V_{ik}^{xz}$. In our view, this is an important distinction, as it pertains to the x-dependent behavior of the diluted variant $LiHo_{x}Y_{1-x}F_4$; namely, it explains the relatively mild reduction of T_{c} with decreasing x, as explored in Ref. [30].

2. Expand on the content and justification of the phenomenological scaling and renormalisation procedures, indicating clearly how the results would be different if they were not included (and ideally presenting results without the scalings / renormalisations implemented, so that the reader can see their effects).

The two scaling and renormalization procedures referred to are indeed not optimal. However, we would like to note that both adjustments are not related to the main result of the paper, namely the reduction of $T_c$ at low fields. For this reason, we purposefully chose procedures that were already used in previous papers (particularly Refs. [28,37]), so that comparison is straightforward. To make the manuscript more self-contained, and in view of the referee's comment, we added Appendices D and E, which discuss the temperature-dependent renormalization that accounts for hyperfine interactions and the rescaling of the longitudinal interaction, respectively. Each appendix includes a new figure (Figs. 6 and 7) that presents the current results with and without these adjustments.

Furthermore, to more directly relate to the main result of the paper, we have included Appendix F, which discusses the analytical behavior of $T_c(B_x)$ at low fields.

3. Ensure that sufficient information, including all parameters (e.g. J_{ex}), is provided for the reader to reproduce the results presented. Give references for the numbers quoted, e.g. the values of alpha and rho given below (4).

The omission of the numerical value of $J_{ex}$ is an oversight on our part, which has been amended (below Eq. (2)). We thank the referee for spotting this issue. The other parameters mentioned, $\alpha$, $\rho$, and $\Delta$, are derived in this work by numerical diagonalization of $V_{c}$ and their values are given below Eq. (4). We use the same notations as in Ref. [31], though the values differ slightly because of the use of updated crystal-field parameters. We add a footnote pointing out this difference and its origin.

4. Discuss the validity of the truncation of the Schrieffer-Wolff transformation. (This arises especially because of the observation on page 4 that the three-body terms are long-range because of the additional summations. Doesn't that raise the issue that higher orders in the Schrieffer-Wolff transformation might yield four-, five-, etc.-body interactions that are equally important?)

We acknowledge the referee's point regarding higher-order terms in the Schrieffer-Wolff transformation and have given it considerable thought. Indeed, higher-order terms in the Schrieffer-Wolff transformation should yield interactions involving a growing number of spins, e.g., the third order in the Schrieffer-Wolff transformation will yield 4-body terms, as well as possibly additional two- and three-body terms. The latter are easily discounted, as their relative contribution compared to their respective counterparts from the second Schrieffer-Wolff order is diminished by a factor of $\langle V \rangle / \Delta \sim 1/10$, where $\langle V \rangle$ denotes, symbolically, the energy scale associated with the dipolar interactions.

However, as the referee rightly points out, four-body interactions cannot be naively discounted based solely on the relative smallness of their coefficient. To emphasize this point, we mention a reference that has recently come to our attention (Rau et al., Phys. Rev. B 93, 184408 (2016)), wherein the authors perform a similar perturbation procedure on the pyrochlore XY antiferromagnet $Er_{2}Ti_{2}O_{7}$, and make the case that due to the combinatoric factor associated with four- and six-spin terms, their contribution is comparable to that of lower-order terms in the perturbation series. The combinatoric factor is the number of ways an open chain of nearest-neighbor spins can be chosen in the pyrochlore system (in that work, only nearest-neighbor interactions are considered), which gives an increase of two orders of magnitude because of the high coordination number of the pyrochlore system (compared to $LiHoF_4$). In the present case, instead of a combinatoric factor, it is a sum over presumably long-range four-,five-,etc.-body interactions that seemingly may result in a significant enlargement of the otherwise inherently small higher-order terms.

Despite all of the above, we assert that higher-order terms are negligible in this case. First, we make the case that for an order-of-magnitude estimate, it is enough to consider only nearest-neighbors, despite the long-range nature of the interaction. Consider, for example, the three spins depicted in Fig. 1 of the manuscript, assuming spin 1 is a nearest-neighbor of spin 2 and spin 3 is a (different) nearest-neighbor of spin 2. The energy associated with such a triplet from the three-body term is $ \frac{\alpha^2 \rho^2}{\Delta} E_D^2 V_{12}^{xz} V_{23}^{xz} \approx 0.01 K $. Considering that we have 4*3=12 ways to choose such a chain of neighboring spins in the $LiHoF_4$ crystal, we can estimate the energetic contribution to be $ \approx 0.1 K $, which turns out to be a good estimate of the actual MF correction (the last term in Eq. (7)), which is ~0.1 K. In essence, nearest neighbors dominate the total sum despite its long-range nature because the dipolar-derived emergent interactions provide both positive and negative contributions that partially offset. This would be the case for four-spin (and higher) interactions as well.

Let us apply the same reasoning to a hypothetical four-spin interaction that emerges from the next order in the Schrieffer-Wolff transformation. The energetic contribution of such a term would be $ 4\times 3 \times 3 \times \alpha^2 \rho^4 E_D^3 [V_{12}^{xz}]^3 / \Delta^2 \approx 5 mK $, so its effect on $T_{c}$ would negligible in the context of the present work.

We add the above explanation as a closing paragraph to Appendix A, where H_{eff} is derived.

Even worse, data with two different values of the longitudinal-interaction scaling factor are presented in the paper.

As an additional point, under "Weaknesses," the referee mentions that two longitudinal-interaction scaling factors are used. This is correct in that Fig. 2 shows results without ODD terms with two rescaling factors (0.785 and 0.805). However, we would like to clarify that the former is just a reproduction of results presented in Ref. [37] (specifically, the solid curve labeled "x=1" in Figure 7 therein). To clarify this point, we slightly amended the phrasing of the caption of Fig. 2 to emphasize that the solid line as a whole is taken from Ref. [37] and not just the rescaling factor.

Referee 2:

1) I would have liked a better summary justifying various rescaling factors invoked (with refs. to experimental works and "strong c-axis fluctuations" -- i do not recall the issue and not inclined to go digging in the literature), esp. considering the paper is all about fixing relatively small discrepancies

1) better justfiication of rescaling factors

As per our response to the first referee's request #2, we have added Appendices D and E that aim to justify the different rescaling and renormalization schemes used and present results without them, showcasing their effect.

2) it might have helped to state the symmetry constraints on the 3-body terms, e.g. I assume they are required to be even in \sigma_z operators so as not to break two-fold symmetry. These would limit the first of these terms, which apparently only contains 2 Pauli operators (both \sigma_z^s), to be a two-spin operator and thus be considered properly as E_D/\Delta order reduction of the bare 2-spin dipolar term rather than a 3-spin term. if this is the case, the authors should not call it a 3-spin term AND also address whether this term by itself (without the true 3-spin terms) can reproduced the reduction of the critical temperature or not.

2) symmetry constraints on and possibly correction/reclassification of 3-spin terms

The comment regarding symmetry constraints on emergent terms is, of course, correct. We partially address this fact in Appendix A, right before Eqs. (10)-(13). In the revised manuscript, we rephrase that statement to make that point more explicitly and add a similar statement right after Eq. (6).

As for the reclassification of the first term in $H_{3B}$, as we mention in our response to comment #1 by the first referee, we reaffirm that its classification as a three-body interaction is justified. The interaction fundamentally involves three spins, with the third spin $i$ mediating the emergent interaction—an observation that is crucial to explaining the behavior of the system under mild dilution, as considered in Ref. [30].

As to whether this term alone can reproduce the reduction in the critical temperature—that is precisely what we find in this work. Namely, that this term dominates the reduction, at least in MF at zero-field, where the other two terms in $H_{3B}$, which we denote "quantum terms," completely vanish due to lattice symmetry. To emphasize this point, we slightly change the wording following Eq. (7).

Not only is the term in question dominant at zero-field, but we also show in a new appendix (Appendix F) that its behavior under non-zero transverse field is inline with experimental observations.

Outside of MF, we believe these terms may contribute to the reduction, albeit in a sub-dominant manner, possibly accounting for the remaining discrepancy compared to the experimental $T_{c}$. We raise this possibility in the Conclution section.

3) it seems the 3-spin terms in the Hamiltonian imply the existence of non-classical expectation values induced by classical 2-point correlators, i.e. finite values of \<X> and \<Z> (possibly with shape-dependent spatial profiles?). If true, the authors may want to discuss these as possible predictions for future experiments.

Off-diagonal dipolar terms indeed lead to a non-trivial distribution of internal magnetic fields, which subsequently induce finite local expectation values \<X> and \<Y>. This phenomenon was partly addressed in Ref. [30] by examining the local x-direction fields in a Monte Carlo simulation (see Fig. 4 in that reference). While a similar analysis is possible within the analytical low-energy effective Hamiltonian approach undertaken in the current paper, it does not follow trivially from the three-spin terms in $H_{3B}$ and is thus outside the scope of the present work.

- Added missing summation sum in Eqs. (14)-(15).
- New appendix: D The hyperfine interaction
- includes new Figure 5
- New appendix: E Rescaling the longitudinal interaction
- includes new Figure 6
- New appendix: F Extended analytical characterization of $T_{c}$ in non-zero $B_x$
- includes new Figure 7
- Previous appendices "D Simulation details" and "E Effective low-energy description of the Fe8 molecular magnet" moved to G and H, respectively
- Minor correction: in the paragraph preceding Figure 3, removed redundant k index.
- Added an index to the first term in Eq. (1)
- Amended the paragraph following Eq. (2) to include the value of $J_{ex}$
- Added a footnote on page 4 remarking on the values of the parameters $\rho$, $\alpha$, and $\Delta$
- Added a footnote on page 4 explaining the presence of the first term in $H_{3B}$
- Added a sentence mentioning the symmetry requirements on $H_{eff}$ in the paragraph following Eq. (6)
- Added a clarification at the end of the same paragraph (the one following Eq. (6)) regarding the naming distinction between the first term and last two terms in $H_{3B}$
- Changed "the result of the three-body interaction" to "the result of the first of the three three-body terms in $H_{3B}$" after Eq.(7)
- Added brief explanation and reference to Appendix F at the end of Section 4
- Exlicitly refer to the absence of terms with an odd number of Pauli z operators right before Eq. (10)
- slightly amended the phrasing of the caption of Fig. 2 to emphasize that the solid line as a whole is taken from Ref. \[37\]
- Refined the definition in Eq. (25) and modified the surrounding text (Appendix C) accordingly, emphasizing that $J^{\mu}_{\text{eff}}$ is defined within the low-energy subspace. The referring text (first paragraph of Section 4) has been slightly changed to accordingly
- Two paragraphs added to the end of Appendix A justifying the cutoff of $H_{eff}$ at second-order in $H_{T}$