Dynamics of string-like states of a hole in a quantum antiferromagnet: a diagrammatic Monte Carlo simulation

The dynamics of a hole motion in a quantum antiferromagnet has been studied in the past three decades because of its relationship to models related to superconductivity in cuprates. The same problem has received significant attention because of its connection to very recent experiments of the dynamics of ultra-cold atoms in optical lattices where models of strongly correlated electrons can be simulated. In this paper we apply the diagrammatic Monte Carlo method to calculate the single-hole Green’s function in the t–J model, where the J term is linearized, in a wide range of imaginary-time with the aim to examine the polaron formation and in particular the details of the contribution of the so-called string excitations found in such recent experiments. We calculate the single-hole spectral function by analytic continuation from imaginary to real time and study the various aspects that constitute the string picture, such as, the energy–momentum dependence of the main quasiparticle peak and its residue, the internal excitations of the string which appear as multiple peaks in the spectral function as well as their momentum dependence. We find that the earlier analysis of the spectral function based on a mobile-hole connected with a string of overturn spins and the contribution of the internal string excitations as obtained from the non-crossing approximation is accurate.


Introduction
It has been three decades ago when the high-T c superconductivity in cuprates was discovered and shortly afterward the t-J model emerged [1] as a simple non-trivial model to describe certain basic features of the environment experienced by the correlated electrons in these materials. The Hamiltonian of the t-J model contains the Heisenberg antiferromagnetic term [2] and a hole-hopping part where the double occupancy is prohibited due to strong on-site Coulomb repulsion. While one could question the degree of its applicability to such a complex problem as the superconductivity in the cuprates, when one attempts to understand the interplay of the charge and spin degrees of freedom in any strongly correlated system, this model is as basic as the Ising model for understanding phase transitions. Yet, even in the simpler regime of the t-J model, regarding the single-hole motion in the presence of quantum spin fluctuations in an otherwise ordered but fluctuating antiferromagnetic background, an exact solution is absent and, numerically, there are still very interesting unresolved issues.
The motion of a hole inside a quantum antiferromagnet as described by the t-J and related models has been extensively studied using a number of analytical, semi-analytical and numerical techniques [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. One of the main and well-known conclusions of these studies is that the hole becomes a well-defined low-energy quasiparticle with the minimum of the quasiparticle band at (±π/2, ±π/2), which is an unusual position in the Brillouin zone for a minimum to occur, and, thus, at sufficiently light doping, a hole pocket develops near (±π/2, ±π/2). These predictions were confirmed later by means of angle-resolved photoelectron studies (ARPES) of the copper-oxide based undoped insulator [19]. These ARPES studies find good agreement with the theoretical predictions [13] for the bandwidth of approximately 2.2J, where J is the antiferromagnetic coupling by using a value [2] for J 0.13 eV and the empirically known value of t 0.4 eV. By adding relatively small direct next-nearest-neighbor hopping terms (t and t ) in the t-J model [20][21][22][23][24][25][26][27], one can also obtain agreement with ARPES along the direction (π, 0) → (0, π).
Recently, novel theoretical [28][29][30][31][32][33][34][35][36] and experimental activity in atomic physics and specifically in many-body systems of ultra-cold atoms in optical lattices [37][38][39][40] has rekindled the interest in the problem of the hole motion in a quantum antiferromagnet and the formation of the magnetic polaron [41] as well as the nature of hole-spin correlations in strongly correlated systems. Some of these studies use a cold-atom quantum simulator to directly observe the formation dynamics and subsequent spreading of individual magnetic polarons. It becomes possible to directly measure the density and spin-resolved evolution of a single hole in a 2D Hubbard insulator with short-range antiferromagnetic correlations. They reveal fast initial delocalization and a dressing of the spin background, indicating polaron formation. More specifically in an experimental realization of the Hubbard model in systems of ultra-cold atoms in an optical lattice [39], it was found consistency with geometric strings, entities that may explain the relationship between hole motion and spin order. This discovery is in contrast to the one-dimensional case where in a time-resolved study of fermionic Hubbard chains using ultra-cold atoms one finds spin-charge deconfinement [40], as expected.
Such string features associated with the charge motion in an antiferromagnetically correlated spin-background, which we will refer to in the present paper as 'string excitations', which are internal excitations of the hole polaron, are well-known [3-6, 8, 9, 12, 13, 18]. In particular, in references [12,13], where using the non-crossing approximation (NCA) [7] (the NCA is the same as the self-consistent Born-approximation) a detailed calculation of the hole spectral function was reported, it was found that these excitations appear as several peaks in the hole spectral function which, from their scaling behavior as a function of J/t, were identified to be the internal excitations of the hole bound by the string to its 'birth' place by a linear 'string' potential. This string picture as revealed by the t-J model was used later [42] to explain the waterfall-like features of some later findings of ARPES [43]. It has been argued [42] that the string carried by the hole as it moves inside an antiferromagnet can be the source of pairing. In addition, if verified by a more rigorous approach, the emergence of the string potential and the excitations of the hole in this linear potential and possibly its pairing to another hole due to the string induced interaction, is an interesting phenomenon in its own right with an interesting analogy in quark confinement in quantum chromodynamics [44].
However, the presence of the string states and its internal excitations has not been clearly and unequivocally revealed by means of other reliable methods. More importantly, one should be skeptical about a picture which emerges from the NCA, which ignores vertex corrections that are especially important in low-dimensions and in the case of low-spin. The reasons other computational methods have not clearly justified the presence of the internal excitations of the string are the following: (a) first, the exact diagonalization technique can only be done on small-size systems. The spectral function of such finite-size systems is a sum of discrete δ-function peaks (broadened by an artificial finite resolution). However, these string-excitation peaks are the result of adding a continuum of such δ-function peaks in the spectral function of an infinite system. (b) Second, quantum Monte Carlo can only be done in imaginary-time and analytic continuation to real-time lacks the resolution required to reveal such a detailed structure and its scaling with t/J and it only yields a broad high-energy peak [18].
In this paper we will closely examine the validity of the details of the string picture, as it arises within the t-J model when the antiferromagnetic interaction is linearized, with a rigorous computational approach, the diagrammatic Monte Carlo (diag-MC) [45][46][47][48][49] as been implemented in reference [50] in conjunction with the flat-histogram method. An advantage of the diag-MC for our problem is that it automatically yields the results for the infinite lattice and, thus, we do not have to worry about finite-size effects. The diag-MC method has been very successfully applied to similar problems and, in fact, to the single-hole Green's function of the t-J model and related models [47][48][49], albeit with a slightly different focus. Namely, the analytic continuation to real-time did not provide the resolution needed to distinguish the multiple peaks which are the internal string excitations. However, the implementation of the flat-histogram approach [50] in conjunction with the diag-MC, carried out in the present work, allows us to extract the single-hole Green's function G k (τ ) with greater accuracy in a wide-range of imaginary time and, in addition, at long imaginary-time scales. This improves the accuracy of the spectral function at relatively low frequencies which enables us to resolve these multiple peaks in the spectral function due to the internal string excitations.
The paper is outlined as follows. In section 2 we discuss how we have applied the flat-histogram idea to the diag-MC method both when we treat the imaginary time τ as a parameter at a fixed value, and when we include τ as a dimension of the sampling space of the Markov process. In section 3 we compare the results obtained with NCA with the results obtained with diag-MC. In section 4 we compare the results obtained with diag-MC for the single-hole dispersion, the bandwidth and the quasiparticle residue as a function of J/t to those obtained with NCA. In section 5 we carry our analytic continuation to real-time to find the spectral functions. Section 6 presents our conclusions.

Model and method
We will use the diag-MC method introduced in references [45,46] as was implemented in reference [50]. Here we will briefly discuss the basic approach as we apply it to the single hole motion in the linearized version of the t-J model within the spin-wave approximation to obtain a polaron-like Hamiltonian [7,12,13] where b † q is the Bogoliubov spin-wave creation operator, ω(k) is the spin-wave dispersion of the square lattice quantum antiferromagnet and a † k is the hole creation operator. Also the g( k, q) is the coupling of the hole to spin waves and u q and υ q are the coefficients of the Bogoliubov transformation. For details of the derivation of this simplified version of the t-J model, which remains a non-trivial quantum many-body problem, as well as for the definitions of the operators and the expression for ω k , u q , v q and g( k, q), the reader is referred to references [7,12,13]. Furthermore, we do not expect that there are any significant numerical differences between the t-J model and that of equation (1), because we know that the spin-wave approximation is remarkably accurate for the antiferromagnetic Heisenberg model on the square lattice [2]. However, it is possible that this exclusion of the non-linear spin-wave interaction terms may miss some qualitative new physics and it is a totally inadequate approximation for frustrating lattices. Nevertheless, the model defined by equation (1) is a non-trivial model, which captures the essence of the physics that we attempt to understand, and, therefore, it deserves scrutiny.
It is rather well-known that the diag-MC technique [45,46] is a Markov process in a space defined by all the terms (or diagrams) which appear in perturbation theory. The perturbation expansion of the imaginary-time single-particle Green's function can be schematically written as follows: Here O ( k) n represents the sum of all the diagrams of order 2n, and O ( k) 0 (τ ) = G (0) k (τ ) is the Green's function in zeroth order, and λ is a variable which labels the contribution T ( k,λ) n of a particular diagram (term) of order 2n. As n increases, the number of integration variables increases in a similar manner. Notice that the nth term in the case of the t-J model is actually the 2n order, because there are no odd orders in the expansion and n is the number of spin-wave-propagators contained in the diagrams.
In the present paper, we will apply the diag-MC method in two different ways: (a) First, in our Markov chain, we will use the imaginary-time τ as a parameter with fixed given value. In diag-MC the random walk makes a series of transitions between states {n, λ, Through such a Markov process the entire series of terms is sampled. This process generates a histogram which represents the number of times N n (τ ) the order n appeared in the Markov process. Since we can compute a low order diagram analytically, say for example the zeroth order O ( k) 0 (τ ), the absolute value of all other orders is computed as follows: (b) In the second approach, τ becomes one of the dimensions of the sampling space of the Markov process as is the case in the usual implementation of the diag-MC method [45,46]. In this case the random walk makes a series of transitions between states {n, λ,  is the perturbation order, λ is a particular diagram, and i (or i ) is the label of a particular imaginary time interval (τ i − δτ /2, τ i + δτ /2) (where δτ = τ i+1 − τ i ). Namely, in this case we sample the histogram of G k (τ ) by including the imaginary time as an extra dimension of the sampling space. In this case, the value G k (i) of the histogram G k (τ i ) in the ith τ -interval is found by using the known value of G k (0) in the first interval, namely, where M i is the number of occurrences of the MC random walk in the ith interval.
As discussed in reference [50] we can obtain significant improvement when calculating the imaginary-time Green's function by applying flat-histogram techniques. These improvements allows us to extend the range of the imaginary time τ to larger values. We find that essential in order to obtain accurate values of the low lying hole-excitations. In the simulations of the present paper where we keep τ fixed, i.e. case (a) above, we will make the histogram of O ( k) n flat (equation (6)). To make the histogram flat in these calculations we used the multicanonical algorithm [51]. While in the simulations of case (b) above, we make to the histogram of G k (τ ) flat (equation (7)). To make the histogram flat in these calculations we used the Wang-Landau algorithm [52,53].
We have carried out calculations in which all diagrams are allowed to be sampled by the Markov process. In order to compare with NCA we have also carried out calculations in which only the diagrams which are included in the NCA are sampled. The computer code for the latter calculations is actually more complicated than that in which we allow all possible diagrams to be sampled. In figure 1 we compare our results obtained by the diag-MC by restricting the sampling space to only NCA diagrams to those obtained  by the real-time perturbation expansion incorporated in reference [13]. This provides a very strict test that our computer program is correct.   spin-wave propatators peaks at some value of n = n max (τ ). As the value of τ increases the value of n max (τ ) also increases almost linearly with τ . As can be seen in figure  has a peak at a very different value of n = n max and, in fact, n max agrees very well with the value where O ( k) n has the peak when only NCA diagrams are sampled. In addition, the amplitude of the distribution also agrees very well with the NCA (blue line). Figure 3 compares the results obtained for J/t = 0.4 and J/t = 0.7 and both correspond to k = (0, 0). Notice that while the agreement is not as perfect as in the case of k = (π/2, π/2) shown in figure 2 (the case of J/t = 0.7 for this case is not shown because the agreement is also excellent), it is still very good.

Accuracy of NCA
Therefore, we conclude that in all cases the positive and negative contributions largely cancel each other and the total contribution (red circles) which peaks at a smaller value of n agrees remarkably well with the results of NCA (blue solid line). Figure 4 presents our results for the single hole Green's function G k (τ ) as a function of τ for k = (0, 0) and k = (π/2, π/2) and J/t = 0.4 and J/t = 0.7 and compares them with those obtained within NCA. Notice that the positive and negative contributions G (±) k (τ ) to G k (τ ) are largely canceling each other at large τ yielding a net G k (τ ) in close agreement with the results of NCA. This is a strong indication that the low energy excitations of the t-J model are correctly described by NCA.

Quasiparticle properties
To obtain the quasiparticle energy-wavenumber dispersion E k and its residue Z k we fit the imaginary-time Green's function G k (τ ) at very long imaginary time as The quasiparticle properties such as the energy dispersion, the quasiparticle residue and the bandwidth W = E (0,0) − E (π/2,π/2) as a function of J/t are shown in figures 5-7 respectively. Notice that the agreement

Spectral functions-string excitations
We have applied the stochastic analytical inference technique discussed in reference [54], a technique similar to the maximum entropy method [55] to carry out analytic continuation of our diag-MC results on G k (τ ) from imaginary-time to real-time. This approach requires a default model D k (ω) as input. In addition, it requires a data table G  A( k, ω) is determined by the method described in references [54,56]. We have carried out two different calculations: (a) one using as default model the result of the NCA, and (b) one using as default model a distribution which is flat everywhere in an interval of ω min < ω < ω max . The default model summarizes any any information that we might have about the features of A( k, ω), which we want to test. If we use the flat distribution as the default model, we implicitly declare no such a priori knowledge. When we use the result of the NCA as the default model, it means that we use the NCA to define the prior probability of the spectrum A( k, ω). This bias will influence the posterior probability of A( k, ω) given the input data in the table of G (j) k (i). The Bayes' theorem relates these probabilities and the likelihood function which is used to sample A( k, ω). The reader is referred to reference [54] for details of this method.
The results of the analytic continuation are presented in figure 8 for J/t = 0.3 (which is close to the physical region of the t-J model) and in figure 9 for J/t = 0.7. Notice that, in general, the spectral function is in excellent agreement with NCA when the result of the NCA is used as default model. This means that the result of the NCA is consistent with the correct A( k, ω). The fact that the results of the analytic-continuation is very close to the results of the NCA, when used as the default model, it implies that the multiple peaks due to the internal excitations of the string, as likely to be present in the t-J model.
On the other hand, when using as default model the flat distribution, namely, assuming that we have no a priori information on what the correct A( k, ω) is, we obtain the results shown in figures 8(b) and (d) for J/t = 0.3 and in figures 9(b) and (d) for J/t = 0.7, where the lowest energy quasiparticle peak is in excellent agreement with the NCA with respect to its location and its spectral weight. In addition, at higher values of ω there is a broad distribution centered around the average value of the string excitations. This is also in agreement with the results of the NCA. The accuracy of the analytic continuation procedure is not good enough to resolve the particular peaks when such a flat distribution is used as the default model.
Let us summarize the meaning of these results. If we use the flat distribution as the default model, we implicitly declare no such a priori knowledge. In this case the accuracy of the inversion is such that only a broad peak can be resolved above the main peak. If we use the NCA (or SCBA) as a default model, the result of the inversion is that the default model is very close to the output. The fact that, in this case, the χ-square of the fit to the data-table obtained by diag-MC is lower when the NCA is used as the default model than that obtained when the flat-distribution is used as the default model, supports the NCA prediction, and, therefore, the suggested string-picture.
Given that the analytic continuation and inversion from G k (τ ) to A( k, ω) (for real values of ω) requires high-degree of accuracy, we would like to discuss the relative accuracy of the present method as compared to previous MC methods. We have published previous articles (see for example our reference [50] and references therein) about the accuracy of the present MC method for the one-particle Green's function in related problems. We have shown that our method which implements flat histogram ideas [51][52][53], because of the re-weighting used in the multicanonical and in the Wang-Landau algorithm, can yield much more accurate results than previous methods for G k (τ ) at long imaginary-time. This is crucial when calculating the low-energy features of the spectral function. When t/J is large, the string excitations are relatively low-energy excitations contributing with relatively small weight in the spectral function. Therefore, our method by more accurately capturing the long imaginary-time behavior of G k (τ ) provides the information needed by the SAI method to yield an accurate spectral function.

Conclusions
We have applied the diag-MC method [18,45,46] using the flat histogram method [50] to obtain the imaginary-time single-hole Green's function in the simplified version of the t-J model as described by the Hamiltonian [7] given by equation (1) in which the hole moves in a fluctuating but antiferromagnetically ordered background of spins. Motivated by recent experimental investigations of ultra-cold atoms in optical lattices [37][38][39][40][41] where strings were found to play a significant role in the dynamics of the Hubbard model simulation, our goal was to examine whether or not the 'string'-picture [3-6, 8, 9, 12, 13] successfully describes the motion of a hole in a quantum antiferromagnet.
We find that the NCA is qualitatively and quantitatively very accurate. First, the contribution of diagrams O n (τ ) which contain n spin-wave propagators, and for a given value of τ , peaks at some value of n = n max (τ ). As the value of τ increases the value of n max (τ ) also increases almost linearly with τ . When the diag-MC process is allowed to sample all the diagrams, there are positive and negative contributions to any given O n which, for large τ peak at approximately the same value of n = n ± max (figures 2 and 3). For large values of τ , the positive and negative contributions to O n largely cancel each other out producing a net contribution which peaks at a very different value of n = n max . In fact, n max agrees very well with the value where O n has its peak when only NCA diagrams are sampled. In addition, the amplitude of the total O n obtained by diag-MC is in very good agreement to that obtained when the Markov process is restricted to sum only those diagrams included in the NCA. Furthermore, the Green's function G k (τ ) obtained by diag-MC (figure 4) has positive and negative contributions G ± k (τ ) which for large enough values of τ cancel each other out to yield a net value of G k (τ ), which is several orders of magnitude smaller than either G ± k (τ ), and which is in very good agreement with the results of a simulation where we only sample the NCA diagrams (illustrated in figure 4).
The physical results and conclusions agree very well with the details of the string-picture. In the physical region of the t-J model the single-hole dispersion obtained using diag-MC is within less than a few percent in agreement with that obtained by NCA. A similar statement of accuracy can be made for the bandwidth as a function of J/t and the quasiparticle residue or spectral weight.
The accuracy of the NCA found by carrying out the calculation by means of the diag-MC method, strongly suggests that the results found using NCA regarding the string-excitations [13] and the string picture are also correct. Finally, in order to provide more direct evidence for the validity of the string-picture which emerges within the NCA, we used the results obtained for the single-hole Green's function in imaginary time and carried out analytic continuation to real time. We found that our results of the analytic continuation are consistent with those obtained by the NCA and the contribution to the spectral weight of the string excitations was clearly found to be in agreement with the results of the NCA, thus, validating further the 'string' picture.