Best Way to Compute Structure Factors

DMRG and Numerical Methods
1

I’m interested in computing structure factors to look for various quantum phases of matter. These take the following form, where O_i is our local order parameter (sometimes 2-site operators):

M(Q)=\sum_{i,j} \langle \psi | O_i^\dagger O_j |\psi\rangle e^{Q(i-j)\sqrt{-1}}

Would it be most efficient to loop through i and j, computing the correlation term \langle \psi | O_i^\dagger O_j |\psi\rangle each time and performing the sum? In which case I either use the correlation matrix function for single site order parameters, or compute explicitly with the inner product function for two site order parameters.

Or would it be faster to first define the following summed operator and take its expectation value for the ground state found by DMRG?

O_{\text{tot}}(Q)=\sum_{i,j} O_i^\dagger O_j e^{\mathcal{i}Q(i-j)}

Thank you!

2

If you’re only interested in a specific value of Q, or a small number of Q’s, I think it might be faster to make the MPO you described, once for each Q you want.

If instead you are interested in all values of Q, then it would be faster to compute C_{ij} = \langle O^\dagger_i O_j \rangle then Fourier transform the matrix C_{ij}.

If the operators act on a single site, then using the correlation_matrix function is a great way to go and it does efficient caching of partially formed state-operator overlaps. I don’t believe that function supports multi-site operators though – what is the type of operators you have?

3

Since I’m mostly approaching this from a numerical standpoint without analytical predictions about the phases we’d expect to find or at what Q it would be present, I’d like to evaluate it for all Q (these correspond to crystal momentum, so there’s finitely many).

As for the types of multisite operators, they’re composed of electron creation/annihilation operators defined only on adjacent sites. For example, I have operators like the following:

O_{i,\text{extended singlet}} = c_{i,\uparrow}c_{i+1,\downarrow}+c_{i+1,\uparrow}c_{i,\downarrow}

O_{i,\text{current}} = \sum_\sigma (c_{i,\sigma}^\dagger c_{i+1,\sigma} - c_{i,+1\sigma}^\dagger c_{i,\sigma})