S = 5 #sites
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
from scipy.linalg import logm
%config InlineBackend.figure_format = "svg"
import matplotlib.colors as clr
my_cmap = clr.LinearSegmentedColormap.from_list("",[(0,"#011F5B"),(0.5,"#EEEEEE"),(1,"#990000")])
adj = lambda u : np.conj(u.T)

Set Up Spin Operators

s0 = np.array([[1,0],[0,1]],dtype=complex)
sx = np.array([[0,1],[1,0]],dtype=complex)
sy = np.array([[0,-1j],[1j,0]],dtype=complex)
sz = np.array([[1,0],[0,-1]],dtype=complex)

pauli = {"0" : s0,
         "x" : sx,
         "y" : sy,
         "z" : sz,
         "+" : (sx+1j*sy)/2,
         "-" : (sx-1j*sy)/2}

def spin(polarization,s,S):
    if(s==0):
        temp = pauli[polarization]
    else:
        temp = s0
    for i in range(1,S,1):
        if(i==s):
            temp = np.kron(temp,pauli[polarization])
        else:
            temp = np.kron(temp,s0)
    return temp

Functions to Create States

def vacuum(S):
    vac = np.zeros(2**S,dtype=complex)
    vac[-1] = 1
    return vac
def random_state(S): #random state--not in a particular number sector!
    state = vacuum(S)
    b = bin(2**S+np.random.randint(0,2**S))
    for i in range(1,S+1,1):
        if(int(b[-i])==1):
            state = spin("+",i-1,S)@state
    return state

Set Up Hamiltonian

def H_s(S,dir): #open boundaries
    H = np.zeros((2**S,2**S),dtype=complex)
    for i in range(0,S-1,1):
        H += spin(dir,i,S)@spin(dir,i+1,S)
    return H

HX = H_s(S,"x")
HY = H_s(S,"y")
HZ = H_s(S,"z")

def h_s(S,dir):
    #create a dictionary for the Sz operators
    numbers = {}
    for i in range(0,S,1):
        numbers[str(i)] = spin(dir,i,S)
    return numbers

hX = h_s(S,"x")
hY = h_s(S,"y")
hZ = h_s(S,"z")

def Ham(S,Jx=1,Jy=0,Jz=0,hx=0,hy=0,hz=0):
    #build the Hamiltonian
    H = Jx*HX + Jy*HY + Jz*HZ
    for i in range(0,S,1): #we could do a random transverse field here
        H += hx*hX[str(i)] + hy*hY[str(i)] + hz*hZ[str(i)]
    return H
tmax = 20
steps = 10001 #more steps makes the naive method better
dt = tmax/(steps-1)
dirs = ["x","y","z"]

psi0 = vacuum(S) #or make another state

H = Ham(S,Jx=1,hz=0.5)
J1 = spin("+",0,S) #raising operator at left boundary
J2 = spin("-",S-1,S) #lowering operator at right boundary
G1 = 0.1
G2 = 0.1

Vectorized

Vectorize the density matrix by columns $$ \rho = \begin{pmatrix} \rho_{11}&\rho_{21}\\ \rho_{12}&\rho_{22} \end{pmatrix} \quad \mapsto \quad \vec{\rho} = \begin{pmatrix} \rho_{11}\\ \rho_{12}\\ \rho_{21}\\ \rho_{22} \end{pmatrix} $$

Time evolution is now an eigenvalue operator (not superoperator) equation Lvecρ⃗ = i∂tρ⃗ with

$$ L_\text{vec} = 1\otimes H - H^\top \otimes 1 -i\sum_m \frac{\gamma_m}{2} ((J_m^\dagger J_m)^\top\otimes 1 + 1 \otimes J_m^\dagger J_m - 2(J_m^\dagger)^\top \otimes J_m) $$ where we note the switched order of the dagger term versus the master equation. With U(δt) = e−iδtLvec we have ρ⃗(t+δt) = U(δt)ρ⃗(t) where here δt need not be small.

%%time
one = np.identity(2**S)

L_vec = np.kron(one,H) - np.kron(H.T,one)\
        - 1j*G1/2 * (np.kron(one,adj(J1)@J1) + np.kron((adj(J1)@J1).T,one) - 2*np.kron((adj(J1)).T,J1))\
        - 1j*G2/2 * (np.kron(one,adj(J2)@J2) + np.kron((adj(J2)@J2).T,one) - 2*np.kron((adj(J2)).T,J2))

rho = np.einsum("a,b->ab",adj(psi0),psi0)
rho_vec = np.zeros(2**S * 2**S,dtype=complex)
for c in range(2**S):
    rho_vec[c*2**S:(c+1)*2**S] = rho[:,c]

mags = np.zeros((steps,3))
U = expm(-1j*dt*L_vec)

for t in range(steps):
    rho_vec = U@rho_vec
    for c in range(2**S):
        rho[:,c] = rho_vec[c*2**S:(c+1)*2**S]

    #expectation values
    for s in range(S):
        for d in range(3):
            mags[t,d] += np.real(np.trace(rho@spin(dirs[d],s,S))) / S
mags_vec = mags
CPU times: user 2min 28s, sys: 26.6 s, total: 2min 55s
Wall time: 15.2 s

Checks: Tr[ρ] = 1? ρ = ρ†? ρ ≽ 0?

print(np.trace(rho))
print(np.linalg.norm(rho-adj(rho)))
print(np.linalg.eigvalsh(rho))
(1.0000000000000329-1.1616479830397412e-15j)
3.934527700058645e-14
[0.00090453 0.00165106 0.00197228 0.00224695 0.00289197 0.00360004
 0.00410139 0.00489935 0.00527876 0.00630579 0.00670563 0.00718395
 0.00894287 0.01151006 0.01223989 0.01311299 0.01462125 0.01566423
 0.01665747 0.02143925 0.02668842 0.02859218 0.03040516 0.03632071
 0.03913343 0.04674714 0.05325733 0.06629681 0.08532835 0.09721151
 0.11612475 0.21196447]

z-magnetization

plt.plot(np.linspace(0,tmax,steps),mags[:,2])
plt.ylim(-1.1,1.1)
plt.xlabel("t")
plt.ylabel(r"$\langle S^z\rangle$")
Text(0, 0.5, '$\\langle S^z\\rangle$')

Exact determination of steady state

We can numerically solve Lvecρ⃗ss = 0 for ρ⃗ss.

plt.scatter(np.real(np.linalg.eigvals(L_vec)),np.imag(np.linalg.eigvals(L_vec)),marker=".")
plt.xlabel("Re(eig)")
plt.ylabel("Im(eig)")
Text(0, 0.5, 'Im(eig)')

vals = np.linalg.eig(L_vec)[0]
idx = np.argmin(np.abs(vals))
plt.plot(np.sort(np.abs(vals)),marker=".")
plt.xlim(-1,10)
plt.ylim(-0.01,0.1)
plt.xlabel("eigenvalue index")
plt.ylabel("|eig|")
Text(0, 0.5, '|eig|')

Find and construct ρss

vals = np.linalg.eig(L_vec)[0]
idx = np.argmin(np.abs(vals)) #careful--check for uniqueness
print(vals[idx])
vec_ss = np.linalg.eig(L_vec)[1][:,idx]
rho_ss = np.zeros((2**S,2**S),dtype=complex)
for c in range(2**S):
    rho_ss[:,c] = vec_ss[c*2**S:(c+1)*2**S]
print(np.trace(rho_ss))
rho_ss = rho_ss/np.trace(rho_ss)
print(np.linalg.norm(rho_ss-adj(rho_ss)))
print(np.linalg.eigh(rho_ss)[0])
(3.210023821369188e-16-9.727698635837807e-16j)
(5.60115637700194+1.0076661727254077e-13j)
1.1424715863040574e-14
[0.02355502 0.02355502 0.02602636 0.02602636 0.02602636 0.02602636
 0.02800245 0.02800245 0.02800245 0.02800245 0.02875699 0.02875699
 0.03094041 0.03094041 0.03094041 0.03094041 0.03094041 0.03094041
 0.03094041 0.03094041 0.0332896  0.0332896  0.03418661 0.03418661
 0.03418661 0.03418661 0.03678227 0.03678227 0.03678227 0.03678227
 0.04064139 0.04064139]

Steady state magnetization ⟨Sz⟩ and connected correlation function ⟨(Sz)2⟩ − ⟨Sz⟩2

mz = 0
for i in range(0,S,1):
    mz += np.trace(rho_ss@spin("z",i,S)) / S
print(np.real(mz))
1.1311481265541268e-15

mz2 = 0
for i in range(0,S,1):
    mz2 += np.trace(rho_ss@spin("z",i,S))**2 / S
print(np.real(mz2-mz**2))
2.481374156352871e-06

Naive Implementation

Numerically integrate the master equation in time using the Euler method $$ \rho(t+\delta t) = \rho(t) - i\, \delta t \bigg([H,\rho(t)]-i\sum_m \frac{\gamma_m}{2} (J_m^\dagger J_m\rho + \rho J_m^\dagger J_m - 2J_m\rho J_m^\dagger)\bigg) $$

Making the timestep smaller seems to help.

%%time
rho = np.einsum("a,b->ab",adj(psi0),psi0)
mags = np.zeros((steps,3))
trace = np.zeros(steps)
for t in range(steps):
    rho = rho + -1j*dt * ((H@rho-rho@H)\
                       -1j*G1/2*(adj(J1)@J1@rho + rho@adj(J1)@J1 - 2*J1@rho@adj(J1))\
                       -1j*G2/2*(adj(J2)@J2@rho + rho@adj(J2)@J2 - 2*J2@rho@adj(J2)))
    for s in range(S):
        for d in range(3):
            mags[t,d] += np.real(np.trace(rho@spin(dirs[d],s,S))) / S
CPU times: user 8.99 s, sys: 2.12 s, total: 11.1 s
Wall time: 8.86 s

Checks: Tr[ρ] = 1? ρ = ρ†? ρ ≽ 0?

print(np.trace(rho))
print(np.linalg.norm(rho-adj(rho)))
print(np.linalg.eigvalsh(rho))
(0.9999999999999996+0j)
0.0
[-0.02554943 -0.01387799 -0.0023756  -0.00128043 -0.00126105 -0.00071998
  0.00099114  0.00129039  0.0015456   0.00256821  0.0028678   0.00484358
  0.00543601  0.00681364  0.00999824  0.01104068  0.01243214  0.01567089
  0.01695945  0.02029302  0.02106511  0.03038169  0.03169766  0.03605861
  0.03941125  0.04855301  0.05924942  0.07005057  0.09096301  0.10952534
  0.13837507  0.25698296]

We can make this a little less naive by imposing hermiticity, positivity, and trace preservation at each time step

Fixing the trace ρ → ρ/Tr[ρ] and hemiticity ρ → (ρ+ρ†)/2 is straightforward.

Throwing out the negative eigenvalues can be done using a spectral decomposition ρ = ∑ipi|ui⟩⟨ui| so ρ ≈ ∑i : pi > 0pi|ui⟩⟨ui| where we have that ρ|ui⟩=pi|ui⟩.

%%time
rho = np.einsum("a,b->ab",adj(psi0),psi0)
mags_positive = np.zeros((steps,3))
trace = np.zeros(steps)
for t in range(steps):
    rho = rho + -1j*dt * ((H@rho-rho@H)\
                       -1j*G1/2*(adj(J1)@J1@rho + rho@adj(J1)@J1 - 2*J1@rho@adj(J1))\
                       -1j*G2/2*(adj(J2)@J2@rho + rho@adj(J2)@J2 - 2*J2@rho@adj(J2)))
    rho = (rho + adj(rho))/2
    vals, vecs = np.linalg.eigh(rho)
    temp = np.zeros(np.shape(rho),dtype=complex)
    for i in range(np.shape(vals)[0]):
        if(vals[i] > 0):
            temp += vals[i] * np.outer(vecs[:,i],adj(vecs[:,i]))
    rho = temp
    rho /= np.trace(rho)
    for s in range(S):
        for d in range(3):
            mags_positive[t,d] += np.real(np.trace(rho@spin(dirs[d],s,S))) / S
    
CPU times: user 22.6 s, sys: 97.3 ms, total: 22.7 s
Wall time: 11.4 s

Checks: Tr[ρ] = 1? ρ = ρ†? ρ ≽ 0?

print(np.trace(rho))
print(np.linalg.norm(rho-adj(rho)))
print(np.linalg.eigvalsh(rho))
(1+0j)
0.0
[4.02368706e-18 9.94160060e-18 4.09618525e-04 4.39644257e-04
 6.94372323e-04 1.11037775e-03 1.31602362e-03 1.33881574e-03
 1.80087585e-03 3.38834871e-03 3.54193264e-03 5.51638934e-03
 6.54383536e-03 6.86132610e-03 1.01220575e-02 1.06866405e-02
 1.24388192e-02 1.52986279e-02 1.69047108e-02 1.95855626e-02
 2.03256785e-02 2.94157579e-02 3.11941022e-02 3.49363718e-02
 3.78293450e-02 4.69500689e-02 5.58588975e-02 6.79659109e-02
 8.76700049e-02 1.03419199e-01 1.28444810e-01 2.37991874e-01]

z-magnetization

plt.plot(np.linspace(0,tmax,steps),mags_vec[:,2],label="exact vectorized")
plt.plot(np.linspace(0,tmax,steps),mags[:,2],ls="dotted",label="naive integrate")
plt.plot(np.linspace(0,tmax,steps),mags_positive[:,2],ls="dashed",label="integrate (positive)")
plt.ylim(-1.1,0.4)
plt.legend()
plt.xlabel("t")
plt.ylabel(r"$\langle S^z\rangle$")
Text(0, 0.5, '$\\langle S^z\\rangle$')