Confusion about the central charge extracted from Entanglement Entropy

DMRG and Numerical Methods
,
1

Hey,

I have a question about the central charge one can extract from the entanglement entropy of some critical model. From my understanding, and from the paper by Calabrese and Cardy
( Pasquale Calabrese and John Cardy 2009 J. Phys. A: Math. Theor. 42 504005) the entanglement entropy of a CFT on a finite slab of length l with open boundary conditions has the following asymptotic behavior

S(l) = \frac{c}{6}\log\left( \frac{2L}{\pi a} \sin \frac{\pi l}{L}\right) + \tilde{c}^\prime \, ,

taken from equation (39) from the mentioned paper. In this equation, c is the central charge, L is the
length of the system, a is some short range cutoff and \tilde{c}^\prime is some uninteresting constant.

One of the most simplest CFT is given by a spinless free fermion:

H = -t\sum_{j = 1}^N c_{j}^\dagger c_{j + 1}^{\vphantom \dagger} + \text{h.c.}

which should have a central charge of c = 1. The above Hamiltonian is defined on a lattice with the convention:

N = \frac{L}{a}\, \quad l = \frac{j}{a}\,,

i.e. all lengths are measured in units of the short distance cutoff a which is simply the lattice constant in that model.

The central charge can also be obtained by going to the bosonized version of this model, where it is simply given by a free compactified boson:

H \approx \frac{u}{2\pi}\int_0^L\! \text{d}x \,K(\pi\Pi(x))^2 + \frac{1}{K} (\nabla\varphi(x))^2

with K = 1 for the non-interacting spinless fermion.

The next non-trivial extension would be to consider two or more flavors of free Fermions, all described by some free Hamiltonian:

H = -t\sum_{f = 1}^{n_f}\sum_{j = 1}^N c_{j,f}^\dagger c_{j + 1, f}^{\vphantom \dagger} + \text{h.c.}

which should simply have a central charge of c_{n_f} = n_fc = n_f.

Since this is one of the most natural quantities to calculate in DMRG, I tested the spinless Fermion and the two flavor Fermion using ITensor. The code I used is given by:

sing ITensors, Plots

function chrd_dist(j,L)
    return 2*L/Ď€ * sin(Ď€/L * j) 
end

function entropy(psi)
    psic = copy(psi)
    N = length(psi)
    ent = Vector{Float64}(undef, N-1)
    for jj in 1:N-1
        orthogonalize!(psic, jj)
        U, S, V,spec = svd(psic[jj], (linkind(psic, jj - 1), siteind(psic, jj)))

        ent[jj] = mapreduce(+, 1:dim(S,1), init = 0) do n
            p = S[n,n]^2
            return -p*log(p)
        end
    end
    return ent, map(b -> chrd_dist(b, N), 1:N-1)
end

This I tested using the following model:

N = 200
sites = siteinds("Fermion", N)
ampo = OpSum()

t = 1
ÎĽ = 0.0
V = 0.0
for j=1:N-1
    ampo += -t,"Cdag", j, "C", j+1
    ampo += -t,"Cdag", j+1, "C", j
    ampo += V/4,"N + N - Id",j,"N + N - Id",j+1
end
for j=1:N
    ampo += -ÎĽ,"N",j
end

H = MPO(ampo, sites)

psi0  = randomMPS(sites, 10)

bnd = [10,10,20,20,50,50,100,100,200,200,200,200,200,200,250,250,250]
nsweeps = length(bnd)
noise  = vcat(repeat([0.1], nsweeps - 2), [0,0,0])
sweeps = Sweeps(nsweeps)
setmaxdim!(sweeps, bnd...)
setnoise!(sweeps, noise...)
setcutoff!(sweeps, 1E-10)

energy, psiof = dmrg(H,psi0, sweeps)

entf, chrdf = entropy(psiof)
pl = plot(log.(chrdf), log.(entf), label = false, seriestype = :scatter)
alg = (x, δ,β) -> δ * x + β 
c = 1
plot!(pl, x -> alg(x, c/6, -0.6),log.(chrdf), label = "c = $c", legend = :topleft)

xlabel!(pl, "d(l)")
ylabel!(pl, "S(l)")
title!("One free Fermion")


savefig("scaling_one_flavor.png")

display(pl)

which shows the correct expected behavior of a c=1 scaling, see left picture.

This behavior also shows by varying the interaction away from V = 0 as expected in the Luttinger Liquid regime.

However, doing the same calculation for a two flavor Fermion system shows some strange result.
Using the following code:

N = 100
t = 1.0


sites = siteinds("Electron", N)

ampo = OpSum()
for b in 1:(N - 1)
  ampo += -t, "Cdagup", b, "Cup", b + 1
  ampo += -t, "Cdagup", b + 1, "Cup", b
  ampo += -t, "Cdagdn", b, "Cdn", b + 1
  ampo += -t, "Cdagdn", b + 1, "Cdn", b
end
H = MPO(ampo, sites)


psi0 = randomMPS(sites, 10)

bnd = [10,10,20,20,50,50,100,100,200,200,400,400]
nsweeps = length(bnd)
noise  = vcat(repeat([0.1], nsweeps - 2), [0,0,0])
sweeps = Sweeps(nsweeps)
setmaxdim!(sweeps, bnd...)
setnoise!(sweeps, noise...)
setcutoff!(sweeps, 1E-10)

# Start DMRG calculation:
energy, psi2f = dmrg(H, psi0, sweeps)

entf, chrdf = entropy(psi2f)
pl = plot(log.(chrdf), log.(entf), label = false, seriestype = :scatter)
alg = (x, δ,β) -> δ * x + β
c = 2
plot!(pl, x -> alg(x, c/6, -0.2),log.(chrdf), label = "c = $c", legend = :topleft)
c = 1
plot!(pl, x -> alg(x, c/6, -0.0),log.(chrdf), label = "c = $c", legend = :topleft)

xlabel!(pl, "d(l)")
ylabel!(pl, "S(l)")
title!("Two free Fermion")

savefig("scaling_two_flavor.png")
display(pl)

shows that in that case the scaling is also given by an c=1 instead of the expected c=2 CFT, see right piture.

both

Is there something wrong or do I miss some fundamental stuff here?

Thanks for all replies in advance!

EDIT:

Sorry the labels of the axis are wrong, it should be \log(d(l)) on the x-axis, where d(l) = \frac{2L}{\pi a} \sin \frac{\pi l}{L} .

2

It’s an interesting question because I would also have guessed that your two-component model would show a central charge of c=2. Nothing about your code looks wrong that I can see – I would encourage you to check that other observables look ok, such as the density, since a system like this with two interleaved 1D chains can be really challenging to converge well in DMRG.

Something you might want to do is to use complementary, purely free-fermion techniques to check your entanglement results. It’s surprisingly easy to get the entanglement of a free fermion system just from the eigenvectors of the hopping Hamiltonian matrix. You just get these eigenvectors, select as many as the number of particles you have (view this as a rectangular matrix now), and square the resulting matrix to get the correlation matrix \langle c^\dagger_i c_j \rangle. Then if you carve out the upper [1:j,1:j] block and compute its eigenvalues – call these n_m – then the entanglement entropy is given by

S = \sum_n S_1(n_m)

where

S_1(n) = - n \log(n) - (1-n) \log(1-n)

So you could construct a two-component fermion chain as just a “hopping matrix” with second-neighbor hops only, and get its entanglement this way. Might be a good way to check your intuition about the central charge.

1 Like
3

Hey Miles,

thanks for the reply. Yes, I checked the convergence of the energy, density and correlation functions with the exact results. And I know that this is a super hard problem for DMRG, therefore I also checked that the result for the Entanglement entropy does not change by increasing the bond dimension, i.e. there was no qualitative change going from 200 to 400 bond dimension.
The code I used for calculating the exact properties is given by:

using LinearAlgebra

function hamiltonianNf(t, N,nf)
    tar = fill(-t, nf*N-nf)
    H   = zeros(Float64,nf*N, nf*N)
    idx = map(1:nf*N-nf) do jj
        CartesianIndex(jj, jj + nf)
    end
    @. H[idx] = tar
    return Symmetric(H)
end

function correlation_matrix(H)
    vals, vecs = eigen(H)
    idx_filled = vals .≤ 0
    idx_empty  = map(!, idx_filled)
    Eexact = sum(vals[idx_filled])
    W = vecs[:,idx_empty]
    
    return W*conj(transpose(W)), Eexact
end

function entanglement_entropy(corr_mat, jj,nf)
    cdagc = I - transpose(corr_mat)
    red_mat = @view cdagc[1:jj*nf, 1:jj*nf]

    S = mapreduce(+, eigvals(red_mat), init = 0) do n
        return abs(n) < 1E-13 || abs(1-n) < 1E-13 ? 0 : -n*log(n)-(1 - n)*log(1 - n)
    end
    return S
end

which also includes your formula for the entanglement entropy (the correlation matrix I was computing was \langle c_j^{\vphantom \dagger} c_k^\dagger \rangle, therefore the 1\!\!1 - Corr^T in the beginning of the entanglement_entropy function). nf denotes the number of fermionic flavors.

For the spinless fermion I find a ground-state energy difference of 10^{-8} with my choice of parameters
And a maximal relative difference in the density given by:

\text{max}\left( \frac{|n_{j, dmrg} - n_{j,exact}|}{n_{j,exact}} \right) \approx 10^{-9}

For the case of the two falvor fermion I find for the energy difference 10^{-4} and the error in the density (only comparing the total density):

\text{max}\left( \frac{|n_{j, dmrg} - n_{j,exact}|}{n_{j,exact}} \right) \approx 10^{-7}

Correlation functions are similar, however I just computed the maximal error instead of the relative error. In the spinless fermion case I got an absoulte error of 10^{-5} and for the two flavor case 10^{-4}.

So in terms of these observables I would say the algorithm converges, However the entropy calculated from your algorithm produces the expected c = 1 and c = 2 result.

exact_ent

4

That all sounds really good. So if I’m understanding correctly, the entanglement obtained from the correlation matrix eigenvalues is giving the (presumably correct) results you expected.

Then here is the puzzle: shouldn’t the entanglement from your MPS be exactly the same as obtained from the correlation matrix? If it isn’t, perhaps it’s just due to a bug in your code somewhere?

5

Yes, the correlation Matrix method is giving me the correct expected results for models. And I am with you, the MPS should give the same results. However, since all observables are correct, the only place to have a bug is in the calculation of the entanglement entropy itself. But I quite confident that that is the correct way of calculating the entropy for some MPS.

Now I am even more confused:D

6

Yes, exactly. I would suggest as a next step making a small test code that does more arbitrary free-fermion systems on some small lattices (just say 10 sites) and where you put some more arbitrary hoppings or on-site potentials, and check whether you can get your two methods of computing entanglement to agree. Then you might find the bug that way.

7

Mhm sounds reasonable. Will try that tomorrow or next week. Will keep you updated!

1 Like
8

Hey Miles,
sorry for the long time since the last reply. I re-implemented the code and it seems that in the old version was some stupid error. However, looking at the new pictures, everything works as expected… Thank you very much for the time!

In this picture I plotted the entanglement extracted from MPS for a chain of length 100, once with one
spinless fermion (on the left) and once with two spinless fermions on the right. I also compare this
with the exact result obtained using the single particle correlation matrix, and also plotted a
straight line according to the formula

y(x) = \frac{c}{6}x + \beta

against the logarithm of the chord distance:

d(l) = \frac{2L}{\pi}\sin\frac{\pi l}{L}

entanglement_scaling
entanglement_scaling1600Ă—800 60.3 KB

1 Like
9

Thanks for the nice followup! Glad it is working well now and that you found the issue.

1 Like