はじめに
イオントラップは量子コンピュータの方式の一つです．量子性を維持できる時間であるコヒーレンス時間が非常に長く，量子ゲートの精度が非常に高いといった利点を持ちます．量子コンピュータの実装方式の中では歴史も長く，近年もIonQやQuantinuumといったベンチャーも他方式と遜色ない性能のマシンを作っています．
2量子ビットゲートは量子コンピュータでも最も基本的な操作ですが，イオントラップ版の日本語記事がwikipediaすらなかったので勉強がてら書くことにしました．詳細な内容は占部伸二「個別量子系の物理」朝倉書店（2017）がおすすめです．またイオントラップ全般のレビューとしてC. D. Bruzewicz, et al. Appl. Phys. Rev. 6, 021314 (2019)も参考になりました（arXiv版）．
シミュレーションにはQuTiPを使用しました．インストール方などは日本語の記事もたくさんあるので適当にググってみてください．
フォノンと2量子ビットゲート
イオントラップ量子コンピュータではその名の通り交流電場によって捕捉されたイオンを量子ビットとして用います．トラップ中では電場が時間に依存して変化するため，イオンもそれに引っ張られて振動します．しかし，この振動はレーザー冷却という技術によって極限まで抑えられており，理想的には振動基底状態まで冷却されます．このレベルで冷却されると，イオンの振動の量子性を見ることができます．この量子化された振動をフォノンと呼びます．
トラップ中に複数のイオンがある場合を考えます．イオンは電荷を持っているので，クーロン相互作用によって互いに影響を及ぼします．イオンを捕捉している電場を上手く選ぶと，複数のイオンが協調して振動するようになります．このとき複数のイオンが1つの振動（フォノン）を共有しているので，フォノンを媒介してイオン間に2量子ビットゲートをかけることができます．
フォノンを利用した方法にはCiracとZollerによる方法（論文, wikipedia），MølmerとSørensenによる方法（論文1, 論文2, 論文3, wikipedia），そしてLeibfriedによる方法（論文）があります．本記事では現在の標準である後者2つについて触れます．
状態依存力
フォノンの振動状態を利用してイオン間に量子ゲートをかけるわけですから，イオン（量子ビット）の状態に依存してフォノンの状態が変わらなければなりません．つまり，量子ビットの状態によってイオンに加わる力が変わる状況を実現しなければいけません．このような相互作用は，次のハミルトニアンによって実現します．
H(t)=(f(t)a+f^*(t)a^\dagger)\sigma_\mu
ここで，
f(t)=Ae^{-i\delta t}
であり，$A$はイオンに加える力によって決まる複素数，$\delta$はイオンに加える力の周波数とフォノンの周波数の差であり，離調と呼ばれます．
$\sigma_\mu$はPauli演算子です．Leibfriedゲートでは$\sigma_\mu=\sigma_z$，Mølmer-Sørensenゲートでは
\sigma_\mu=\sigma_{\phi}=\sigma_x \cos\phi + \sigma_y \sin \phi
です．$\phi_a,\phi_v$はイオンを操作するレーザー光の位相によって決定します．
$a^\dagger$と$a$はそれぞれフォノンの生成・消滅演算子です．生成・消滅演算子はフォノンの数状態$\ket{n}$に対して以下のような働きを持ちます．
a\ket{n}=\sqrt{n}\ket{n-1}, a^\dagger\ket{n}=\sqrt{n+1}\ket{n+1}
がさつに言えばフォノンの数を増やしたり減らしたりする演算子です．
驚くべきことに，このハミルトニアンによる時間発展演算子（ユニタリ演算子）は解析的に計算でき，
U(t)=e^{i\hat{\Phi}(t)}\exp\{(\alpha(t)a^\dagger-\alpha^*(t)a)\sigma_\mu\}
となります．ここで$\hat\Phi(t)$は，イオンが位相空間上で動いた経路と始点と終点を結んだ直線で囲んだ領域の面積に比例する幾何的位相で，
\hat\Phi(t)=\left|\frac{A}{\delta}\right|^2\sigma_\mu^2(\delta t-\sin \delta t)，
$\alpha(t)$はイオンに力を加えて強制振動が起きた結果生成されるコヒーレント状態の固有値で
\alpha(t)=\frac{A}{\delta}(1-e^{i\delta t})
です．
さて，色々と面倒そうな用語が導入されましたが全部説明するときりがないので，ページの下部にまとめました．位相空間，コヒーレント状態などは最後の捕捉を参照してください．
長いですがあと少しです（書く方も流石に疲れてきました）．こんなに長く書く予定ではなかったのですが......．LeibfriedゲートとMølmer-Sørensenゲートの違いも書こうと思いましたが，あまりにも長いので興味がある方は「はじめに」に書いた参考文献を参照してください．基本的にはレーザーの当て方が違うだけです．細かいことを言うと，その違いに伴った利点欠点の違い1があったりしますが省略します．
状態依存力のシミュレーション
さて，ほとんど道具もそろったので2量子ビットゲートをやりたいところですが，その前に1個のイオンとフォノンの相互作用を見ていきましょう．まずは必要なライブラリをインポートします．
import numpy as np
import qutip as qp
import matplotlib.pyplot as plt
import matplotlib.colors as color
from matplotlib.cm import ScalarMappable
from tqdm.notebook import trange #notebook形式でない場合は削除．trange->rangeに置き換える．
次に状態依存力による時間発展を計算する関数を定義します．
def f(t, args): #上の解説f(t)と同じ
return args['amplitude']*np.exp(-1j*args['detuning']*t)
def f_conj(t, args): #f(t)の複素共役
return args['amplitude']*np.exp(1j*args['detuning']*t)
def SDF(Atom0, phonon0, detuning, int_strn, time, phi, method): #状態依存力を計算する関数
a=qp.destroy(phonon0.dims[0][0])
ad=qp.create(phonon0.dims[0][0])
if method=='L': #Leibfried
H0=qp.tensor(qp.Qobj(np.zeros((2,2))),qp.Qobj(np.zeros((phonon0.dims[0][0],phonon0.dims[0][0]))))
H1=qp.tensor(qp.sigmaz(),a)
H2=qp.tensor(qp.sigmaz(),ad)
if method=='MS': #Molmer-Sorensen
H0=qp.tensor(qp.Qobj(np.zeros((2,2))),qp.Qobj(np.zeros((phonon0.dims[0][0],phonon0.dims[0][0]))))
H1=qp.tensor(qp.sigmax()*np.cos(phi)+qp.sigmay()*np.sin(phi),a)
H2=qp.tensor(qp.sigmax()*np.cos(phi)+qp.sigmay()*np.sin(phi),ad)
H=[H0,[H1,f],[H2,f_conj]]
return qp.sesolve(H=H,psi0=qp.tensor(Atom0,phonon0),tlist=time,args={'amplitude':int_strn,'detuning':detuning})
qutipでは時間依存のハミルトニアンと時間無依存のハミルトニアンを分けてリストにするので，その都合で関数を分けています．
さて一通り準備が整ったので状態依存力によるイオンとフォノンの状態変化を見ていきましょう．以下のコードと出力は量子回路風に書くと以下の図のようになっています．
Nv=50 #考える最大のフォノン数
g=qp.basis(2,0) #基底状態
e=qp.basis(2,1) #励起状態
v=qp.basis(Nv) #フォノンの真空状態
delta=0.5 #離調
A=0.5 #相互作用の大きさ
tlist=np.linspace(0,2*np.pi/delta,100)
res=SDF((e+g)/np.sqrt(2),v,delta,A,time=tlist,phi=0,method='L')
パラメータは適当ですが，そんなに遠からずだと思います（多分......）．状態依存力の効果を見やすくするために$A$を少し大きめにしています．
この計算結果全体を見るのは大変なので，まずはイオンを$\frac{1}{\sqrt{2}}(\ket{g}+\ket{e})$に射影してトレースアウトし，フォノンのウィグナー関数を見てみましょう．カラーバーまわりはこちらを参考にしました．
M=qp.tensor((e+g)*(e+g).dag()/2,qp.identity(Nv)) #測定演算子
phonon_M=[qp.ptrace(M*r*r.dag()*M,1)/(M*r*r.dag()).tr() for r in res.states] #結果の各要素を測定
xvec=np.linspace(-7,7,200) #wigner関数の各リストを指定
W_lst=[] #wigner関数格納用リスト
for i in trange(0,np.size(tlist),11): #いい感じの時間間隔でwignerを取得
W_lst.append(qp.wigner(phonon_M[i],xvec,xvec))
plt.rcParams["font.size"] = 14 #ここから下はプロット用の設定
fig,ax=plt.subplots(2,5,figsize=(29,12))
for j in range(len(W_lst)):
ax[j//5,j%5].set_aspect('equal')
ax[j//5,j%5].pcolor(xvec,xvec,W_lst[j],cmap='bwr',vmin=np.min(W_lst),vmax=np.max(W_lst))
ax[j//5,j%5].text(4,6,f't={tlist[j*11]:.2f}')
norm=color.Normalize(vmin=np.min(W_lst),vmax=np.max(W_lst))
mappable = ScalarMappable(cmap='bwr',norm=norm)
mappable._A=[]
cbar_pos=ax[1,4].get_position()
cbar_ax=fig.add_axes([0.87,cbar_pos.y0,0.02,cbar_pos.height*2.2])
fig.colorbar(mappable,cax=cbar_ax)
plt.subplots_adjust(right=0.86)
plt.subplots_adjust(wspace=0.1)
fig.text(0.5,0.08,'Position',fontsize='x-large')
fig.text(0.1,0.45,'Momentum',rotation=90,fontsize='x-large')
plt.show()
出力結果
横軸が位相空間上での位置で縦軸が運動量です．正確には質量とか色々かけて次元を持たせなければいけませんが，正規化していると思ってください．時間を追うごとに，中心にあったピークが2つに分裂しながら回転し，最終的に元にもどっているのがわかります．分裂した2つの山の間には干渉が見えます．これはシュレディンガーの猫状態2，あるいは単純に猫状態と呼ばれるもので，逆方向に変位したコヒーレント状態の重ね合わせとして以下のように表されます．
\ket{\rm cat}=\frac{1}{\sqrt{N}}(\ket{\alpha(t)}+\ket{-\alpha(t)}).
解析的な結果では$\alpha(t)$が0になるのは$t=2\pi/\delta$のときで，上の数値計算でも確かにそうなっていることがわかります．
次に量子ビットであるイオンの内部状態を見てみましょう．フォノンをトレースアウトして，混合状態になったイオンの密度行列の各要素と非対角成分の実部と虚部を見てみましょう．
ion_M=[qp.ptrace(r,0) for r in res.states]
fig,ax=plt.subplots(2,3,figsize=(12,6))
ax[0,0].plot(tlist,qp.expect(g*g.dag(),ion_M),label=r'$|g\rangle\langle g|$')
ax[0,0].legend()
ax[0,0].set_ylim(0,1)
ax[1,0].plot(tlist,qp.expect(e*e.dag(),ion_M),label=r'$|e\rangle\langle e|$')
ax[1,0].legend()
ax[1,0].set_ylim(0,1)
ax[0,1].plot(tlist,np.real(qp.expect(g*e.dag(),ion_M)),label=r'$|g\rangle\langle e|$, (Re)')
ax[0,1].legend()
ax[0,2].plot(tlist,np.imag(qp.expect(g*e.dag(),ion_M)),label=r'$|g\rangle\langle e|$, (Im)')
ax[0,2].legend()
ax[1,1].plot(tlist,np.real(qp.expect(e*g.dag(),ion_M)),label=r'$|e\rangle\langle g|$, (Re)')
ax[1,1].legend()
ax[1,2].plot(tlist,np.imag(qp.expect(e*g.dag(),ion_M)),label=r'$|e\rangle\langle g|$, (Im)')
ax[1,2].legend()
plt.subplots_adjust(bottom=0.2,left=0.17,right=1, top=1)
plt.subplots_adjust(wspace=0.3)
fig.text(0.54,0.08,'Time',fontsize='x-large')
fig.text(0.1,0.4,'Expectation value',rotation=90,fontsize='x-large')
plt.show()
これらの結果から1つのイオンによる状態依存力は，フォノンが位相空間上で一周する前後でイオンに対して何の変化も起こしていないことがわかります．また，相互作用の過程では非対角成分の実部がほぼ0まで減ってから再び元に戻っており，イオンの内部状態とフォノンの間に量子もつれが生成されていることがわかります．上で幾何的位相をわざわざ導入したにも関わらず何の影響もないのは，グローバル位相としてかかってしまっているためです．
この幾何的位相は複数のイオンを考慮したときに極めて大きな役割を果たします．
2量子ビットゲート
2つのイオンに同じ力を加えて強制振動を起こすと，上で求めたユニタリ変換中の$\sigma_\mu$は以下のように置き換わります．
\sigma_\mu\rightarrow \sigma_{\mu,1}+\sigma_{\mu,2}.
添え字の1と2はそれぞれどのイオンに対する演算子かを示しています．幾何的位相$\hat\Phi(t)$が$\sigma_\mu^2$に比例することを思い出すと，上の置き換えにより
\hat\Phi(t)\propto (\sigma_{\mu,1}+\sigma_{\mu,2})^2=\sigma_{\mu,1}^2+\sigma_{\mu,1}^2+\sigma_{\mu,1}\sigma_{\mu,2}+\sigma_{\mu,2}\sigma_{\mu,1}.
が得られます．これは2つのイオンの間の相互作用に他なりません．
以上の置き換えを考慮した2量子ビットゲートのコードは以下になります．
def SDFGate(Atom0, phonon0, detuning, int_strn, time, phi, method):
a=qp.destroy(phonon0.dims[0][0])
ad=qp.create(phonon0.dims[0][0])
if method=='L':
H0=qp.tensor(qp.tensor(qp.Qobj(np.zeros((2,2))),qp.Qobj(np.zeros((2,2)))),qp.Qobj(np.zeros((phonon0.dims[0][0],phonon0.dims[0][0]))))
sigma=qp.tensor(qp.sigmaz(),qp.identity(2))+qp.tensor(qp.identity(2),qp.sigmaz())
H1=qp.tensor(sigma,a)
H2=qp.tensor(sigma,ad)
if method=='MS':
H0=qp.tensor(qp.tensor(qp.Qobj(np.zeros((2,2))),qp.Qobj(np.zeros((2,2)))),qp.Qobj(np.zeros((phonon0.dims[0][0],phonon0.dims[0][0]))))
sigma=qp.tensor(qp.sigmax()*np.cos(phi)+qp.sigmay()*np.sin(phi),qp.identity(2))+qp.tensor(qp.identity(2),qp.sigmax()*np.cos(phi)+qp.sigmay()*np.sin(phi))
H1=qp.tensor(sigma,a)
H2=qp.tensor(sigma,ad)
H=[H0,[H1,f],[H2,f_conj]]
return qp.sesolve(H=H,psi0=qp.tensor(Atom0,phonon0),tlist=time,args={'amplitude':int_strn,'detuning':detuning})
最大エンタングルメント状態を作るためには2つの状態ベクトルが均等に混ざり合っている状態を作る必要があります．この条件は
\frac{A}{\delta}=\frac{1}{4}
とすれば満たされます．今回は$\delta=0.5$, $A=0.125$とします．その他のパラメータは上と同じにして，初期状態$\ket{g}\ket{g}$のイオンに対するMølmer-Sørensenゲートを見てみましょう．
A=0.125
delta=0.5
res2=SDFGate(qp.tensor(g,g),v,delta,A,time=tlist,phi=0,method='MS')
ion_DM=[qp.ptrace(r,(0,1)) for r in res2.states]
C=np.array([qp.concurrence(i) for i in ion_DM])
plt.plot(tlist,C)
plt.xlabel("Time")
plt.ylabel("Concurrence")
Concurrenceは2量子ビットに対するエンタングルメント測度で，0でseparable，1で最大エンタングルメント状態を示します．詳細な解説は省略します．上のプロットから確かに最大量子もつれ状態が生成されていることがわかります．では具体的にどのような状態かを見てみましょう．
print(ion_DM[-1])
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[ 5.00001777e-01+0.j 0.00000000e+00+0.j 0.00000000e+00+0.j
-1.19996653e-05-0.5j]
[ 0.00000000e+00+0.j 1.11293787e-11+0.j 1.11293787e-11+0.j
0.00000000e+00+0.j ]
[ 0.00000000e+00+0.j 1.11293787e-11+0.j 1.11293787e-11+0.j
0.00000000e+00+0.j ]
[-1.19996653e-05+0.5j 0.00000000e+00+0.j 0.00000000e+00+0.j
4.99998223e-01+0.j ]]
小さい要素を0に近似すると
\rho\approx\frac{1}{2}\left[\matrix{1&0&0&-i\\0&0&0&0\\0&0&0&0\\i&0&0&1}\right]
状態ベクトルで書くと
\ket{\psi}=\frac{1}{\sqrt{2}}(\ket{gg}-i\ket{ee}).
となります．1量子ビットゲートでこの状態はベル状態に置き換わります．
以上がイオントラップ上での2量子ビットゲートです．
最後に，状態依存力を使ったゲートが実験的にロバストである理由を見てみましょう．Mølmer-SørensenゲートやLeibfriedゲートでは，最終的にイオンとフォノンはほぼ切り離されるので，位相空間上での幾何的位相のみが重要になるのは今まで見てきた通りです．これはすなわち，位相空間上での始点にそれなりの自由度があることを示しています．いままでの計算では真空状態を初期状態として計算してきましたが，この条件は必ずしも必要ありません．加えて言えば，必ずしも純粋状態である必要もありません．フォノンが常に加熱させられ続けるイオンにおいて，この特性は非常に魅力的です．
フォノンの初期状態を熱状態としてこの性質を確認してみましょう．SDFGateで用いているsolverをmesolveに置き換え同様の計算を行います．
A=0.125
delta=0.5
th=qp.thermal_dm(Nv,5)
res2=SDFGate(qp.tensor(g,g)*qp.tensor(g,g).dag(),th,delta,A,time=tlist,phi=0,method='MS')
ion_DM=[qp.ptrace(r,(0,1)) for r in res2.states]
C=np.array([qp.concurrence(i) for i in ion_DM])
plt.plot(tlist,C)
plt.xlabel("Time")
plt.ylabel("Concurrence")
このとき最大のconcurrenceは0.99954...と非常に高い値になり，真空状態を初期状態としたときと同様に高い精度で量子もつれが実現できています．
終わりに
本記事ではイオントラップの2量子ビットゲートをpythonでシミュレーションしてみました．イオントラップの2量子ビットゲートは32qubitsのマシンにおいて平均Fidelityが99.8%という桁外れの精度を持っています．一方でゲートスピードは離調（イオンに加える力の周波数とフォノンの周波数の差）で基本的に律速されます．なぜなら，上で見た通り位相空間上を一周する時間が$2\pi/\delta$だからです．そして離調を大きくすることは難しいです．イオンが2つ以上あると，同じ方向で周波数の異なるフォノンが複数存在します．そのため離調を大きくしすぎると望ましくないフォノンと相互作用してしまい，fidelityを下げる要因となります．したがってAriaなどのプロセッサでは2量子ビットゲートの時間は数百μs程度かかってしまいます．なお，研究レベルではレーザーパルスを工夫して1.6μsで99.8%のfidelityを達成しています．
状態依存力は様々な応用が存在します．例えば更に多くのイオンで同時にこのゲートを使うとGHZ状態をロバストに作ることができます．フォノンでcontinuous variableの誤り訂正符号であるGKP stateを作る際にもこの状態依存力が用いられます（本記事で用いた猫状態の生成回路）．
状態依存力を用いたゲートは非常に洗練された方法ですが，今でも多くの研究が為されており，目が離せません．今回レーザーをどのように当てるかなどの詳細な物理的実装には触れませんでしたが，興味がある方は「はじめに」で紹介した個別量子系の物理がおすすめです．
捕捉
位相空間
まず，位相空間とは（今回の場合）位置と運動量からなる空間のことを指します．例えば調和振動子では以下の図（wikipediaより. public domainです）のように位相空間上を動いていると言えます3．フォノンはイオンの運動に対応しているので，このような空間を考えると取り扱いがしやすいです．以下の図は古典的な調和振動子なので，位相空間上での点が一意に定まっています．しかし今回取り扱っているのは量子状態なので，状態が一点には定まりません4．したがって位相空間において量子状態は定まった点ではなく，確率密度関数（のようなもの）として表現されます．
すなわち，イオンに外部から力を加えてフォノンの状態を変えることは，この位相空間上での状態の変化として計算することができます．
これを考えると幾何的位相$\hat\Phi(t)$は以下のようにイメージすることができます．位相空間上での状態は確率分布じゃなかったのかというツッコミが入りそうですが，今回の場合は期待値の軌跡を取っているのでこのように定まります．
コヒーレント状態
コヒーレント状態の最も直感的なイメージは，いわゆる高校物理で出てくる波のイメージです．すなわち単純なばねの振動であり，レーザー光などもこれにあたります．実際には量子的なゆらぎを考慮する必要がありますが，これは位相空間上においては$x$と$p$の分散が同じガウス分布になります．
堅苦しい言い方をすると，コヒーレント状態$\ket{\alpha}$は消滅演算子の固有状態であり以下のように表されます．
\ket{\alpha}=e^{-\frac{|\alpha|^2}{2}}\sum_{n=0}^{\infty}\frac{\alpha^n}{\sqrt{n!}}\ket{n}
a\ket{\alpha}=\alpha\ket{\alpha}.