\(S=\frac12\) chain: Symmetries and Ground State Correlators
This example demonstrates how ground state expectation values can be computed for the spin-\(\frac{1}{2}\) XXX antiferromagnetic chain
with periodic boundary conditions.
This is done in three different ways, showcasing how the conservation of the overall number of up-spins as well as translation symmetry can be used to lower memory requirements.
-
The function
gs_correlator_simple()is the most basic way of obtaining the ground state and computing a correlator. It does not use any of the aforementioned symmetries of the system and therefore needs to store a full Hilbert space vector. Although being fast, this can cause memory issues. -
The function
gs_correlator_Sz_sym()exploits that \(H\) conserves the number of up-spins, decomposing the Hamiltonian into \(N+1\) separate blocks with a fixed number \(n_{up} = 0, 1, \ldots, N\) of spins pointing up. The Lanczos algorithm is then performed on each of the sub-blocks which may take longer but yields a ground state vector whose size is bounded by the size of the largest sub-block, substantially reducing the memory footprint. -
The function
gs_correlator_translation_sym()exploits that all nearest-neighbor bonds are equal and periodic boundary conditions are used. Thus, the Hamiltonian is invariant under translations and decomposes into \(N\) irreducible representations (irreps) labelled by lattice momentum \(2\pi/N \times k\) where \(k=0, 1, \ldots, N-1\). Again, the Lanczos algorithm is run on all momentum sub-blocks, effectively exchanging a longer runtime for less memory consumption. An important caveat when working with representations is that the operators inside expectation values must be symmetrized with respect to the symmetry group.
#include <xdiag/all.hpp>
#include <vector>
#include <complex>
#include <math.h>
using namespace xdiag;
// treat Hilbert space as a whole (memory consuming)
void gs_correlator_simple(int N, OpSum& H, OpSum& corr_op){
auto full_block = Spinhalf(N);
auto lanczos_res = eigs_lanczos(H, full_block);
auto e0 = lanczos_res.eigenvalues[0];
auto psi0 = lanczos_res.eigenvectors;
auto corr = inner(corr_op, psi0);
Log("Ground state energy: {:.12f}", e0);
Log("Ground state correlator: {:.12f}", corr);
Log("State vector length: {:d}", size(psi0));
};
// use Sz_tot conservation (less memory consuming)
void gs_correlator_Sz_sym(int N, OpSum& H, OpSum& corr_op){
double e0 = 0.0;
double gs;
State psi0;
// check blocks with fixed number of up-spins
for (int nup = 0; nup < N+1; nup++){
auto sym_block = Spinhalf(N, nup);
auto lanczos_res = eigs_lanczos(H, sym_block);
gs = lanczos_res.eigenvalues[0];
if (gs < e0){
e0 = gs;
psi0 = lanczos_res.eigenvectors;
};
};
auto corr = inner(corr_op, psi0);
Log("Ground state energy: {:.12f}", e0);
Log("Ground state correlator: {:.12f}", corr);
Log("State vector length: {:d}", size(psi0));
};
std::complex<double> C_N_character(int N, int k, int p){
return std::exp( ( (2 * k * p) * M_PI / (double)N ) * std::complex<double>(0,1.) );
};
// use translation symmetry (less memory consuming)
void gs_correlator_translation_sym(int N, OpSum& H, OpSum& corr_op){
// define shift by one site
std::vector<int64_t> T_perm (N, 0);
for (int i=0; i<N-1; i++){
T_perm[i] = i + 1;
};
auto T = Permutation(T_perm);
// define cyclic group
std::vector<Permutation> perms (N);
for (int i=0; i<N; i++){
perms[i] = pow(T, i);
};
auto C_N = PermutationGroup(perms);
// define irreps from character tables, labelled by momentum (2pi/N ×) k
std::vector<Representation> irreps (N);
std::vector<std::complex <double> > characters (N, 0.0);
for (int k = 0; k < N; k++){
for (int p=0; p<N; p++){
characters[p] = C_N_character(N, k, p);
};
irreps[k] = Representation(C_N, characters);
};
// check all translation symmetry blocks
double e0 = 0.0;
double gs;
State psi0;
for (int k = 0; k < N; k++){
auto sym_block = Spinhalf(N, irreps[k]);
auto lanczos_res = eigs_lanczos(H, sym_block);
gs = lanczos_res.eigenvalues[0];
if (gs < e0){
e0 = gs;
psi0 = lanczos_res.eigenvectors;
};
};
// now the operator must be symmetrized!
auto corr = inner(symmetrize(corr_op, C_N), psi0);
Log("Ground state energy: {:.12f}", e0);
Log("Ground state correlator: {:.12f}", corr);
Log("State vector length: {:d}", size(psi0));
};
int main() try {
// periodic N-site Heisenberg antiferromagnet
int N = 14;
auto H = OpSum();
for (int i = 0; i < N; i++) {
H += "J" * Op("SdotS", {i, (i+1) % N});
};
H["J"] = 1.0;
// define spin correlator at "half-chain" distance
auto corr_op = OpSum();
corr_op = Op("SdotS", {0, N/2 - 1});
// run different implementations
gs_correlator_simple(N, H, corr_op);
gs_correlator_Sz_sym(N, H, corr_op);
gs_correlator_translation_sym(N, H, corr_op);
} catch (Error e) {
error_trace(e);
};
using XDiag
function main()
# periodic N-site Heisenberg antiferromagnet
N = 14
H = OpSum()
for i in 1:N
H += "J" * Op("SdotS", [i, mod1(i+1, N)])
end
H["J"] = 1.0
# define spin correlator at "half-chain" distance
corr_op = OpSum()
corr_op += Op("SdotS", [1, N ÷ 2])
# run different implementations
gs_correlator_simple(N, H, corr_op)
gs_correlator_Sz_sym(N, H, corr_op)
gs_correlator_translation_sym(N, H, corr_op)
end
# treat Hilbert space as a whole (memory consuming)
function gs_correlator_simple(N::Int, H::OpSum, corr_op::OpSum)
full_block = Spinhalf(N)
lanczos_res = eigs_lanczos(H, full_block)
e0 = lanczos_res.eigenvalues[1]
psi0 = lanczos_res.eigenvectors
corr = inner(corr_op, psi0)
println("Ground state energy: $e0")
println("Ground state correlator: $corr")
println("State vector length = ", size(psi0))
end
# use Sz_tot conservation (less memory consuming)
function gs_correlator_Sz_sym(N::Int, H::OpSum, corr_op::OpSum)
e0 = Inf
psi0 = nothing
Nup_vals = 0:N
# check blocks with fixed number of up-spins
for nup in Nup_vals
sym_block = Spinhalf(N, nup)
lanczos_result = eigs_lanczos(H, sym_block)
gs = lanczos_result.eigenvalues[1]
if gs < e0
e0 = gs
psi0 = lanczos_result.eigenvectors
end
end
corr = inner(corr_op, psi0)
println("Ground state energy: $e0")
println("Ground state correlator: $corr")
println("State vector length = ", size(psi0))
end
# use translation symmetry (less memory consuming)
function gs_correlator_translation_sym(N::Int, H::OpSum, corr_op::OpSum)
# define shift by one site
T = Permutation(circshift(1:N, -1))
# define cyclic group
C_N = PermutationGroup([T^k for k in 0:(N-1)])
# define irreps from character tables, labelled by momentum (2pi/N ×) k
character_table = [ [C_N_character(N, k, p) for p in 0:(N-1)] for k in 0:(N-1)]
irreps = [ Representation(C_N, character_table[k]) for k in 1:N ]
# check all translation symmetry blocks
e0 = Inf
psi0 = nothing
for k in eachindex(irreps)
sym_block = Spinhalf(N, irreps[k])
lanczos_result = eigs_lanczos(H, sym_block)
gs = lanczos_result.eigenvalues[1]
if gs < e0
e0 = gs
psi0 = lanczos_result.eigenvectors
end
end
# now the operator must be symmetrized!
corr = inner(symmetrize(corr_op, C_N), psi0)
println("Ground state energy: $e0")
println("Ground state correlator: $corr")
println("State vector length = ", size(psi0))
end
function C_N_character(N::Int, k::Int, p::Int)
exp(im * 2 * pi * k * p / N)
end
main()