# Time evolution of the entanglement entroy of a product state

This project is part of the Computational Physics course held at Perimeter, Fall 2018. 


```julia
using LinearAlgebra
using Arpack
using Plots
```

The following functions are used to compute the entanglement entropy, given a state separated in sub-states A and B.


```julia
function compute_rdmA(PsiAB)
    @assert ndims(PsiAB) == 2 # did we input a matrix?
    @assert size(PsiAB,1) < 2^13 # make sure the resulting matrix is not too large!
    rdm = PsiAB*PsiAB'
end

function compute_entropy(dm)
    p = eigvals(dm) # we diagonalize and use the eigenvalues. the more Julia way
    # has a problem when the log is very large
    entropy = 0.0
    #entropy = - p'*log2.(p) # problems when p[n] is of order 0 (possibly negative)
    for n in 1:size(p,1)
        if abs(p[n]) > 1e-12
            entropy = entropy - p[n]*log2(p[n])
        end
    end
    entropy
end
```




    compute_entropy (generic function with 1 method)



## Product state

We use the outer product of 1 spin states to construct a separable state:

\begin{equation}
| \Psi_{\text{sep}} > = |\Psi>_1 \otimes |\Psi>_2 \otimes \cdots \otimes |\Psi>_N
\end{equation}


```julia
N=2*5 # > 2

function normalizeVec(Vec)
    Vec = Vec/sqrt(abs(Vec'*Vec)) 
end

# here we implement the outer product
function separable(N)
    PsiInitial= reshape( (randn(2) + im*randn(2))*(randn(2) + im*randn(2))', 2^2 )
    for i in 3:N
        PsiInitial = reshape( PsiInitial*(randn(2) + im*randn(2))', 2^i )
    end
    return normalizeVec(PsiInitial)
end

PsiInitial = separable(N);
```


```julia
# checking separability
nA = Int(N/2) # number spins/qubits in A
nB = Int(N/2) # number of spins/qubits in B

PsiAB = reshape(PsiInitial, (2^(nA),2^(nB)))

rdmA = compute_rdmA(PsiAB)

isapprox(compute_entropy(rdmA), 0, atol=10^-14)  # checking that entropy is zero to 10^-14 accuracy
```




    true



## Time evolution 

To obtain $|\Psi(t)\rangle \equiv \exp (-i t H) |\Psi(0)\rangle$, we can use
 $|\Psi(t+\delta t)\rangle = |\Psi(t)\rangle - i \delta t H |\Psi(t)\rangle$. We want to keep the accuracy within $10^{-5}$ at each time step.
 
We want the code to be able to deal with large spin chains, so we do not want to record the full Hamiltonian matrix.


```julia
# Ising model
function buildIsing(theta=pi/4) # Ising model with transverse magnetic field h (critical theta=pi/4 by default)
    X = [0. 1; 1 0]
    Z = [1. 0; 0 -1]
    E = diagm(0=>ones(2))
    XX = kron(X,X)
    HXX = XX
    HZ = kron(Z,E) + kron(E,Z)
    H2 = -(cos(theta)*XX + sin(theta)/2*HZ)
    return H2
end

# Energy shift
theta = pi/4 # critical magnetic field
H2 = buildIsing(theta) # Hamiltonian
D,U = eigen(H2)
shiftE = D[end] 
H2 = H2 - shiftE*diagm(0=>ones(4))

# Sparse multiplication
function multiplyHPsi(Psi, H2) 
    L = length(Psi)  # Dimension of the vector space
    N = convert(Int64,log2(L))  # Number of spins
    HPsi = zeros(L)
    for n=1:N                              # This multiplies by the Hamiltonian H
        Psi = reshape(Psi,(4,2^(N-2)))
        HPsi = reshape(HPsi,(4,2^(N-2)))
        HPsi += H2*Psi                 
        Psi = reshape(Psi,(2,2^(N-1)))
        HPsi = reshape(HPsi,(2,2^(N-1)))
        Psi = permutedims(Psi,(2,1))
        HPsi = permutedims(HPsi,(2,1))
    end
    Psi = reshape(Psi,L)
    HPsi = reshape(HPsi, L)
    return HPsi
end
```




    multiplyHPsi (generic function with 1 method)



## Example: N=12


```julia
# Initialization
init_steps = 1
entropy = zeros(1)
normdiff = zeros(1)
initial_step=1
Nsteps = 0;

N = 2 * 6 
Psi = separable(N);
deltat = 0.0001
```




    0.0001



The for loop implements $|\Psi(t+\delta t)\rangle = |\Psi(t)\rangle - i \delta t H |\Psi(t)\rangle$.


```julia
# Loop over the time steps
initial_step += Nsteps
Nsteps = 80000  # change this number if the entanglement entropy has not reached its plateau 
final_step = initial_step + Nsteps-1

for n=initial_step:final_step
    Psi = Psi .- im * deltat .* multiplyHPsi(Psi,H2)
    if n%400==0
    normdiff = [normdiff; abs(1 - sqrt(abs(Psi'*Psi)))]
    #Psi = normalizeVec(Psinew)
    PsiAB = reshape(Psi, (2^(Int(N/2)),2^(Int(N/2))))
    rdmA = compute_rdmA(PsiAB)
    entropynew = compute_entropy(rdmA)
    entropy = [entropy; entropynew]
    end
    if n%4000==0 print(n-initial_step+1, ":", Nsteps, " ") end
end
```

    4000:80000 8000:80000 12000:80000 16000:80000 20000:80000 24000:80000 28000:80000 32000:80000 36000:80000 40000:80000 44000:80000 48000:80000 52000:80000 56000:80000 60000:80000 64000:80000 68000:80000 72000:80000 76000:80000 80000:80000 

We plot the difference from 1 of the normalization, to check the precision of our time integration.


```julia
tarray = collect(range(0, stop = deltat * Nsteps, length = length(normdiff)))
standard = [10^-5 for i in 1:length(tarray)]
plot(tarray, standard, label = "goal", xlabel = "t", ylabel = "normalization loss")
plot!(tarray, normdiff, label = "numerical", markershape = :hexagon, markerstrokealpha = 0,
    markersize = 2)
```




    

    




```julia
savefig("12spins_Ising_norm.png")
```

## Evolution of the entanglement entropy

We can finally plot the evolution of the entanglement entropy of the half chain with time.


```julia
tarray = collect(range(0, stop = deltat * Nsteps, length = length(normdiff)))
asympt = [N/2 for i in 1:length(tarray)]
plot(tarray, asympt, xlabel = "t", ylabel = "entanglement entropy", label = "N/2")
plot!(tarray, entropy, label = "numerical", markershape = :hexagon, markerstrokealpha = 0,
    markersize = 2)
#(length(tarray),length(entropy))
```




    

    



For the Ising model, with similar parameters as the ones used for the modified Hamiltonian, we do not saturate at the same point.


```julia
savefig("12spins_Ising_ee.png")
```

27feaf7e3550fb240a34a5ce44ccfc5f40d70e50
ipynb
Jupyter Notebook
EntanglementEntropy_Ising.ipynb
laurasberna/EntanglementEntropy.jl
e1299890cea105a5cb23227ac796947d9b05b059
Apache-2.0
2021-05-21T12:28:00.000Z
2021-05-21T12:28:00.000Z
EntanglementEntropy_Ising.ipynb
laurasberna/EntanglementEntropy.jl
e1299890cea105a5cb23227ac796947d9b05b059
Apache-2.0
EntanglementEntropy_Ising.ipynb
laurasberna/EntanglementEntropy.jl
e1299890cea105a5cb23227ac796947d9b05b059
Apache-2.0
2022-02-28T14:14:54.000Z
2022-02-28T14:14:54.000Z
Qwen/Qwen-72B
1. YES
2. YES