Entanglement entropy of the diamond chain

[HamidArianZad]

Dear Jie,

I would like to investigate the bipartite entanglement entropy of the spin-1 diamond chain in the absence of the magnetic field when conserved total sz and k=0 are assumed. I have encountered some issues and I would be grateful to you if help me to solve them. I extended the code as:

# Hamiltonian parameters:
J1 = 1.0
J2 = J1
J3 = 2.8
J4 = J1
J5 = J1
L = 2
N = 3*L
sz_min = 0
sz_max = 2*N
# -------------------------------------------------------
all_energies = []
EntEntropy_arr = []
Size_EigVals = []
Size_EntEnt = []
mun_EigVals = 0
for i in sz_min:sz_max
    Size_EigVals_k = 0
    Size_EntEnt_k = 0
        # onsite Hamiltonian matrix
        O1 = spin((J1, "xx1"), (J1, "yy1"), (J1, "zz1"), D=3)
        O2 = spin((J2, "x1x"), (J2, "y1y"), (J2, "z1z"), D=3)
        O3 = spin((J3, "1xx"), (J3, "1yy"), (J3, "1zz"), D=3)
        h1 = O1 + O2 + O3
        C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1z1z11"), D=3)
        C2 = spin((J5, "11xx11"), (J5, "11yy11"), (J5, "11zz11"), D=3)
        h2 = C1 + C2
        Nj = diag( spin((1,"z11"), (1,"1z1"), (1,"11z"), D=3) + 3*I )
        total_z(s) = sum(Nj[sj+1] for sj in s)
        basis = TranslationalBasis(L=L, f = s -> total_z(s) == i, k=0, base=3^3)
        H = trans_inv_operator(h1, 1:1, basis) + trans_inv_operator(h2, 1:2, basis)
        vals, vecs = Array(H) |> Hermitian |> eigen
        EE = [ent_S(vecs[:,i], 1:2, basis) for i=1:size(basis,1)]
        Size_EigVals_k += size(vals, 1)
        Size_EntEnt_k += size(EE, 1)
        for sublist in vals
            for element in sublist
                push!(all_energies, element)
            end
        end
        for sublist in EE
            for element in sublist
                push!(EntEntropy_arr, element)
            end
        end
    push!(Size_EigVals, Size_EigVals_k)
    push!(Size_EntEnt, Size_EntEnt_k)
end

1. After running the code, for all values of the total sz I get correct eigenenergies, but all entropy values are zero.
2. I can not compute bipartite EE for 1:3 or 1:4 and so on. The code gives only bipartite EE for the case 1:2.

As I understood, the above code gives EE between two neighboring blocks 1 and 2 not two spins, right? If so, how can I compute EE of arbitrary pairs of spins?

[jayren3996]

Hi Hamid,

Apologies for the delayed response, as I've been occupied lately. Regarding the questions. Since the blocked system comprises only two sites, it's only capable of providing bipartite entanglement entropy between two blocks. This can be achieved using:

ent_S(vecs[:,i], 1:1, basis)

At present, unfortunately, I do not have a straightforward workaround within the EDKit functions to compute entanglement within a block.

[HamidArianZad]

Dear Jie,

Thank you very much, I could calculate entanglement entropy between two blocks by using expression ent_S(vecs[:, i], 1:1, basis). Please know that the code works only for the case k=0. could you please comment on the physical meaning of the ent_S(vecs[:, i], 1:1, basis) and why it works only for the case k=0? Please also know that in QuSpin, I need to consider the case here is a part of my QuSpin code:

for mz in range(sz_min, sz_max+1):
    basis = spin_basis_1d(L, pauli=False, S="1", m=mz / L, kblock=0)
    H = hamiltonian(static, dynamic, basis=basis, dtype=np.float64,
                    check_symm=False, check_herm=False, check_pcon=False,
                    )
    E, V = H.eigh()
    Energies.append(E)  # Append only the ground state energy
    ent = basis.ent_entropy(V, sub_sys_A=(0, 1, 2, 3, 4, 5), density=False, enforce_pure=True)["Sent_A"]
    EnEnt.append(ent)

I am not sure that it is possible to take in to account other values of kblock symmetry in QuSpin to calculate EE for the diamond chain, as well. I would be thankful to you if comment in this case as well.

[jayren3996]

Hi Hamid, The entanglement for k≠0 is possible. For this L=2 case, the other possible k is pi. In EDKit's convention, you can set k=1, and the following commands should return the entanglement entropy:

basis = TranslationalBasis(L=L, f = s -> total_z(s) == i, k=1, base=3^3)
H = trans_inv_operator(h1, 1:1, basis) + trans_inv_operator(h2, 1:2, basis)
vals, vecs = Array(H) |> Hermitian |> eigen
EE = [ent_S(vecs[:,i], 1:1, basis) for i=1:size(basis,1)]