Tutorial 1 - End to end

using ITensors
using ITensorsOpenSystems
using Plots
using LinearAlgebra

The density matrix $\rho$ of a dissipative spin system evolves as

\[\frac{d}{dt}\rho(t) = \mathcal{L}\{\rho(t)\}\]

defined by the Liouvillian superoperator

\[\mathcal{L}\{\rho(t)\} = -i[H,\rho(t)] + \sum_{q=1}^J\mathcal{D}_{A_j}\{\rho\},\]

where $H$ is the system Hamiltonian and

\[\mathcal{D}_A\{\rho\} = A\rho A^\dagger - \frac{1}{2}A^\dagger A \rho - \frac{1}{2}\rho A^\dagger A\]

is the Lindblad dissipator for a jump operator $A$. In this tutorial we will consider a system comprising a chain of $N$ spin-1/2's governed by the XXZ Hamiltonian with a transverse field

\[H = \sum_{i=1}^{N-1}\left(\sigma^x_i\sigma^x_{i+1} + \sigma^y_i\sigma^y_{i+1} + \Delta\sigma^z_i\sigma^z_{i+1}\right) + h\sum_{i=1}^N\sigma^x_i\]

and apply $J=3$ jump operators

\[\begin{aligned} A_1 = \sqrt{\gamma_{\rm in}}\sigma^+_1 A_2 = \sqrt{\gamma_{\rm out}}\sigma^-_N A_3 = \sqrt{\gamma_{\rm dp}}\sum_{i=1}^N\sigma^z_i \end{aligned}\]

describing the injection of spin $\uparrow$ excitations at the first spin, their ejection at the last spin, and dephasing of the collective $z$ magnetization. We will consider small system sizes $N$ to allow comparison to exact diagonalization.

begin
	N    = 3    					# Number of spins
	Δ    = 1.5  					# Anisotropy of the spin-spin coupling
	h    = 0.6  					# Transverse field strength
	γin  = 0.4  					# Spin-injection noise rate
	γout = 0.2  					# Spin-ejection noise rate 
	γdp  = 0.8  					# Collective dephasing noise rate
end;

system = siteinds("S=1/2",N); 		# Define system spin-1/2 Hilbert space

fatsys = fatsiteinds(system); 		# Define the vectorized Hilbert space

Define Hamiltonian and jump operators

We first construct the system Hamiltonian OpSum

begin
	H = OpSum()

	# Open boundary chain:
	for i = 1:(N-1)
    	H += h, "X", i 				# Transverse field
    
    	H += 1, "X", i, "X", i+1 	# XX term
    	H += 1, "Y", i, "Y", i+1 	# YY term
    	H += Δ, "Z", i, "Z", i+1 	# ZZ term
	end
	H += h, "X", N  				# Transverse field on last spin
end;

Thenn define Lindblad jump operators

begin
	A = fill(OpSum(), 3)         	# An array of OpSum's for jump operators

	A[1] += sqrt(γin), "Sp", 1    	# Spin injection
	A[2] += sqrt(γout), "Sm", N   	# Spin ejection

	for i = 1:N
    	A[3] += sqrt(γdp), "Z", i 	# Collective dephasing
	end
end

Construct MPOs

Hmpo = MPO(H,system);           	# Build the system Hamiltonian MPO

LH = im * commutator(H,fatsys); 	# Build the commutator MPO for the fat-system

LD = dissipator(A,fatsys);      	# Build the dissipator MPO for the fat-system

L = LH + LD;                    	# Full Liouvillian MPO for the fat-system

LdagL = apply(dagprimes(L),L);   	# Construct Liouvillian "Hamiltonian"

Compute stationary state

Perform standard DMRG sweeps to find the "ground state" of the Liovillian Hamiltonian.

begin
	sweeps = Sweeps(10)				# Initialize sweeps with number 
    maxdim!(sweeps,50,100,200,400,800,800,1000,1000,2000,2000) # Bond dimension      
    cutoff!(sweeps,1E-10)  			# Desired truncation error
	ρ0 = randomMPS(fatsys, "↑"; linkdims=10) 	# Create a random initial state
end

val, ρ = dmrg(LdagL,ρ0,sweeps);		# Find the steady-state 

Now compute observables from the steady state:

magz = expect(ρ,"Sz"); 				# Compute z-magnetization profile

zz = correlation_matrix(ρ,"Sp","Sm"); 	# Compute S+S- correlation matrix 

Exact diagonalization

Compute the exact result for comparison. Start by constructing the full $2^{2N} \times 2^{2N}$ matrix for the Liouvillian Hamiltonian:

if N > 10
	throw("System size too large for exact diagonalization") 	# Don't continue
end

begin
	LH_ed = im * full_commutator(H,fatsys); 
	LD_ed = full_dissipator(A,fatsys); 
	L_ed = LH_ed + LD_ed;
end

Now fully diagonalize to find the steady-state

begin
	F = eigen(L_ed)				# Diagonalize matrix
	E = F.values 				# Extract eigenvalues
	V = F.vectors 				# Extract eigenvectors
	ind = findmin(E) 			# Locate the smallest eigenvalue (should be zero)
	ρ_ed = V[:,ind] 			# Extract the corresponding eigenvector
	ρ_ed = reshape(ρ_ed, (2^N,2^N)); 	# Reshape vector into a matrix
end

Compute observables from the steady-state

magz_ed = expect(ρ_ed,"Sz"); 				# Compute z-magnetization profile

zz_ed = correlation_matrix(ρ_ed,"Sp","Sm"); # Compute S+S- correlation matrix

Compare results

begin
	plot(1:N, magz_ed, label="exact")  
	scatter!(1:N, magz, label="MPO") 
	heatmap(1:N, 1:N, real(zz), c = :heat)
	fontsize = 10
	nrow, ncol = size(zz)
	ann = [(i,j, text(round(real(zz[i,j]), digits=3), fontsize, :white, :center))
            for i in 1:nrow for j in 1:ncol]
			annotate!(ann, linecolor=:white)
	heatmap(1:N, 1:N, real(zz_ed), c = :heat)
	fontsize = 10
	nrow, ncol = size(zz)
	ann = [(i,j, text(round(real(zz_ed[i,j]), digits=3), fontsize, :white, :center))
            for i in 1:nrow for j in 1:ncol]
			annotate!(ann, linecolor=:white)
end