Entanglement Entropy Ground-State XXZ
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:
-
Obtain the ground state using the Lanczos algorithm.
-
Construct the reduced density matrix for the region with the firs \(\ell\) spins by tracing out the complementary degrees of freedom:
- Compute the entanglement entropy from the reduced density matrix:
which is done by diagonalizing \(\rho_\ell\).
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 <iostream>
#include <vector>
#include <string>
#include <xdiag/all.hpp>
using namespace xdiag;
void push_to_ps(xdiag::ProductState &ps1, xdiag::ProductState &psA, xdiag::ProductState &psB);
void Create_RDM(arma::cx_mat &RDM, arma::vec &vstate, xdiag::Spinhalf &Block, xdiag::Spinhalf &BlockA, xdiag::Spinhalf &BlockB);
double getVnE(arma::cx_mat &RDM);
int main()
try
{
unsigned N = 16; // chain length
double Delta = 1.0; // anisotropy
auto ops = OpSum();
auto block = Spinhalf(N); // Spin1/2 Block full
for (unsigned 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 (unsigned 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).
{
unsigned 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(xdiag::ProductState &ps1, xdiag::ProductState &psA, xdiag::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, xdiag::Spinhalf &Block, xdiag::Spinhalf &BlockA, xdiag::Spinhalf &BlockB)
{
unsigned i = 0;
unsigned 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);
unsigned id1 = index(Block, p1aux);
unsigned 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 (unsigned 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)