Anticommutation and Jordan Wigner strings in connected correlation function

ITensor Julia Questions
1

ITensor team,

I’m curious about the necessity of Jordan Wigner strings, or a custom operator containing them, for the “connected” correlation function of a fermionic system, i.e. (AB) - (A)(B) where A and B are local operators and () denotes expectation value with respect to the ground state. I know that the correlation_matrix function can handle the calculation of (AB) correctly, but I’m curious wether using the expect function for the other two terms in the above expression with built in local operators could be leading to incorrect results.

Specifically, I’m trying to compute the correlation functions inside the square brackets in (2) and (4) in this paper: https://arxiv.org/pdf/1503.08188 and then track how they decay over distance, fixing the first site and then tracking the correlations over the lattice. Basic physics of the Extended Fermi-Hubbard Model say that with a Hamiltonian that indicates say, a Spin-Density Wave ordered ground state, the corresponding correlation function (equation (2) in the paper) should show power law decay over the lattice. The method above returns an exponential decay. Here’s a code that can reproduce the error:

using ITensors, ITensorMPS, Plots


let
    N = 100
    t1 = 1.0
    U = 8.0
    V = -1.0

    sites = siteinds("Electron", N; conserve_qns=true)

    os = OpSum()
    for b in 1:(N - 1)
        os -= t1, "Cdagup", b, "Cup", b + 1
        os -= t1, "Cdagup", b + 1, "Cup", b
        os -= t1, "Cdagdn", b, "Cdn", b + 1
        os -= t1, "Cdagdn", b + 1, "Cdn", b
        os += V, "Ntot", b, "Ntot", b + 1
    end
    os += V, "Ntot", N, "Ntot", 1
    os -= t1, "Cdagup", N, "Cup", 1
    os -= t1, "Cdagup", 1, "Cup", N
    os -= t1, "Cdagdn", N, "Cdn", 1
    os -= t1, "Cdagdn", 1, "Cdn", N
    for i in 1:N
        os += U, "Nupdn", i
    end
    H = MPO(os, sites)

    nsweeps = 5
    maxdim = [50, 100, 200, 400]
    cutoff = [1E-8]

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi0 = randomMPS(sites, state; linkdims=10)
    energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)

    corr_mat_sdw = correlation_matrix(psi, "Sz", "Sz")
    sdw_decay = []
    Sz1 = expect(psi, "Sz"; sites = 1)
    for i in 1:N
        value = corr_mat_sdw[1, i] - Sz1 * expect(psi, "Sz"; sites = i)
        if value != 0
            push!(sdw_decay, abs(value))
        end
    end

    scatter!(1:length(sdw_decay), sdw_decay, xscale = :log10, yscale = :log10)

end

where a log-log plot is used to show exponential decay more clearly.

2

At a conceptual level, the answer to your question about fermionic, connected correlators is that the operators you are considering (whether \hat{\sigma}^z or \hat{n}) are ‘bosonic’ meaning they do not carry any Jordan-Wigner string (because they are made up of an even number of \hat{c} and \hat{c}^\dagger operators). So there cannot be any issue due to the expect function, as you will be putting bosonic operators into expect anyway in both of your cases (both cases meaning equations (2) and (4) in the linked paper).

I would guess, more likely, that your wavefunction isn’t sufficiently accurate to reproduce the power-law behavior you expect. I notice that you are doing just 5 sweeps, which can be plenty for a gapped system with a short correlation length. But for a system with power-law correlations, these can take a while longer for DMRG to “figure out” meaning that more sweeps are called for. I would suggest a sweep plan more like the following:

nsweeps = 25
cutoff = 1E-8
maxdim = [2,2,2,2,2,8,8,8,8,8,20,20,20,20,20,50,50,50,50,50,100,100,200,200,400]
1 Like
3

Thanks so much, the the sweeping plan is making order parameter (2) work much better and display power-law decay in the correct region. However, I’m still having issues with other parameters, where the charge-density wave correlation function inside (4) is showing power-law decay inside the superconducting region with a periodicity of two sites (the phase diagram is given just above these equations in the linked paper). Again, my understanding of the Fermi-Hubbard model is that the value of this correlation function should decay exponentially in superconducting regions. Here’s another code that can reproduce that issue.

Excuse me if it’s inappropriate to ask a question this general in the same post, my intention is not to have anyone debug my code for me, but rather point out issues that result from a misunderstanding of the ITensor package.

using ITensors, ITensorMPS, Plots


let
    N = 100
    t1 = 1.0
    U = -2.0
    V = -0.3

    sites = siteinds("Electron", N; conserve_qns=true)

    os = OpSum()
    for b in 1:(N - 1)
        os -= t1, "Cdagup", b, "Cup", b + 1
        os -= t1, "Cdagup", b + 1, "Cup", b
        os -= t1, "Cdagdn", b, "Cdn", b + 1
        os -= t1, "Cdagdn", b + 1, "Cdn", b
        os += V, "Ntot", b, "Ntot", b + 1
    end
    #os += V, "Ntot", N, "Ntot", 1
    #os -= t1, "Cdagup", N, "Cup", 1
    #os -= t1, "Cdagup", 1, "Cup", N
    #os -= t1, "Cdagdn", N, "Cdn", 1
    #os -= t1, "Cdagdn", 1, "Cdn", N
    for i in 1:N
        os += U, "Nupdn", i
    end
    H = MPO(os, sites)

    nsweeps = 25
    cutoff = 1E-8
    maxdim = [2,2,2,2,2,8,8,8,8,8,20,20,20,20,20,50,50,50,50,50,100,100,200,200,400]

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi0 = randomMPS(sites, state; linkdims=10)
    energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)

    corr_mat_cdw = correlation_matrix(psi, "Ntot", "Ntot")
    cdw_decay = []
    Ntot1 = expect(psi, "Ntot"; sites = 1)
    for i in 1:N
        value = corr_mat_cdw[1, i] - Ntot1 * expect(psi, "Ntot"; sites = i)
        if value != 0
            push!(cdw_decay, abs(value))
        end
    end

    scatter!(1:length(cdw_decay), cdw_decay, xscale = :log10, yscale = :log10)

end