量子优化控制解决非协性多能级系统的能级泄漏问题

主要阅读文献:

Motzoi, F., Gambetta, J. M., Rebentrost, P., & Wilhelm, F. K. (2009). Simple pulses for elimination of leakage in weakly nonlinear qubits. Physical review letters, 103(11), 110501.

非谐性

非谐性是指系统偏离谐振子的状态,往往会出现其它频率的振荡出现。而超导量子比特就是使用这种非协性量子比特实现的。

在本片文章的背景中,所使用的物理方法比如 phase qubit,具有弱的非协性能级结构。称$|0\rangle, |1\rangle$张成的空间是 qubit 子空间,其它的子空间为 nonqubit 子空间,那么 nonqubit 子空间的频率与 qubit 子空间的频率有较小的差异,因此成为弱非协性。这种系统的特点是, qubit 能级之间的跃迁频率与 nonqubit 能级的跃迁频率相差了一个$\Delta=\omega_2-\omega_1-\omega_1$,如果驱动脉冲的带宽接近 $\Delta$,即它的时常太短接近$\Delta$,就会产生很大的能级泄漏。

总体思想

对于多能级的系统,一般会把其中两个能级单独选出来作为量子比特。但以前的方法不能有效地避免量子比特泄漏到其它的能级,因此本文提出了一种方法可以有效解决这个问题。

在量子优化控制中,我们的目的是使用外加电磁场,对超导电路中的量子比特进行操作,以实现量子门的功能,其中电磁波的磁场分量与自旋耦合,可以改变自旋的动力学演化。

本篇论文主要提出了使用两个相互独立、正交的控制来实现,而驱动电磁波可能是由高斯波形或是正切波形组成的,这两种波形都比较容易产生,但仍然会有较大的误差,因此作者引入了DRAG技术来进一步降低错误率。使用这种方法,可以确定控制参数 $\varepsilon^x, \varepsilon^y$ 以及 $\delta_1$ 的值。

后面作者还使用了 GRAPE 进行优化,其原理是把DRAG中的结果作为初始值再来进行优。

DRAG

在我的理解看来,这里引入 DRAG 的原因在于,我们可以定义一种随时间变化的酉变换 $V(t)$ ,来得到一种近似的有效哈密顿量 $H^V(t)$:

$$H_{eff}(t)=V(t)H(t)V^{\dagger}(t)+i\hbar\dot{V}(t)V^{\dagger}(t)$$

其中 $V(t)$ 就是含时的绝热过程,可以看作很多个子过程的乘积。这个过程的意义在于将量子比特的子空间和泄漏的子空间解耦合,由于这里的泄漏的子空间是处于高能级态,而量子比特的子空间是低能级态,因此实际上也就是把高能级算符和低能级算符解耦合,以方便我们找到导致泄漏的算符,并控制其参数。

这里实际上是把系统的哈密顿量 $H^R$ 中跃迁至高能级的部分看作是微扰 $H_T$,而其它的部分则是我们需要的 $H_0$,因此有 $H^R = H_0 + H_T$,具体来说就是:

$$H_0=\hbar(\delta_1\Pi_1+\frac{\epsilon^x(t)}{2}\sigma^x_{0,1}+\frac{\varepsilon^y(t)\lambda}{2}\sigma^y_{0,1}+\delta_2\Pi_2)$$

$$H_T=\hbar(\frac{\epsilon^x(t)}{2}\sigma^x_{1,2}+\frac{\varepsilon^y(t)\lambda}{2}\sigma^y_{1,2})$$

为什么是这么区分呢?因为只有升降算子 $\sigma_{1,2}^{x(y)}$ 才会将原本低能级的态转换到第二激发态去,因此我认为将这两项看作是微扰是合理的。

现在来计算有效哈密顿量,这实际上是一种近似并且把哈密顿量表示成了 $H_{qubit}(t) \oplus H_{rest}(t)$ 的形式,这将有助于我们对泄漏的能量进行估算,从而找到合理的参数,根据相关文献,其中的 $V(t)H(t)V^{\dagger}(t)$,其中 $V(t) = e^{S(t)}$:

$$VHV^{\dagger}= H + [S, H] + \frac{1}{2}[S, [S, H]]=H_0 + \frac{1}{2}[S, V] + O(V^5)$$

现在来计算一下:

因此:

上式即为 $H^V/\hbar$,但与论文中的结果并不完全吻合,目前尚未找到具体原因,还有待进一步计算。如下文献

Poudel, A., & Vavilov, M. G. (2010). Effect of an Ohmic environment on an optimally controlled flux-biased phase qubit. Physical Review B, 82(14), 144528.

中也对这个式子进行了计算,结果也不太一样,关于这个问题还需要进一步计算。

如果按照论文中的计算结果来考虑:

这里我们需要消掉虚的惯性项(不知道为什么),从而使得:

$$\varepsilon^y = -\frac{\dot{\varepsilon}^x}{\Delta}$$

同时倒数第二项也需要消掉,因为其为投影算子,不能进行转换,因而不是 NOT 门需要的,因此我们得到:

$$\delta_1 = \frac{(\lambda^2 – 4) \varepsilon_x^2}{4\Delta}$$

对于 $H^V/\hbar$,我们还可以进行进一步的展开,得到高阶项,最终我们可以对 $\varepsilon^x(t), \varepsilon^y(t)$ 进行修正,得到:

这些式子即是我们需要在数值仿真中应用的。

数值模拟

参考资料:https://jebej.github.io/Schrodinger.jl/latest/examples/DRAG.html。这里面的代码是用 Julia 写的不是很友好(或者说是过于友好了)。所以还是自己用python些吧

这里数值模拟使用了 python+qutip。我们首先看看参数设定和一些函数定义:

def gaussian(sigma, x):
    return np.exp(- x**2 / 2 / sigma**2)

def epsilon_x(t, args):
    A = args['A']
    sigma = args['sigma']
    t_g = args['t_g']
    
    B = A * gaussian(sigma, -t_g/2)
    normal = (A - B) / A
    epsilon_pi = (A * gaussian(sigma, t) - B) / normal
    
    return epsilon_pi

# 门执行时间
t_g = 6e-9
sigma = t_g / 2

# 高斯波形的振幅
A = np.pi

# 一个谜之常数
UNK = 1 / sigma / (A - gaussian(sigma, -t_g/2) * A) / 2

这里的迷之常数 UNK 是我推导出来的,其作用在于,由于门的时间非常短,所以如果不乘以这样一个因子,就根本看不到结果。但是似乎不是非常精确,因为即使对于二维的理想情况也有一定的误差:

由于论文中对 NOT门进行的模拟,因此我们先来看看三能级的情况,这个时候就应该考虑泄漏到高能级的情况:

# 初始态
psi0 = basis(3, 0)

#投影算符
PI0 = Qobj([[1, 0, 0],[0, 0, 0],[0, 0, 0]])
PI1 = Qobj([[0, 0, 0],[0, 1, 0],[0, 0, 0]])
PI2 = Qobj([[0, 0, 0],[0, 0, 0],[0, 0, 1]])

# 这里取的是δ1=0,δ2=Δ
H2 = Delta * PI2

# 这是升降算符,包括了σ^x_{01(12)}
H  = operators.create(3) + operators.destroy(3)
H += operators.create(3) * (-1j) + operators.destroy(3) * 1j

# 计算哈密顿演化结果
tlist = np.linspace(-t_g/2, t_g/2, 200)
args = {'A':A, 'sigma':sigma, 't_g':t_g}
exp = [PI0, PI1, PI2] # 需要计算期望的算符
output = sesolve([H2, [H * UNK, epsilon_x(tlist, args)]], psi0, tlist, exp)

# 画图
fig, axes = plt.subplots(1,1)
axes.plot(tlist, output.expect[0], label=r'$Prob(|0\rangle)$')
axes.plot(tlist, output.expect[1], label=r'$Prob(|1\rangle)$')
axes.plot(tlist, output.expect[2], label=r'$Prob(|2\rangle)$')
axes.grid()
axes.legend()

# 输出最终的期望值
print("最终的期望值(可以理解为误差):", output.expect[0][-1])

这里需要解释的在于,哈密顿量的构造,即第20行,其哈姆顿量写出来即为:

$$\Delta\Pi_2 + \varepsilon^x\sigma^x_{0,1} + \varepsilon^x\sigma^x_{1,2}+ \varepsilon^y\sigma^y_{0,1} + \varepsilon^y\sigma^y_{1,2}$$

上方代码的输出为,可以看出仅使用高斯脉冲的错误率如下:

那么我现在现在使用 DRAG 技术来试一试,我们首先需要定义 DRAG 的参数(5阶):

# 根据论文中的公式 (10) 得来
    
def drag_x(t, args):
    Delta = args['Delta']
    lambd = args['lambd']
    epsilon_pi = epsilon_x(t, args)
    
    epsilon_x_ = epsilon_pi + (lambd**2 - 4)*(epsilon_pi**3)/(8 * Delta**2)
    epsilon_x_ -= (13*(lambd**4) - 76*(lambd**2) + 112) * (epsilon_pi**5)/(128 * Delta**4)
    
    return epsilon_x_
    
def drag_y(t, args):
    A = args['A']
    Delta = args['Delta']
    lambd = args['lambd']
    sigma = args['sigma']
    epsilon_pi = epsilon_x(t, args)
    
    epsilon_pi_diff = A * epsilon_x(t, args) * t / (sigma ** 2)
    
    epsilon_y = - epsilon_pi_diff/Delta
    epsilon_y += 33*(lambd**2 - 2) * (epsilon_pi ** 2) * epsilon_pi_diff / (24 * Delta**3)
    
    return epsilon_y
    
    
def drag_delta1(t, args):
    epsilon_pi = epsilon_x(t, args)
    
    drag_delta = ((lambd**2) - 4) * (epsilon_pi**2) / (4 * Delta)
    drag_delta -= ((lambd**4) - 7*(lambd**2) + 12) * (epsilon_pi**4) / (16 * (Delta**3))
    
    return drag_delta

随后进行具体的计算:

#初始态
psi0 = basis(3, 0)

# 这里取的是δ2=Δ
H2 = Delta * PI2

# 这里 δ1 不再等于0
H1 = PI1

# 包括了σ^x_{01(12)}和σ^y_{01(12)}
sigma_jk_x = operators.create(3) + operators.destroy(3)
sigma_jk_y = operators.create(3) * (-1j) + operators.destroy(3) * 1j

# 计算哈密顿演化结果
tlist = np.linspace(-t_g/2, t_g/2, 200)
args = {'A':A, 'sigma':sigma, 't_g':t_g, 'lambd':lambd, 'Delta':Delta}
exp = [PI0, PI1, PI2] # 需要计算期望的算符
output = sesolve([H2, [H1 * UNK , drag_delta1(tlist, args)], [sigma_jk_y * UNK , drag_y(tlist, args)], [sigma_jk_x * UNK , drag_x(tlist, args)]], psi0, tlist, exp)

# 画图
fig, axes = plt.subplots(1,1)
axes.plot(tlist, output.expect[0], label=r'$Prob(|0\rangle)$')
axes.plot(tlist, output.expect[1], label=r'$Prob(|1\rangle)$')
axes.plot(tlist, output.expect[2], label=r'$Prob(|2\rangle)$')
axes.grid()
axes.legend()

# 输出最终的期望值
print("最终的期望值(可以理解为误差):", output.expect[0][-1])

仍然把上面的哈密顿量写出来:

$$\delta_1\Pi_1 + \Delta\Pi_2 + \varepsilon^x\sigma^x_{0,1} + \varepsilon^x\sigma^x_{1,2}+ \varepsilon^y\sigma^y_{0,1} + \varepsilon^y\sigma^y_{1,2}$$

输出的结果如下:

可以看出,错误率确实有明显降低,但是却没有达到论文中的效果,目前还不清楚是什么原因。最终,在高斯情况下的曲线为:

可以说,与论文中的差别还挺大的。原因怀疑是前面的自己定义的 UNK 的值有问题,目前还有待进一步研究。

Leave a Reply

Your email address will not be published. Required fields are marked *