2023 KIAS-APCTP 통계물리 겨울학교 조재윤 교수님 프로젝트 내용 정리입니다.
두 개의 에너지 준위를 가지는 입자 두개로 이루어진 계를 생각해보자. 이 때, 각 입자의 두 에너지 준위 사이의 에너지 차이가 ℏω이고 두 입자가 J라는 세기로 상호작용하고 있다면 다음과 같은 해밀토니안을 생각할 수 있다. H=ℏω2(σAz+σBz)+ℏJσAxσBx 그리고 계의 입자와 외부 환경과의 상호작용을 다음과 같은 항으로 도입하자. L0=σA−, L1=σBA, γi=γ 여기서 외부 환경은 고려하고자 하는 계에 비해 매우 커서 이러한 상호작용을 도입하더라도 외부 환경은 변하지 않는다고 가정한다.
이를 가지고 린드블라드 방정식을 다음과 같이 쓸 수 있다. ˙ρ(t)=−iℏ[Hf+Hint,ρ(t)]+1∑i=0γi(Liρ(t)L†i−12L†iLiρ(t)−12ρ(t)L†iLi) 여기서 A와 B에 작용하는 연산자가 곱해져 있을 때 이는 실제로 텐서곱을 의미하고, 둘 중 하나에만 작용하는 연산자는 항등연산자와의 텐서곱을 의미한다는 것을 주의해야 한다. 예를 들어, σA−=σA−⊗I,σAxσBx=σAx⊗σBx 이다.
위 린드블라드 방정식을 풀 때는 밀도행렬을 열벡터로 펼치고, 각 항에 걸려있는 해밀토니안과 결어긋남 항들을 이 밀도행렬에 걸리는 연산자로 취급해서 푸는 것이 편리하다. 즉, 변환 XρY=X⊗YT→ρ 를 이용해서 ˙→ρ(t)=M→ρ(t) 의 형태로 만들자. 이 때 행렬 M을 리우빌리안이라고 부른다. 이렇게 하면, 어떤 양자계의 동역학은 위 리우빌리안 행렬의 고윳값과 고유벡터를 구하는 문제가 된다.
위 변환을 이용해서 우리 계의 리우빌리안을 구해보면 M=(−2γ00iJ00000000−iJ0000−3γ2−iωiJ0000000000−iJ000iJ−3γ2−iω00000000000−iJ0iJ00−γ−2iω00000000000−iJ0000−3γ2+iω00iJ−iJ0000000γ0000−γiJ00−iJ00000000000iJ−γ000−iJ0000000γ0iJ00−γ2−iω000−iJ00000000−iJ000−3γ2+iω00iJ000000000−iJ000−γiJ00000γ00000−iJ00iJ−γ000000γ00000−iJiJ00−γ2−iω0000−iJ00000000000−γ+2iω00iJ0−iJ000000γ0000−γ2+iωiJ000−iJ0γ00000000iJ−γ2+iω0000−iJ0γ0000γ0iJ000) 가 된다.
이론적으로는 위 행렬의 고윳값과 고유벡터를 구하면 계의 모든 특성을 알아낼 수 있지만, 식의 형태가 매우 복잡해지므로 여기서는 수치해를 구하는 것에 집중하도록 한다. 적분 방법으로는 4차 룽게-쿠타 방법을 이용했고 초기 조건으로는 다음과 같은 상태를 고려했다. ρ(0)=1−p4I+p|ψ+⟩⟨ψ+| 여기서 |ψ+⟩=1/√2(|eg⟩+|ge⟩)이다. 즉, p=1이면 벨 상태이고, p=0이면 분리가능한 상태이다.
그리고 물리량으로는 얽힘 엔트로피 SE(t)=TrAρA(t)logρ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