Effect of next-nearest neighbor hopping on the single-particle excitations at finite temperature

We are thankful to the referee for various constructive comments. Below, we respond to those comments pointwise while the corresponding changes in the revised manuscript is also mentioned.

Comment: A detailed comparison with the existing literature at half-filling is missing. Only stating that the findings are similar to those of cluster perturbation theory is not sufficient. It is unclear which features are in (qualitative or quantitative) agreement with CPT. Furthermore, there are many techniques which provided results that allow for comparison, some of them going beyond the limitations of the mean-field technique used here, some being at a comparable level. These differences need to be discussed. In particular, the authors might want to compare their results to numerically exact results on the PG phase of the Hubbard model, e.g. using lattice Quantum Monte Carlo (QMC), determinant Quantum Monte Carlo (detQMC) or Variational quantum Monte Carlo (VMC).

Reply: We would like to clarify that the exact-diagonalization + Monte-Carlo (ED+MC) based method used in the current work incorporates all types of thermal and spatial fluctuations in the order-parameter fields associated with the antiferromagnetic order. This is in contrast with the standard meanfield theories which ignore both thermal and spatial quantum fluctuations. Thus the magnetic moments, unlike the static mean-field theory, does not melt in the method used here when the system approaches $T_N$. Similarly, in the Hartree-Fock approximation based approaches, the spectral function can develop a gap only below $T_N$. On the contrary, we find that the gap persists beyond $T_N$, which we refer to pseudogap-like gap. Another advantage of the current approach is that as the temperature increases, thermal fluctuations start dominating over the quantum fluctuations, therefore, the accuracy of ED + MC increases in capturing essential physics.

A variety of techniques have been utilized which goes beyond standard meanfield theories, which mainly include density-matrix renormalization group (DMRG), dynamical-meanfield approximation (DMFT), cluster DMFT, QMC, detQMC, VMC, etc. The applicability of DMRG has largely been restricted to the quasi-one dimensional systems. Although, the DMFT-based methods do capture the Mott transition, the self-energy correlation incorporated is independent of momentum rendering it not suitable for the study of momentum-dependent spectral features. The QMC and detQMC were used extensively to study finite-size system with the difficulty that, at low temperatures, the correlations length is greater than the lattice size. Therefore, the correlations are overestimated for smaller clusters because they are artificially closer to criticality than a system in the thermodynamic limit. This may result in the disentanglement of the MI and AFM transitions. Another important consequence of the finite-size effect is a good momentum resolution for the single-particle spectral function, which has remained challenging till now.

Another important issue with QMC and DetQMC is the sign problem, which severly restricts its range of applicability especially when the hopping beyond nearest-neighbors is considered, which introduces particle-hole asymmetry. Even when doping away from half-filling is considered, the DetQMC can used only in a limited temperature range.
For these reasons, despite a significant volume of work at half filling, most of them were restricted to only nearest-neighbor hopping, particularly, the studies based QMC, detQMC, etc. The effect of $t^{\prime} \ne 0$ has been mainly studied using cluster-perturbation theory (CPT). CPT combines the solutions of small individual clusters of an infinite lattice system with Block theory of conventional band description to provide an approximation for the Green's function in the thermodynamic limit. In particular, at intermediate interaction strength, it may be challenging to make accurate prediction for the Hubbard gap mainly constrained by finite-size induced level splitting. We have incorporated the summary of all the points discussed above in the revised manuscript.

Comment:Also the technique itself should be discussed in view of other techniques, in particular of mean-field type. How does it differ from standard mean-field techniques, slave bosons, slave spins, Hubbard-1 approximation etc.

Reply: Slave-bosons or -spins mean-field theoretic approaches have been quite successful in studying the Mott transition. Recent developments also show their uses in studying the magnetically ordered phases. However, again inclusion of magnetic order ignores the spatial and thermal fluctuation, which is contrastingly different from the method we adopt in this work.

Comment: It is unclear why the discussion focuses on the high-symmetry path $(\pi/2,\pi/2)-(\pi,0)$. Changing $t^\prime$ results in a change of the dispersion on the non-interacting level, $U=0$, and reshapes the Fermi surface. It is thereby clear that some changes of the spectral function are simply due to the modifications of the dispersion, the shift of van-Hove singularities etc. In order to assess the effects of interactions, it would be useful to provide the non-interacting dispersion or the Fermi surface. To make statements on the 'pseudogap-like' features of the spectral function, its behavior should be studied along the Fermi surface with respect to the non-interacting reference for each $t^\prime$.

Reply: We were primarily focused on how the gap opens along the Fermi surface. While the Fig. 4 and 5 showed the dispersion with quasiparticle momentum-resolved amplitude, Fig. 6 and 7 showed how the gap evolves along the Fermi surface. In Fig. 7, there is typo, the high-symmetry direction $(\pi/2,\pi/2)-(\pi,0)$ should be replaced with Fermi surface along $(\pi/2,\pi/2)-(\pi,0)$. Along the Fermi surface, gapless excitation is present, therefore, we did not plot it. We have corrected the typo, that is, replaced the phrase ``high-symmetry direction''
by ``Fermi surface $(\pi/2,\pi/2) \rightarrow (\pi,0)$''.

Comment: Referring to the dip in the density of states (DOS) as a 'pseudogap-like feature' can be misleading. The pseudogap is a $k$-dependent phenomenon that does not necessarily need to be related to a dip in the DOS. Such a dip can in principle also stem from a (fully) gapped system at low temperature, whose gap is progressively filled by thermal excitations at higher temperature. The authors' technique includes spatial fluctuations in the spin sector, thereby they should have access to the $k$-dependent self-energy which can give them a precise (and direct) answer on the existence or absence of a pseudogap.

Reply: In using the phrase 'pseudogap-like feature' for the current context, we have followed the convention used for the hole-doped cuprates. In other words, with rising temperature, the gap should be filled up. Conventionally, it is expected that the gap should disappear as the temperatures increases up to $T_N$. As seen from Fig. 2, the gap continues to persist despite the loss of long-range order, which is also reflected in the momentum-resolved spectral function. In such a scenario, we needed a terminology to describe this phenomenon. For this reason, we have used the phrase.

Comment:Figure 1: A legend should be added to make clear which quantities are plotted in panels a) and b) - the spin structure factor and the NN spin-spin correlator.

Reply: In the revised manuscript, we have added the legends separately to the Fig. 1 (a) and (b) to avoid any confusion.

Comment: Figures 4 \& 5 seem to show the very same data, plotted in a color map and a line plot respectively. Figure 5 should therefore be removed or replaced by a plot which shows the amplitude of the spectral function along the high symmetry path for the different temperatures.

Reply: We have removed Fig. 5 and included a Figure showing the amplitude of the spectral function along the high symmetry path at different temperatures.

Comment:The presentation of Figures 6 \& 7 is slightly confusing and needs to be improved, in particular depicting the $k$-path more clearly. Also, it is not clear why the path $(\pi/2,\pi/2)-(\pi,0)$ can be compared in this way as a function of $t^\prime$, which even modifies the dispersion of the non-interacting model.

Reply: The confusion here arises because of a typo. Actually, the ${\bf k}$ path in both the figures are along the Fermi surface along the path $(\pi/2,\pi/2)-(\pi,0)$ instead along the high-symmetry direction along $(\pi/2,\pi/2)-(\pi,0)$. We have modified the caption to remove this confusion.

Comment: The whole study is done at $U=4t$, which is less than half the bandwidth of the system. Nevertheless, the authors seem to suggest that they can safely interpret the broadening of the spectral function in therms of the Heisenberg spin exchange $J = 4t^2\/U$. This is problematic for different reasons: i) $U=4t$ is far from the Heisenberg limit; ii) the NNN spin exchange $J^\prime$ is not taken into account, despite the presence of the NNN hopping $t^\prime$; iii) the broadening of the spectral function should be rather described by (the imaginary part of) the self-energy, which the authors do not show.

Reply: It is true that $U=4t$ is far from the strong coupling limit where the Hubbard model can be mapped to the Heisenberg model. Here, we would like to clarify that though the results are presented in the manuscript only for $U = 4t$, we did check the calculations for higher value of $U$ also and found that the thermal broadening increases with increasing $U$. For this reason, we wanted to comment as to what is expected when $U$ is increased. When incorporated NNN spin exchange $J^\prime$ because of NNN hopping $t^\prime$, the broadening will be enhanced further because of the frustration introduced.

Comment: Claims on agreement with experiment need to be substantiated. Which cuprates, which measurement techniques and which studies do the authors have in mind? What does it mean to be "in good agreement with experiment"?

Reply: In the current work, we are primarily interested in the effect of next-nearest neighbor hopping on the single-particle excitation of half-filled Hubbard model. Therefore, the phrase "in good agreement with experiment" was used in reference to the ARPES measurements carried out in the undoped cuprates, particularly, the momentum-dependent single-particle gap sructures (PRL 74, 964 (1995), PRL 80, 4245 (1998), PRB 70, 092503
(2004)).


Comment: In the discussion and conclusion, the authors should be careful in assessing the transferability of their results to doped cuprates: 1) Their technique cannot be applied in a straight-forward way to doped systems since it neglects several types of fluctuations, which are known to be important in these systems. 2) The pseudogap of doped systems is not necessarily the same pseudogap that the authors study here at half-filling with a mean-field approach that is tailored to capture the physics in the strong-coupling limit.

Reply: We agree that our results are not transferable to the doped cuprates because then several types of order-parameter fields will come into picture and the $d$-wave superconducting order parameter is prominent amongst them. It is indeed true that the pseudogap feature that we discuss may be entirely different from the one arising as a result of multiple competing orders in the doped cuprates. Precisely, for this reason, we have used the phrase ``pseudogap-like'' instead of psuedogap at various points in the manuscript.

Comment: Given the vast literature on the $t-t^\prime$ Hubbard model at half-filling, the authors should explain more carefully what they mean when saying that their study fills a 'long-standing gap'.

Reply: As discussed in the reply to earlier comments, to the best of our understanding, most of the earlier studies, which go beyond mean-field theory, have focused on the half-filled Hubbard model without next-nearest neighbor hopping mainly because of particle-hole asymmetry induced sign problem. Furthermore, CDMFT or CPT may provide a picture corresponding only to small clusters, thus suffering from finite-size induced level splitting. On the other hand, the method that we have used for the simulation, beside being free from sign problem, has the advantage of accessing a large-sized system, thus able to obtain a momentum-resolution never obtained before. A very good momentum resolution, which is free from finite-size effect, is absolutely necessary to conclusively establish the existence of small gap as found in the pseudogap or pseudogap-like phases along the Fermi surface. It is in this respect, we used the phrase 'long-standing gap'. We have modified the phrase in the revised manuscript to avoid any emphasis on the ``long-standing gap''.

General comment: In the present manuscript the authors provide an analysis of the effects induced by a nearest-neighbor hopping amplitude on the half filled Hubbard model. In particular, they focus on the spectral properties along the node and antinode in the pseudogap regime. The present studies extends a previous work of the same authors (Ref. 34) to a finite $t^\prime$. They conclude that large values of $t^\prime$ are unfavorable to the pseudogap behavior.

Reply: We are thankful to the referee for various constructive comments. In the following, we reply to specific comments by listing them pointwise.

Comment: The method they used combines several approaches, and involves a series of limitations. Specifically, the method can not be extended to finite doping because the charge degrees of freedom are frozen (both spatial and thermal fluctuations are neglected).

Reply: (i) We would like to clarify that the method can be extended to finite doping, however, the effective model should be modified accordingly. An appropriate effective model will be $t-J$ model, which is a derived from the Hubbard model. The $d$-wave superconductivity can arise from the next-nearest neighbor attractive interactions present in the $t-J$ model and magnetic order will arise from the nearest neighbor exchange interaction. Alternatively, one can use $t-U-V$ model as well to investigate the consequence of competing $d$wave superconductivity, where $U$ and $V$ are on-site repulsive and next-nearest neighbor attractive interaction parameters.

(ii) In the current work, we were primarily interested in investigating the role of next-nearest neighbor hopping on the single-particle spectrum at half filling, because, to the best of our knowledge, this problem has remained largely unexplored due to the challenges arising from (a) small cluster size limiting the momentum resolution and (b) famous sign problem in the implementation of simulations based on quantum Monte Carlo (QMC) or determinant QMC. Notably, these methods go beyond the mean-field approximation, and can incorporate the thermal and spatial fluctuations.

(iii) Furthermore, at half filling, the effect of the thermal fluctuations in the order-parameter field, which corresponds to the magnetic moments in the current work, is expected to be significantly large on the single-particle spectral function in comparison to the charge fluctuations. This is because a small rotation of magnetic moment does not require as much energy as involved in the charge fluctuations corresponding to the double occupancy of a site, unless the system is doped with the charge carriers.

Comment: Similarly, the treatment of $d$-wave superconducting fluctuations is not included. In the conclusions, the authors state that additional auxiliary fields can be introduced, details are however missing.

Reply: As we mentioned earlier that we were primarily focused on the role of next-nearest neighbor hopping on the single-particle spectrum. It would be of strong interest to know how the presence of two competing order such as magnetic and $d$-wave superconductivity affect the spectral features, particularly, the pseudogap-like features. In the discussion, we have a added a brief discussion as to how the additional auxiliary fields can be introduced via $t-J$ model, which is a derived from the Hubbard model for small dopings. The $d$-wave superconductivity can arise from the next-nearest neighbor attractive interactions present in the $t-J$ model. Alternatively, one can use $t-U-V$ model as well to investigate the consequence of competing $d$wave superconductivity, where $U$ and $V$ are on-site repulsive and next-nearest neighbor attractive interaction parameters.

Comment: Given the restricted applicability of the approach, I would have appreciated a paragraph or short section about a concrete route that allows to overcome the above restrictions. I also think the limitations of the method should be mentioned in the methods section and not only in the discussion at the end.

Reply: We have incorporated suggested discussion as well as mentioned the limitations of the method in the method section.

Comment: The findings on the $t^\prime$ dependence are interesting and worth to be published, but given the limitations of the approach I would rather recommend SciPost Physics Core as suitable journal, after the points raised have been addressed.

Reply: Our earlier work was mostly focused on the development and testing of the new method which combined three techniques namely, traveling-cluster approximation, parallelization of the update process and twisted-boundary conditions. The sole aim of combining these three techniques is the accessibility to a large-system size critical for a better momentum resolution for the spectral function so that any conclusion on the nature of spectral gap may be free from any finite size effect. Then, we demonstrated this idea with the help of half-filled Hubbard model with only nearest-neighbor hopping.

However, the actual correlated-electron system described by an effective one-band model may involve long-range hopping. The system of high-$T_c$ cuprates is one such standard example. Presence of the long-range hopping can have significant impact not only on the phase diagram but also on the spectral properties etc. which deserves a separate treatment. Several studies in the past have attempted to investigate the role of $t^{\prime}$. However, the approaches mainly include Hartree-Fock approximation based mean-field theory, dynamical mean-field theory (DMFT) etc., which do not incorporate the spatial fluctuations in the order-parameter field. The former does not even take into account the thermal fluctuations in the order-parameter field.

On the other hand, the methods such as QMC, detQMC etc., which go beyond mean-field theories, also use finite cluster size and suffer from famous sign problem in the absence of the particle-hole symmetry, \textit{i. e.} when the long-range hopping is considered. Even if doping is introduced, the use of these methods may be restricted to a certain temperature range.

In our simulation, in the absence of particle-hole symmetry, the chemical potential may fluctuate with temperature, therefore, the simulation requires to check the chemical potential and update the same in accordance with half filling at intermediate steps. This enhances the computational cost in the simulation process. Another factor which is responsible for raising the computational cost is the frustration introduced by the next-nearest neighbor hopping in the system, which raises the equiliberation time.

To the best of our knowledge, the current work may be the first study, which has examined the temperature dependent role of next-nearest neighbor hopping in the half-filled Hubbard model on the spectral features with a momentum resolution which is almost free from finite-size effect. For all the above reasons, we believe that the manuscript may deserve the visibility associated with SciPost Physics.

We are thankful to the referee for the encouraging comments on the current work. We also appreciate specific comments, that have been helpful in improving the manuscript further.

Comment: The introduction is currently organized via a back and forth from theory and experiments to discuss the two main aspects that this work wants to address, i.e. the asymmetry of the phase diagram and the temperature and momentum-dependence of the quasiparticle excitations. I think it would be helpful to separate the two so that the reader will have first a clear idea of the experimental aspects that will be considered followed by an overview of the theoretical state of the art in this context. I invite the authors to consider such a reorganization of the introduction.

Reply: We have now restructured the introduction so that the parts giving the theoretical and experimental backgrounds are now separate, begining with the experimental aspects first, and then followed by an overview of the theoretical status.

Comment: The description of Fig. 7 is misleading. In the text the authors claim that the raising of energy moving from $(\pi/2,\pi/2)$ to $(\pi,0)$ is almost linear, however the plots show a more complex behavior (those curve are not even close straight lines). The only temperature at which the energy curves are lines is at $T~ T_N$ were the grows is also quite small. I think that what describe the plot saying that the increasing is monotonic would generate less confusion.

Reply: We have modified the corresponding sentence in the manuscript so that there is no confusion. Now, instead, it is mentioned that the hole-excitation energy increases monotonically upon moving from $(\pi/2,\pi/2)$ toward $(\pi,0)$ and it varies almost linearly in the vicinity of $T_N$.

Comment: Again in Fig 7. Can the author find a x-label easier to understand looking at the plot without reading the caption? One option would be putting explicitly $(\pi/2,\pi/2)$ and $(\pi,0)$ as stating and ending point.} \\

Reply: We have now used symbols, which are frequently employed to represent the high-symmetry points in the Brillouin zone. In particular, it may be noted that $(\pi/2,\pi/2)$ is also one of the high symmetry points in the reduced or magnetic Brillouin zone. Thus, $(\pi/2,\pi/2)$ and $(\pi,0)$ are represented by X$^{\prime}$ and M$^{\prime}$, respectively. We hope that with this change, the $x$-label will now be easier to understand.