Differences

This shows you the differences between two versions of the page.

--- 물리:열린_양자계_겨울학교_프로젝트 [2023/01/29 20:34] – [모형] jiwon
+++ 물리:열린_양자계_겨울학교_프로젝트 [2023/09/05 15:46] (current) – external edit 127.0.0.1
@@ Line 55: / Line 55: @@
 =====수치해=====
-이론적으로는 위 행렬의 고윳값과 고유벡터를 구하면 계의 모든 특성을 알아낼 수 있지만, 식의 형태가 매우 복잡해지므로 여기서는 수치해를 구하는 것에 집중하도록 한다. 적분 방법으로는 4차 룽게-쿠타 방법을 이용했고 시뮬레이션 코드는 다음과 같다.
+이론적으로는 위 행렬의 고윳값과 고유벡터를 구하면 계의 모든 특성을 알아낼 수 있지만, 식의 형태가 매우 복잡해지므로 여기서는 수치해를 구하는 것에 집중하도록 한다. 적분 방법으로는 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)$$
+를 선택했다.
+=====결과=====
+<code:julia>
+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
+</code>