# EDKit.jl

Julia package for general many-body exact diagonalization calculation. The package provides a general Hamiltonian constructing routine for specific symmetry sectors. The functionalities can be extended with user-defined bases.

## Installation

Run the following script in the Julia Pkg REPL environment:

pkg> add EDKit

Alternatively, you can install `EDKit`

directly from GitHub using the following script.

pkg> add https://github.com/jayren3996/EDKit.jl

This is useful if you want to access the latest version of `EDKit`

, which includes new features not yet registered in the Julia Pkg system but may also contain bugs.

## Examples

Instead of providing documentation, I have chosen to introduce the functionality of this package through practical calculation examples. You can find a collection of Jupyter notebooks in the examples folder, each showcasing various computations.

Here are a few basic examples:

### XXZ Model with Random Field

Consider the Hamiltonian

```
H = \sum_i\left(\sigma_i^x \sigma^x_{i+1} + \sigma^y_i\sigma^y_{i+1} + h_i \sigma^z_i\sigma^z_{i+1}\right).
```

We choose the system size to be `L=10`

. The Hamiltonian needs 3 pieces of information:

- Local operators represented by matrices;
- Site indices where each local operator acts on;
- Basis, if using the default tensor-product basis, only need to provide the system size.

The following script generates the information we need to generate XXZ Hamiltonian:

L = 10
mats = [
fill(spin("XX"), L);
fill(spin("YY"), L);
[randn() * spin("ZZ") for i=1:L]
]
inds = [
[[i, mod(i, L)+1] for i=1:L];
[[i, mod(i, L)+1] for i=1:L];
[[i, mod(i, L)+1] for i=1:L]
]
H = operator(mats, inds, L)

Then we can use the constructor `operator`

to create Hamiltonian:

julia> H = operator(mats, inds, L)
Operator of size (1024, 1024) with 10 terms.

The constructor returns an Operator object, a linear operator that can act on vector/ matrix. For example, we can act `H`

on the ferromagnetic state:

julia> ψ = zeros(2^L); ψ[1] = 1; H * random_state
1024-element Vector{Float64}:
-1.5539463277491536
5.969061189628827
3.439873269795492
1.6217619009059376
0.6101231697221667
6.663735992405236
⋮
5.517409105968883
0.9498121684380652
-0.0004996659995972763
2.6020967735388734
4.99027405325114
-1.4831032210847952

If we need a matrix representation of the Hamitonian, we can convert `H`

to julia array by:

julia> Array(H)
1024×1024 Matrix{Float64}:
-1.55617 0.0 0.0 0.0 … 0.0 0.0 0.0
0.0 4.18381 2.0 0.0 0.0 0.0 0.0
0.0 2.0 -1.42438 0.0 0.0 0.0 0.0
0.0 0.0 0.0 -1.5901 0.0 0.0 0.0
0.0 0.0 2.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 2.0 … 0.0 0.0 0.0
⋮ ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 2.0 0.0 0.0
0.0 0.0 0.0 0.0 … 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -1.42438 2.0 0.0
0.0 0.0 0.0 0.0 2.0 4.18381 0.0
0.0 0.0 0.0 0.0 0.0 0.0 -1.55617

Or use the function `sparse`

to create the sparse matrix (requires the module `SparseArrays`

being imported):

julia> sparse(H)
1024×1024 SparseMatrixCSC{Float64, Int64} with 6144 stored entries:
⠻⣦⣄⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢹⡻⣮⡳⠄⢠⡀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠙⠎⢿⣷⡀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠲⣄⠈⠻⣦⣄⠙⠀⠀⠀⢦⡀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠈⠳⣄⠙⡻⣮⡳⡄⠀⠀⠙⢦⡀⠀⠀⠀⠈⠳⣄⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠙⠮⢻⣶⡄⠀⠀⠀⠙⢦⡀⠀⠀⠀⠈⠳⣄⠀
⢤⡀⠀⠀⠀⠀⠠⣄⠀⠀⠀⠉⠛⣤⣀⠀⠀⠀⠙⠂⠀⠀⠀⠀⠈⠓
⠀⠙⢦⡀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠘⠿⣧⡲⣄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠙⢦⡀⠀⠀⠀⠈⠳⣄⠀⠀⠘⢮⡻⣮⣄⠙⢦⡀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠈⠳⠀⠀⠀⣄⠙⠻⣦⡀⠙⠦⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠈⢿⣷⡰⣄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠈⠃⠐⢮⡻⣮⣇⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠙⠻⣦

### Solving AKLT Model Using Symmetries

Consider the AKLT model

```
H = \sum_i\left[\vec S_i \cdot \vec S_{i+1} + \frac{1}{3}\left(\vec S_i \cdot \vec S_{i+1}\right)^2\right],
```

with the system size chosen to be `L=8`

. The Hamiltonian operator for this translational-invariant Hamiltonian can be constructed using the `trans_inv_operator`

function:

L = 8
SS = spin((1, "xx"), (1, "yy"), (1, "zz"), D=3)
mat = SS + 1/3 * SS^2
H = trans_inv_operator(mat, 1:2, L)

The second input specifies the indices the operators act on.

Because of the translational symmetry, we can simplify the problem by considering the symmetry. We construct a translational-symmetric basis by:

B = TranslationalBasis(L=8, k=0, base=3)

Here, `L`

is the length of the system, and `k`

labels the momentum `k = 0,...,L-1`

(integer multiplies of 2π/L). The function `TranslationalBasis`

returns a basis object containing 834 states. We can obtain the Hamiltonian in this sector by:

julia> H = trans_inv_operator(mat, 1:2, B)
Operator of size (834, 834) with 8 terms.

In addition, we can take into account the total `S^z`

conservation, by constructing the basis

B = TranslationalBasis(L=8, N=8, k=0, base=3)

where the `N`

is the filling number with respect to the all-spin-down state. N=L means we select those states whose total `Sz`

equals 0 (note that we use 0,1,2 to label the `Sz=1,0,-1`

states). This gives a further reduced Hamiltonian matrix:

julia> H = trans_inv_operator(mat, 1:2, B)
Operator of size (142, 142) with 8 terms.

We can go one step further by considering the spatial reflection symmetry.

B = TranslationParityBasis(L=8, N=0, k=0, p=1, base=3)

where the `p`

argument is the parity `p = ±1`

.

julia> H = trans_inv_operator(mat, 1:2, B)
Operator of size (84, 84) with 8 terms.

### PXP Model and Entanglement Entropy

Consider the PXP model

```
H = \sum_i \left(P^0_{i-1} \sigma^x_i P^0_{i+1}\right).
```

Note that the model is defined on the Hilbert space where there is no local `|↑↑⟩`

configuration. For system size `L=20`

and in sector `k=0,p=+1`

, the Hamiltonian is constructed by:

mat = begin
P = [1 0; 0 0]
kron(P, spin("X"), P)
end
pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[mod(i, length(v))+1]==0 for i=1:length(v))
basis = TranslationParityBasis(L=20, f=pxpf, k=0, p=1)
H = trans_inv_operator(mat, 3, basis)

where `f`

argument is the selection function for the basis state that can be user-defined. We can then diagonalize the Hamiltonian. The bipartite entanglement entropy for each eigenstate can be computed by

vals, vecs = Array(H) |> Hermitian |> eigen
EE = [ent_S(vecs[:,i], 1:L÷2, basis) for i=1:size(basis,1)]