# 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.

Here we look at the evolution generated by a modified Ising model Hamiltonian,

\begin{equation}
H = - ( \cos(\theta) \sigma_x \otimes \sigma_x + \sin(\theta)/2 (\sigma_z \otimes I + I \otimes \sigma_z ) ) + 1/2 (\sigma_z \otimes \sigma_y + \sigma_y \otimes \sigma_z)
\end{equation}


```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]
    Y = [0. im; -im 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) + 0.5 * ( kron(Z,Y) + kron(Y,Z) )
    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 = 100000  # 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%10000==0 print(n-initial_step+1, ":", Nsteps, " ") end
end
```

    10000:100000 20000:100000 30000:100000 40000:100000 50000:100000 60000:100000 70000:100000 80000:100000 90000:100000 100000:100000 

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_modIsing_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)#, markersize = 10, markerstrokewidth = 2)
#(length(tarray),length(entropy))
```




    

    




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

The entanglement entropy of the initially separable state saturates at $N/2$, as expected.

6c582fae98b3fdc2b14d0b407b6e0833b2b4ff9c
ipynb
Jupyter Notebook
EntanglementEntropy_ModifiedH.ipynb
laurasberna/EntanglementEntropy.jl
e1299890cea105a5cb23227ac796947d9b05b059
Apache-2.0
2021-05-21T12:28:00.000Z
2021-05-21T12:28:00.000Z
EntanglementEntropy_ModifiedH.ipynb
laurasberna/EntanglementEntropy.jl
e1299890cea105a5cb23227ac796947d9b05b059
Apache-2.0
EntanglementEntropy_ModifiedH.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