【原】Triborg | 简单两能级系统动力学

Naz摘星星 2022-04-24

展开全文

作者介绍：

加州理工学院 · 访问教授

纽约大学 · 博士后

简单两能级系统动力学

已知无相互作用Hamiltonian为：

而包含相互作用的Hamiltonian为：

相互作用绘景下的系统

从Schrödinger绘景转换成相互作用绘景：

相应的Schrödinger方程也发生变化：

把上边的

写成指数形式，

并考虑到

丢掉两频率加和项，也就是运用旋转波近似（RWA），可得到：

可见，相互作用绘景下，系统的Hamiltonian仍然显含时间，研究不便。

适当选取绘景后的系统

我们可以通过适当选取绘景，把系统Hamiltonian中的时间去掉。引入一个算符：

用它生成一个时间演化算符：

则可求得对应的绘景变换后的Hamiltonian为：

可见不含时间。

3. 求解运动方程—— ——含时Schrödinger方程

变换后的Hamiltonian可求得本征值：

本征向量：

对角化后的Hamiltonian的时间演化算符可变换为三个矩阵乘积：

注意，我们的初态是变换前的，所以需要添上绘景变换的俩矩阵：

这样，给定初态，我们就可以愉快的求解整个系统随时间的变化了：

这就是所谓的谱方法，我们使用一套本征值，五个矩阵都是不变的，一次性算好即可。要传播多久就做多少次矩阵乘法。

此问题还可以用split operator方法求解。我第一次学习的时候忘了包括俩绘景变换矩阵，做出来结果不对，一着急就用下边的QuTip求解了。这几天才想明白。

例子：

设初始态为激发态

，

，求得系统处于激发态的概率：

随时间变化如下：

求解这个系统随时间演化，用QuTip更为方便（为了跟上边不同，初始态选取为基态，求的是基态概率随时间演化）：

import qutip as qt
# parameters of the Hamiltonian 
Omega = 1.0 #Rabi freq
delta0 = 0 #detuning 
delta1 = Omega
#the time−indepdendent hamiltonians
H0 = (delta0/2)*qt.sigmaz() - (Omega/2)*qt.sigmax() 
H1 = (delta1/2)*qt.sigmaz() - (Omega/2)*qt.sigmax()

t0 = 0
t1 = 20 
N=1.0E4
tlist = np.linspace(t0, t1, N)

#expectation values to compute over evolution
proj1 = qt.basis(2,0)*qt.basis(2,0).dag()
e_ops = [ proj1]

## Start from ground state:
psi0 = qt.basis(2,0)

result0 = qt.mesolve(H0, psi0 , tlist , c_ops=[], e_ops=e_ops)
result1 = qt.mesolve(H1, psi0 , tlist , c_ops=[], e_ops=e_ops)

plt.plot(tlist, result0.expect[0],'-',lw=3)
plt.plot( tlist , result1.expect[0],'--',lw=3)

plt.xlabel(r'Time ($1/\Omega$)')
plt.ylabel(r'$P(1)$')
plt.legend([ r'$ \delta = 0$', 'delta = Omega'])
plt.title('Ground−state probability in RWA')
plt.savefig('RWA_qutip.png',format='png',dpi=300)
plt.show()

方法二

（前一个处理方法中）使用interaction picture，再加上旋波近似(RWA)，会得到一个新的含时哈密顿和相应的TDSE：

且对于相互作用绘景来说，有：

所以我们仍然可以暴力求解这个TDSE。且比上一节简单，只需要计算三个矩阵乘积即可：

但是需要注意，此时两个用于对角化的Unitary变换矩阵是含时的。态随时间演化就变成：

考虑初条件：

并且考虑一个有趣的含时

：

直接积分，可以得到所谓的(Di-)Adiabatic passages：

跟使用QuTip计算，得到的结果一致。

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import qutip as qt

Omega = 1.0
H0 = 0*qt.identity(2)

def delta_linear(t, d0, t0):
    return d0*(t/t0)

def sigmaPlusCoeff(t, d0, t0):
    return -(Omega/2)*np.exp(1j*delta_linear(t, d0, t0)*t)
def sigmaMinusCoeff(t,d0, t0):
    return np.conj(sigmaPlusCoeff(t, d0, t0))
sPlus = qt.basis(2,0)*qt.basis(2,1).dag()
sMinus= qt.basis(2,1)*qt.basis(2,0).dag()

d0 = 100*Omega
t0 = 10*Omega

spFast = lambda t, args: sigmaPlusCoeff(t, d0, t0)
smFast = lambda t, args: np.conj(spFast(t, args))

hlist = [H0, [sPlus, spFast],[sMinus, smFast]]

tstart = -t0
tstop  = t0
N = 1.0E4
tlist = np.linspace(tstart, tstop, N)

proj1 = qt.basis(2,0)*qt.basis(2,0).dag()

e_ops = [proj1]

psi0 = qt.basis(2,0)

result = qt.mesolve(hlist, psi0, tlist, c_ops=[], e_ops = e_ops)
plt.plot(tlist, result.expect[0],'-',lw=3)

plt.show()