Calculating asymptotic behavior of a correlation function

[conana]

I was hoping to get some assistance in reproducing a calculation from https://arxiv.org/abs/0803.1292 (https://journals.aps.org/pra/abstract/10.1103/PhysRevA.78.012304 for the published version).

A certain four-spin correlation function for Kitaev's spin-1/2 model on the honeycomb lattice is defined in Eq. 28, and from Eq.'s 29 and 30 one arrives at the final expression for the correlation function as

$$C(\mathbf{r}_1, \mathbf{r}_2) = \frac{1}{N^2} \sum_{q,q'} \cos{[(\mathbf{q}-\mathbf{q}')\cdot(\mathbf{r}_1-\mathbf{r}_2)]} \frac{\Delta_q \Delta_{q'} - \epsilon_q \epsilon_{q'}}{E_q E_{q'}},$$

where N is the number of sites in the lattice, the summations are over momenta in the first Brillouin zone, $\mathbf{r}_1, \mathbf{r}_2$ denote unit cell positions, $\Delta_q = \sin{q_1} + \sin{q_2}$, $\epsilon_q = \cos{q_1} + \cos{q_2} + 1$, and $E_{q} = \sqrt{\epsilon_q^2 + \Delta_q^2}$ with $q_1, q_2$ being the components of $\mathbf{q}$ in the basis of reciprocal lattice vectors.

Arriving at this expression is not a problem. However, the authors go on to argue that for large $|\mathbf{r}_1 - \mathbf{r}_2|$, the correlation function decays algebraically as
$$C(\mathbf{r}_1, \mathbf{r}_2) \sim \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|^4}.$$
In arriving at this expression, the authors state "[...] the denominator in Eq. (30) has two zero points, which are of order 1/N in the large N limit. Their contribution causes the summation to be finite in the thermodynamic limit. Then using the stationary phase method, we can evaluate the exponents of the correlation function at long distance to be 4." Eq. 30 refers to (essentially) the above expression for the correlation function involving the summations over the Brillouin zone, and the two zero points of the denominator correspond to the two Dirac points in the fermionic quasiparticle spectrum.

It is the calculation of this asymptotic behavior which I am trying to reproduce, but I seem to be unable to do this with the information given. It seems there are no momenta at which the phase is stationary due to the phase being a linear function of momentum. So then I thought, perhaps the authors intend to imply that the phase is never stationary and, thus, the summation would vanish at all points in the Brillouin zone EXCEPT for at the Dirac points where the denominator vanishes, yielding a singularity in the summand. However, the numerator also vanishes at these points and, in fact, the limit of the summand at these points is finite.

Having run into only dead-ends with the stationary phase method, I though about approximating the sums as integrals in the $N \rightarrow \infty$ limit and performing successive integration by parts to extract the asymptotic behavior. (For sake of simplicity, in the following I take $\mathbf{r} = \mathbf{r}_1 - \mathbf{r}_2$ to be parallel to one of the lattice vectors). However, I seem to be getting terms proportional to
$$\frac{\sin^2{[\pi r]}}{r^{2n}}, \qquad n\in \mathbb{Z}^+,$$
where the lattice spacing is taken to be unity. These terms all vanish as $r$ is an integer (it denotes a distance along one direction of the lattice).

It seems clear to me that I am not understanding the spirit of the calculation/am making math errors as all of my attempts suggest that the correlation function vanishes. The numerics (of the authors and of my own) show, however, that the power law behavior described above is correct.

Any help would be greatly appreciated. Thank you in advance.

 

[DeathbyGreen]

I think they're just using an integral trick in order to approximate the Correlation function over a BZ. The function has a sum over q, which can be written as an integral in some limit. You could then rewrite the trigonometric terms as exponential functions and use this method: https://en.wikipedia.org/wiki/Stationary_phase_approximation. I think the integral will just produce a term proportional to x^-4. In this method you only evaluate the integral in an approximation, not exactly as it looks like you have done.

 

[conana]

Thank you for your response, DeathbyGreen.

Despite the authors claiming to have used the stationary phase approximation, I don't believe it applies as the cosine (or complex exponential) does not possesses a stationary phase. This is why I believe one should develop an asymptotic series with the boundary terms from successive integration by parts as outlined in, e.g., "Advanced Mathematical Methods for Scientists and Engineers" by Carl M. Bender and Steven A. Orszag.

Regarding the two zero points of the denominator mentioned in my original post, although the integrand remains finite at these points, I failed to realize initially that the integrand is discontinuous at these points and, thus, I was not correctly taking them into account when evaluating the boundary terms of the integration by parts.

The idea I then pursued was to express the integrand in cylindrical coordinates relative to the ill-behaved points and integrate by parts as (keeping in mind I also have to integrate everything with respect to q' as well)

$$\begin{align} \int_{BZ} d^2q \exp{[i \mathbf{q}\cdot\mathbf{r}]} \Big( \ldots \Big) &\rightarrow \lim_{\epsilon\rightarrow 0^+} \int_0^{2\pi} d\phi \int_\epsilon^\infty dq~q \exp{[i q r \cos{\phi}]} \Big( \ldots \Big) \nonumber\\ &= \lim_{\epsilon\rightarrow 0^+} \frac{1}{i r} \int_0^{2\pi} \frac{d\phi}{\cos{\phi}} ~q \exp{[i q r \cos{\phi}]} \Big(\ldots\Big) \Bigg|_{q=\epsilon}^\infty + \int d^2q \ldots \nonumber \end{align}$$

I extended my q-integration to infinity and figured I would just ignore that part of the boundary term due to the rapidly oscillating phase. Then my hope was to take the $\epsilon\rightarrow 0^+$ limit before performing the angular integral, yielding a simpler integrand. Here the discontinuity of the integrand at $\epsilon = 0$ should now show up as a $\phi$-dependent limit. Following this reasoning (and assuming my math beyond this point holds), the boundary term due to this first integration by parts vanishes in the limit due to an overall prefactor of $\epsilon$. I would then integrate by parts once more yielding two boundary terms, one of which vanishes when I take the limit, the other of which does not. I then perform this integration by parts procedure twice more for the integral over q' with it's phase factor $\exp{[-i\mathbf{q}'\cdot\mathbf{r}]}$.

Having done all that I wind up with
$$C(\mathbf{r}) \sim \frac{1}{r^4} \int_0^{2\pi} \frac{d\phi}{\cos^2{\phi}} \int_0^{2\pi} \frac{d\phi'}{\cos^2{\phi'}} \lim_{\epsilon\rightarrow 0^+} \lim_{\epsilon'\rightarrow 0^+} \Big( ... \Big)$$

At this stage I'd be happy if I could show that the remaining integration neither vanishes nor diverges. Unfortunately, it does vanish by anti-symmetry under $\phi\rightarrow\phi+\pi$.

Ok, if I'm making some simple mistake during all of my integrating by parts and limiting etc. then I need to be careful and find it. I don't expect anyone to go through all of that math for me hunting down my mistake. My question at this stage is more, is my method for approximating this integral conceptually sound? That is, integrating by parts, extending limits of integration and ignoring $q\rightarrow\infty$ contribution, etc.

Thanks again for reading and any help you all may have.