Skip to content

Entanglement Entropy Ground-State XXZ

Author: Rafael Soares

In this example, we compute the entanglement entropy[1] of the ground state of the spin \(-\frac{1}{2}\) XXZ chain, described by the Hamiltonian: $$ \mathcal{H} = J \sum_{n=0}^{N-1} \left(S^x_n\cdot S^x_{n+1} + S^y_n\cdot S^y_{n+1} + \Delta S^z_n\cdot S^z_{n+1}\right). $$

The algorithm follows these steps:

  1. Obtain the ground state \(|\psi_0\rangle\) using the Lanczos algorithm.

  2. Construct the reduced density matrix for the region with the first \(\ell\) spins by tracing out the complementary degrees of freedom: $$ \rho_{\ell} = \text{Tr}_{\bar{\ell}} \left(|\psi_0\rangle\langle\psi_0| \right). $$

  3. Compute the entanglement entropy from the reduced density matrix: $$ S_\ell = -\text{Tr}\left( \rho_\ell \ln \rho_\ell \right), $$

which is done by diagonalizing \(\rho_\ell\).

Image title

As shown in the Figure above, the entanglement entropy of the XXZ chain ground state for two distinct values of \(\Delta\). In the critical phase, when \(\Delta \in (-J,J)\) [1], the entanglement entropy follows a logarithmic scaling with the subsystem size, as predicted by (1+1)-dimensional conformal field theory. Specifically, we extract the central charge, \(c\), by fitting the entanglement entropy to the expression derived for a finite system[2]: $$ S_{\text{CFT}}\left(\ell,L \right) = \dfrac{c}{3} \ln \left( \dfrac{L}{\pi} \sin\left(\frac{\pi\ell}{L}\right)\right) + b_0, $$ where \(L\) is the total size of the system and \(b_0\) is a constant. Away from this regime the system is in gapped phase, the system is in a gapped phase, where the entanglement entropy does not depend on the subsystem size for \(1\ll\ell\ll L\),i.e, it follows the area-law behaviour.

#include <xdiag/all.hpp>

using namespace xdiag;

void push_to_ps(ProductState &ps1, ProductState &psA, ProductState &psB);
void Create_RDM(arma::cx_mat &RDM, arma::vec &vstate, Spinhalf &Block,
                Spinhalf &BlockA, Spinhalf &BlockB);
double getVnE(arma::cx_mat &RDM);

int main() try {
  int N = 16;         // chain length
  double Delta = 1.0; // anisotropy

  auto ops = OpSum();
  auto block = Spinhalf(N); // Spin1/2 Block full

  for (int i = 0; i < N; i++) {
    ops += Op("Exchange", {i, (i + 1) % N});
    ops += Delta * Op("SzSz", {i, (i + 1) % N});
  }

  auto [e0, gs] = eig0(ops, block); // get GS with Lanczos

  arma::vec vstate = vector(gs);

  arma::mat vne = arma::mat(
      N / 2, 2,
      arma::fill::zeros); // Array to store EE as a function of subsystem size

  for (int na = 0; na < N / 2;
       na++) // Compute the Entanglement entropy for different subsystems from
             // l=1 to l=N/2, as the entanglement entropy satisfies: S(l) =
             // S(N=l).
  {
    int nb = N - (na + 1);
    auto blockA = Spinhalf((na + 1)); // Spin1/2 Block full
    auto blockB = Spinhalf(nb);       // Spin1/2 Block full
    arma::cx_mat RDM =
        arma::cx_mat(dim(blockA), dim(blockA), arma::fill::zeros);
    Create_RDM(RDM, vstate, block, blockA, blockB);
    vne(na, 0) = (na + 1);
    vne(na, 1) = getVnE(RDM);
  }

  // Construct the filename
  std::string flstring = "EE_XXX_model.Nsites." + std::to_string(N) + ".Delta" +
                         std::to_string((double)Delta) + ".outfile.h5";

  auto save_fl = FileH5(flstring, "w!");
  save_fl["entanglement_entropy"] = vne;
  return 0;
} catch (Error e) {
  error_trace(e);
}

void push_to_ps(ProductState &ps1, ProductState &psA, ProductState &psB) {
  for (auto sub : psA) {
    ps1.push_back(sub);
  }
  for (auto sub : psB) {
    ps1.push_back(sub);
  }
}

void Create_RDM(arma::cx_mat &RDM, arma::vec &vstate, Spinhalf &Block,
                Spinhalf &BlockA, Spinhalf &BlockB) {
  int i = 0;
  int j = 0;
  for (auto p1 : BlockA) {
    j = 0;
    for (auto p2 : BlockA) {
      for (auto p3 : BlockB) // loop over all configurations in subspace that we
                             // are interested in tracing.
      {
        auto p1aux = ProductState();
        auto p2aux = ProductState();
        push_to_ps(p1aux, p1, p3);
        push_to_ps(p2aux, p2, p3);
        int id1 = index(Block, p1aux);
        int id2 = index(Block, p2aux);
        RDM(j, i) += vstate(id2) * std::conj(vstate(id1));
      }
      j++;
    }
    i++;
  }
}

double getVnE(arma::cx_mat &RDM) {
  // Compute the Entanglement entropy
  double vne = 0.0;
  arma::vec pvals;
  arma::cx_mat eigvec;
  arma::eig_sym(pvals, eigvec, RDM); // FULL ED in the RDM

  for (int p = 0; p < pvals.n_elem; p++) {
    if (pvals(p) > 1e-14)
      vne -= pvals(p) * log(pvals(p));
  }
  return vne;
};

using LinearAlgebra
using XDiag
using HDF5
using Printf

function push_ProductState!(ps1, psA, psB)
    for sub in psA
        push!(ps1, String(sub))
    end
    for sub in psB
        push!(ps1, String(sub))
    end
end

function Create_RDM(RDM, vstate, Block, BlockA, BlockB) #Function to construct the RDM
    for (i, p1) in enumerate(BlockA)
        for (j, p2) in enumerate(BlockA)
            for (k, p3) in enumerate(BlockB)
                p1aux = ProductState()
                p2aux = ProductState()
                push_ProductState!(p1aux, p1, p3)
                push_ProductState!(p2aux, p2, p3)
                id1 = index(Block, p1aux)
                id2 = index(Block, p2aux)
                RDM[j, i] += vstate[id2] * conj(vstate[id1])
            end
        end
    end
end

function getVnE(RDM) #Function to compute the entanglement entropy.
    vne = 0
    p = eigvals(Hermitian(RDM))#FULL ED in the RDM
    for i in p
        if i>1e-14
            vne -= i*log(i)
        end
    end
    return vne
end

function main()
    N = 16 # chain length
    ops = OpSum()
    block = Spinhalf(N) # Spin1/2 Block full
    Delta = 1.0
    for i in 1:N
        ops += Op("Exchange", [i, mod1(i + 1, N)])
        ops += Delta*Op("SzSz", [i, mod1(i + 1, N)])
    end

    e0, gs = eig0(ops, block)# get GS with Lanczos  
    vstate = vector(gs)
    vne = zeros(Float64,(length(1:div(N,2)),2))
    for Na in 1:div(N,2)
        Nb = N - Na
        blockA = Spinhalf(Na) # Spin1/2 Block full
        blockB = Spinhalf(Nb) # Spin1/2 Block full
        RDM = zeros(ComplexF64, (dim(blockA), dim(blockA)))
        Create_RDM(RDM, vstate, block, blockA, blockB)
        vne[Na,1] = Na
        vne[Na,2] = getVnE(RDM)
    end
    filename = @sprintf("EE_XXX_model.Nsites.%d.Delta.%d.outfile.h5", N,Delta)
    h5open(filename, "w") do file
        write(file, "vne", vne)
    end 
end

main()

references

[1] Laflorencie, Nicolas, Quantum entanglement in condensed matter systems, Physics Reports 646, 1-59 (2016)

[2] Calabrese, Pasquale and Cardy, John, Entanglement entropy and conformal field theory, Journal of Physics A: Mathematical and Theoretical 42,50 (2009)