2023 KIAS-APCTP 통계물리 겨울학교 조재윤 교수님 프로젝트 내용 정리입니다.

두 개의 에너지 준위를 가지는 입자 두개로 이루어진 계를 생각해보자. 이 때, 각 입자의 두 에너지 준위 사이의 에너지 차이가 $\hbar \omega$이고 두 입자가 $J$라는 세기로 상호작용하고 있다면 다음과 같은 해밀토니안을 생각할 수 있다. $$H = \frac{\hbar\omega}2\left(\sigma_z^A+\sigma_z^B\right)+\hbar J\sigma_x^A\sigma_x^B$$ 그리고 계의 입자와 외부 환경과의 상호작용을 다음과 같은 항으로 도입하자. $$L_0 = \sigma_-^A,~~L_1=\sigma_A^B,~~\gamma_i = \gamma$$ 여기서 외부 환경은 고려하고자 하는 계에 비해 매우 커서 이러한 상호작용을 도입하더라도 외부 환경은 변하지 않는다고 가정한다.

이를 가지고 린드블라드 방정식을 다음과 같이 쓸 수 있다. $$\dot\rho(t) = -\frac i\hbar[H_f+H_{\text{int}},\rho(t)]+\sum_{i=0}^1\gamma_i\left(L_i\rho(t)L_i^\dagger-\frac12L_i^\dagger L_i\rho(t)-\frac12\rho(t)L_i^\dagger L_i\right)$$ 여기서 $A$와 $B$에 작용하는 연산자가 곱해져 있을 때 이는 실제로 텐서곱을 의미하고, 둘 중 하나에만 작용하는 연산자는 항등연산자와의 텐서곱을 의미한다는 것을 주의해야 한다. 예를 들어, $$ \sigma_-^A = \sigma_-^A\otimes I,\quad\sigma_x^A\sigma_x^B=\sigma_x^A\otimes\sigma_x^B $$ 이다.

위 린드블라드 방정식을 풀 때는 밀도행렬을 열벡터로 펼치고, 각 항에 걸려있는 해밀토니안과 결어긋남 항들을 이 밀도행렬에 걸리는 연산자로 취급해서 푸는 것이 편리하다. 즉, 변환 $$X\rho Y=X\otimes Y^T\vec\rho$$ 를 이용해서 $$\dot{\vec \rho}(t) = M\vec\rho(t)$$ 의 형태로 만들자. 이 때 행렬 $M$을 리우빌리안이라고 부른다. 이렇게 하면, 어떤 양자계의 동역학은 위 리우빌리안 행렬의 고윳값과 고유벡터를 구하는 문제가 된다.

위 변환을 이용해서 우리 계의 리우빌리안을 구해보면 $$M = \left( \begin{array}{cccccccccccccccc} -2 \gamma & 0 & 0 & i J & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i J & 0 & 0 & 0 \\ 0 & -\frac{3 \gamma }{2}-i \omega & i J & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i J & 0 & 0 \\ 0 & i J & -\frac{3 \gamma }{2}-i \omega & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i J & 0 \\ i J & 0 & 0 & -\gamma -2 i \omega & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i J \\ 0 & 0 & 0 & 0 & -\frac{3 \gamma }{2}+i \omega & 0 & 0 & i J & -i J & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \gamma & 0 & 0 & 0 & 0 & -\gamma & i J & 0 & 0 & -i J & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & i J & -\gamma & 0 & 0 & 0 & -i J & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \gamma & 0 & i J & 0 & 0 & -\frac{\gamma }{2}-i \omega & 0 & 0 & 0 & -i J & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -i J & 0 & 0 & 0 & -\frac{3 \gamma }{2}+i \omega & 0 & 0 & i J & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -i J & 0 & 0 & 0 & -\gamma & i J & 0 & 0 & 0 & 0 & 0 \\ \gamma & 0 & 0 & 0 & 0 & 0 & -i J & 0 & 0 & i J & -\gamma & 0 & 0 & 0 & 0 & 0 \\ 0 & \gamma & 0 & 0 & 0 & 0 & 0 & -i J & i J & 0 & 0 & -\frac{\gamma }{2}-i \omega & 0 & 0 & 0 & 0 \\ -i J & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\gamma +2 i \omega & 0 & 0 & i J \\ 0 & -i J & 0 & 0 & 0 & 0 & 0 & 0 & \gamma & 0 & 0 & 0 & 0 & -\frac{\gamma }{2}+i \omega & i J & 0 \\ 0 & 0 & -i J & 0 & \gamma & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & i J & -\frac{\gamma }{2}+i \omega & 0 \\ 0 & 0 & 0 & -i J & 0 & \gamma & 0 & 0 & 0 & 0 & \gamma & 0 & i J & 0 & 0 & 0 \\ \end{array} \right)$$ 가 된다.

이론적으로는 위 행렬의 고윳값과 고유벡터를 구하면 계의 모든 특성을 알아낼 수 있지만, 식의 형태가 매우 복잡해지므로 여기서는 수치해를 구하는 것에 집중하도록 한다. 적분 방법으로는 4차 룽게-쿠타 방법을 이용했고 초기 조건으로는 다음과 같은 상태를 고려했다. $$\rho(0) = \frac{1-p}4I + p\vert\psi^+\rangle\langle\psi^+\vert$$ 여기서 $\vert\psi^+\rangle = 1/\sqrt2\left(\vert eg\rangle + \vert ge\rangle\right)$이다. 즉, $p=1$이면 벨 상태이고, $p=0$이면 분리가능한 상태이다.

그리고 물리량으로는 얽힘 엔트로피 $$S_E(t) = \text{Tr}_A \rho_A(t)\log\rho_A(t)$$ 를 선택했다.

using LinearAlgebra
using Printf
 
if (length(ARGS)<5)
	println("Usage: julia lindblad.jl J gamma p tf dt")
	exit()
end
 
J = parse(Float32,ARGS[1])
gamma = parse(Float32,ARGS[2])
p=parse(Float32,ARGS[3])
tf = parse(Float32,ARGS[4]) 
dt = parse(Float32,ARGS[5])
omega = 1
 
sx = [0 1; 1 0]
sz = [1 0; 0 -1]
sm = [0 0; 1 0]
id2 = [1 0; 0 1]
id4 = kron(id2,id2)
 
function Hf(omega)
  hf = omega*(kron(sz,id2)+kron(id2,sz))/2.0
  return -im*(kron(hf,id4)-kron(id4,hf))
end
 
function Hint(J)
  hint = J*kron(sx,sx)
  return -im*(kron(hint,id4)-kron(id4,hint))
end
 
function Decoherence(gamma)
  L0 = kron(sm,id2)
  L1 = kron(id2,sm)
 
  decoherence0 = kron(L0,L0)-0.5*(kron(transpose(L0)*L0,id4)+kron(id4,transpose(L0)*L0))
  decoherence1 = kron(L1,L1)-0.5*(kron(transpose(L1)*L1,id4)+kron(id4,transpose(L1)*L1))
  return gamma*(decoherence0+decoherence1)
end
 
function EntanglementEntropy(r)
  r = reshape(r,(4,4))
  c = real((r[1,1]+r[2,2])*(r[3,3]+r[4,4])-(r[3,1]+r[4,2])*(r[1,3]+r[2,4]))
  l1 = 0.5*(1.0+sqrt(1.0-4.0*c))
  l2 = 0.5*(1.0-sqrt(1.0-4.0*c))
  return -l1*log2(l1)-l2*log2(l2)
end
 
 
M = Hf(omega)+Hint(J)+Decoherence(gamma)
dmatrix = [(1-p)/4; 0; (1-p)/4; 0; 0; (1+p)/4; p/2; (1-p)/4; (1-p)/4; p/2; (1+p)/4; 0; 0; (1-p)/4; 0; (1-p)/4]
 
t=0.0
while(t<tf)
  global dmatrix,t,dt
  k1 = dt*M*dmatrix
  k2 = dt*M*(dmatrix+k1/2.0)
  k3 = dt*M*(dmatrix+k2/2.0)
  k4 = dt*M*(dmatrix+k3)
  dmatrix = dmatrix+(k1+2.0*k2+2.0*k3+k4)/6.0
  @printf "%.4f %.4f\n" t EntanglementEntropy(dmatrix)	
  t+=dt
end

STEADY-STATE ENTANGLEMENT IN BIPARTITE SYSTEMS UNDER THE INFLUENCE OF NON-MARKOVIAN ENVIRONMENTS, Ph.D thesis of Joachim Fischbach