How to correctly compute two-site reduced density matrices of fermionic MPS?

ITensor Julia Questions
, ,
1

Hi everyone,

I am new to this amazing package and working on extracting two-site reduced density matrices from MPS. With the help of Is this a correct way to calculate one-site and two-site entanglement spectrum and entropy? - ITensor Julia Questions - ITensor Discourse, I figure out the way and benchmark it with theoretical formula from https://arxiv.org/abs/2310.15273,

\rho_A=\left[\begin{array}{cccc} \frac{1}{4}-m_z+\varrho_r^z & 0 & 0 & \varrho_r^x-\varrho_r^y \\ 0 & \frac{1}{4}-\varrho_r^z & \varrho_r^x+\varrho_r^y & 0 \\ 0 & \varrho_r^x+\varrho_r^y & \frac{1}{4}-\varrho_r^z & 0 \\ \varrho_r^x-\varrho_r^y & 0 & 0 & \frac{1}{4}+m_z+\varrho_r^z \end{array}\right] \tag{1}

where m_z=\left\langle S_{i_1}^z\right\rangle and \varrho_r^\alpha=\left\langle S_{i_1}^\alpha S_{i_i}^\alpha\right\rangle. Note that there are some constraint on the models. The minimal code for transversed-field Ising chain is as follows,

using ITensors, ITensorMPS

function TFIM_MPO(L::Int64, J::Float64, h::Float64)
    sites = siteinds("S=1/2",L)
    os = OpSum()
    for jc ∈ 1:L
        # transverse field
        os += -2.0*h,"Sz",jc
        # Ising interaction
        os += -4.0*J,"Sx",jc,"Sx",mod1(jc+1,L)
    end
    H = MPO(os,sites)
    # initialize the MPS
    state = [isodd(i) ? "Up" : "Up" for i in 1:L]
    psi0 = productMPS(sites, state)
    return sites, H, psi0
end

L = 8
J = 1.0
h = 0.5
sites, H, psi0 = TFIM_MPO(L, J, h)
nsweeps = 40
maxdims = [2,2,8,8,20,20,50,50,100,100,200,200,400]
cutoff1 = [1E-13]
noise = [1E-4, 1E-4, 1E-5, 1E-6, 1E-7, 1E-8, 0.0]
energy, psi = dmrg(H,psi0;nsweeps,maxdim=maxdims,cutoff=cutoff1, noise)

## two-site reduced density matrix
s1 = 2
s2 = 4

begin
    orthogonalize!(psi,s1)
    psidag = prime(dag(psi), linkinds(psi))
    rho_red = prime(psi[s1],linkinds(psi,s1-1))*prime(psidag[s1],sites[s1])
    # @show rho_red
    for i = s1+1:s2-1
        rho_red *= psi[i]
        rho_red *= psidag[i]
    end
    rho_red *= prime(psi[s2],linkinds(psi,s2))*prime(psidag[s2],sites[s2])
    C1 = combiner(sites[s1],sites[s2]; tags="s")
    C2 = combiner(sites[s1]',sites[s2]'; tags="s'")
    rho_red = C2*(C1*rho_red)

    ## Correlation function
    SS = zeros(3)
    for (i,comp) in enumerate(["Sx","Sy","Sz"])
        SS_os = OpSum()
        SS_os += 1,comp,s1,comp,s2
        SS_mpo = MPO(SS_os, sites)
        SS[i] = inner(psi',SS_mpo,psi)
    end
    mzs1_os = OpSum()
    mzs1_os += 1,"Sz",s1
    mzs1_mpo = MPO(mzs1_os, sites)
    mzs1 = inner(psi',mzs1_mpo,psi)
    @show rho_red
    @show (0.25+mzs1+SS[3]) ≈ rho_red[1,1]
    @show (SS[1]-SS[2]) ≈ rho_red[1,4]
    @show (0.25-SS[3]) ≈ rho_red[2,2]
    @show (SS[1]+SS[2]) ≈ rho_red[2,3]
    @show (SS[1]+SS[2]) ≈ rho_red[3,2]
    @show (0.25-SS[3]) ≈ rho_red[3,3]
    @show (SS[1]-SS[2]) ≈ rho_red[4,1]
    @show (0.25-mzs1+SS[3]) ≈ rho_red[4,4]
end

One can change the s1 and s2 and find that all the relations can be verified. The reduced density matrices are also consistent with my exact diagonalization calculations. I hope that this benchmark would be helpful.

However, my question lies in the fermionic scenario. There is also a similar formula for two-site fermionic reduced density matrix,

\rho_{A}=\left[\begin{array}{cccc} 1-2\overline{\mathrm{n}}+\varrho_{r}^{\mathrm{n}} & 0 & 0 & 0\\ 0 & \overline{\mathrm{n}}-\varrho_{r}^{\mathrm{n}} & g_{r} & 0\\ 0 & g_{r}^{*} & \overline{\mathrm{n}}-\varrho_{r}^{\mathrm{n}} & 0\\ 0 & 0 & 0 & \varrho_{r}^{\mathrm{n}} \end{array}\right] \tag{2}

where \overline{\mathrm{n}}=\left\langle\mathrm{n}_{i_1}\right\rangle, g_r=\langle c_{i_1}^{\dagger} c_{i_2}\rangle and \varrho_r^{\mathrm{n}}=\left\langle\mathrm{n}_{i_1} \mathrm{n}_{i_2}\right\rangle. The above method seems to not respect the anti-commutation relationship between fermions. In the following codes for simply free fermions,

using ITensors, ITensorMPS

function free_fermion_MPO(L::Int64, t::Float64, μ::Float64)
    sites = siteinds("Fermion",L)
    os = OpSum()
    for jc ∈ 1:L
        # chemical potential
        os += -μ, "n", jc
        # hopping
        os += -t, "c†", jc, "c", mod1(jc+1,L)
        os += -t, "c†", mod1(jc+1,L), "c", jc
    end
    H = MPO(os,sites)
    # initialize the MPS
    psi0 = random_mps(sites;linkdims=10)
    return sites, H, psi0
end

L = 8
t = 1.0
μ = 0.5
sites, H, psi0 = free_fermion_MPO(L, t, μ)
nsweeps = 10
maxdims = [2,2,8,8,20,20,50,50,100,100,200,200,400]
cutoff1 = [1E-10]
noise = [1E-4, 1E-4, 1E-5, 1E-6, 1E-7, 1E-8, 0.0]
energy, psi = dmrg(H,psi0;nsweeps,maxdim=maxdims,cutoff=cutoff1, noise)

## two-site reduced density matrix
s1 = 1
s2 = 2

begin
    orthogonalize!(psi,s1)
    psidag = prime(dag(psi), linkinds(psi))
    rho_red = prime(psi[s1],linkinds(psi,s1-1))*prime(psidag[s1],sites[s1])
    # @show rho_red
    for i = s1+1:s2-1
        rho_red *= psi[i]
        rho_red *= psidag[i]
    end
    rho_red *= prime(psi[s2],linkinds(psi,s2))*prime(psidag[s2],sites[s2])
    C1 = combiner(sites[s1],sites[s2]; tags="s")
    C2 = combiner(sites[s1]',sites[s2]'; tags="s'")
    rho_red = C2*(C1*rho_red)

    ## Green's function
    ns1_os = OpSum()
    ns1_os += 1, "n", s1
    ns1_mpo = MPO(ns1_os, sites)
    ns1 = inner(psi', ns1_mpo, psi)
    ns2_os = OpSum()
    ns2_os += 1, "n", s2
    ns2_mpo = MPO(ns2_os, sites)
    ns2 = inner(psi', ns2_mpo, psi)
    gr_os = OpSum()
    gr_os += 1, "c†", s1, "c", s2
    gr_mpo = MPO(gr_os, sites)
    gr = inner(psi', gr_mpo, psi)
    nnr_os = OpSum()
    nnr_os += 1, "n", s1, "n", s2
    nnr_mpo = MPO(nnr_os, sites)
    nnr = inner(psi', nnr_mpo, psi)
    @show rho_red
    @show 1-ns1-ns2+nnr ≈ rho_red[1,1]
    @show ns2-nnr ≈ rho_red[2,2]
    @show ns1-nnr ≈ rho_red[3,3]
    @show nnr ≈ rho_red[4,4]
    @show gr ≈ rho_red[2,3]
end

when s2 is not adjacent to s1 (e.g., s1=1,s2=3), the element rho_red[2,3] is not equal to the Green’s function gr, while the diagonal elements are good.

I’m quite certain that the issue lies in neglecting the sign change that occurs when fermions are exchanged. There are two reasons:

  1. When s2=s1±1, all the elements are good.
  2. At the beginning, the reduced density matrix was consistent with my own exact diagonalization codes. But soon I realized that both of them were problematic. I solved the discrepancy with Eq. (2), in my ED code, by adding the sign caused by partial trace. For example, \langle B_1 B_2|A_1B_1A_2B_2\rangle=(-1)^{n_{B_1}*n_{A_2}}\langle B_1 B_2|A_1A_2B_1B_2\rangle (see my last commit, if interested).

So I would like to ask what is the correct way to obtain two-site reduced density matrix of fermionic MPS.

Since the code works for s2=s1±1, I also tried to use ITensors.ITensorMPS.movesite before doing partial trace, like:

s1=1
if abs(s2-s1) == 1
    rho_red = reduced_density_matrix(psi, s1, s2)
else
    psi2 = movesite(psi, s2=>2)
    rho_red = reduced_density_matrix(psi2, s1, 2)
end

but that give me the same result as before. So maybe a related question is whether movesite respects the anti-commutation relation between fermions.

Thank you for reading!
Fo-Hong

1 Like
2

Hi Fo-Hong,
Thanks for the question. I believe the reason your code is not giving the correct results is that you are not including the Jordan-Wigner (J-W) “string” operators which must be attached to any fermionic operators, such as "c" or "c†" when working with them “manually” such as computing correlation functions. By “manually” I mean when simply using tensor contractions versus using one of our provided functions such as expect or correlation_matrix, which put in the J-W strings automatically. (Also OpSum puts in the J-W strings too.)

So my recommendation to you is to use the correlation_matrix function to measure your correlation functions involving "c" or "c†". Would that work for you?

Best,
Miles

3

Hi Miles,
Thanks for your advice. I agree with your opinion. Indeed, I tried to place Jordan-Wigner strings between the two sites s1 and s2 for calculating off-diagonal elements, and the results are correct. However, when calculating the diagonal elements, the string should be avoided. So my current workaround for two-site RDM is to do a “partial trace” for the diagonal and off-diagonal elements separately.

    sites = siteinds(psi)
    orthogonalize!(psi,s1)
    psidag = prime(dag(psi), linkinds(psi))
    rho_red = prime(psi[s1],linkinds(psi,s1-1))*prime(psidag[s1],sites[s1])
    rho_red_diag = rho_red
    # @show rho_red
    for i = s1+1:s2-1
        rho_red *= psi[i]
        JW = op("F", sites, i)
        rho_red *= JW
        rho_red *= prime(psidag[i],sites[i])
        rho_red_diag *= psi[i]
        rho_red_diag *= psidag[i]
    end
    rho_red *= prime(psi[s2],linkinds(psi,s2))*prime(psidag[s2],sites[s2])
    rho_red_diag *= prime(psi[s2],linkinds(psi,s2))*prime(psidag[s2],sites[s2])
    for i ∈ 1:2, j ∈ 1:2
        rho_red[i,i,j,j] = rho_red_diag[i,i,j,j]
    end

In fact, this workaround is no different from directly using Eq. (2) to calculate the RDM via correlation functions. I have not yet figured out a systematic method to perform a partial trace on fermionic MPSs. Furthermore, I wonder if the partial trace for fermions in MPS, in a generic bipartition, is still an open issue.

Yours sincerely,
Fo-Hong