MC-05 ボゾン

ボース・ハバード模型は格子上の相互作用するボゾンを記述します。

H = -t \sum_{\langle i,j \rangle} (a_i^\dagger a_j + \text{h.c.}) + \frac{U}{2} \sum_i n_i(n_i-1) - \mu \sum_i n_i,

ここで $t$ はホッピング振幅、 $U$ はオンサイト斥力、 $\mu$ は化学ポテンシャルです。整数充填かつ $U/t$ が大きいとき、系はモット絶縁体になります。ボゾンは相互作用によって局在し、超流動密度は $\rho_s = 0$ です。 $t/U$ が大きくなると、量子ゆらぎがやがて $\rho_s > 0$ の超流動相への転移を引き起こします。このチュートリアルでは、ALPS の worm QMC コードを用いて、2 次元正方格子上の充填率 $\langle n \rangle = 1$ （ $\mu = U/2 = 0.5$ で設定）におけるこの量子相転移を特定します。

転移をまたぐ超流動密度

まず $4 \times 4$ 格子上でホッピングの値を広い範囲にわたって走査し、超流動密度 $\rho_s$ （ALPS では “Stiffness” と呼ばれます）が転移をまたいでどのように変化するかを観察します。ヒルベルト空間は 1 サイトあたり Nmax=2 個のボゾンで打ち切られます。これは単位充填のモットローブ近傍では良い近似です。

コマンドライン

パラメータファイル parm5a：

LATTICE="square lattice"
L=4
MODEL="boson Hubbard"
NONLOCAL=0
U=1.0
mu=0.5
Nmax=2
T=0.1
SWEEPS=500000
THERMALIZATION=10000
{t=0.01;}
{t=0.02;}
{t=0.03;}
{t=0.04;}
{t=0.05;}
{t=0.06;}
{t=0.07;}
{t=0.08;}
{t=0.09;}
{t=0.1;}

NONLOCAL=0 は非局所測定を無効にし、出力を簡潔に保ちます。

parameter2xml parm5a
worm --Tmin 10 --write-xml parm5a.in.xml

--Tmin 10 はチェックポイント間隔を 10 秒に設定します。--write-xml は XML 出力ファイルを書き出し、pyalps は HDF5 の結果と併せてこれを読み込めます。

Python

スクリプト tutorial5a.py：

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

parms = []
for t in [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1]:
    parms.append(
        {
            'LATTICE'        : "square lattice",
            'L'              : 4,
            'MODEL'          : "boson Hubbard",
            'NONLOCAL'       : 0,
            'U'              : 1.0,
            'mu'             : 0.5,
            'Nmax'           : 2,
            'T'              : 0.1,
            'SWEEPS'         : 500000,
            'THERMALIZATION' : 10000,
            't'              : t
        }
    )

input_file = pyalps.writeInputFiles('parm5a', parms)
pyalps.runApplication('worm', input_file, Tmin=10)

評価と作図

超流動スティフネスを読み込み、ホッピングの関数としてプロットします。

data = pyalps.loadMeasurements(pyalps.getResultFiles(prefix='parm5a'), 'Stiffness')
rhos = pyalps.collectXY(data, x='t', y='Stiffness')

plt.figure()
pyalps.plot.plot(rhos)
plt.xlabel('Hopping $t/U$')
plt.ylabel('Superfluid density $\\rho_s$')
plt.title('Bose-Hubbard model on a $4\\times 4$ lattice')
plt.show()

$\rho_s$ は小さな $t/U$ では小さく（ゼロと矛盾しない値）、大きな $t/U$ では有限の値まで増大し、臨界ホッピング $(t/U)_c \approx 0.060$ の付近でクロスオーバーが見られるはずです。

臨界点の特定

$(t/U)_c$ をより精密に決めるために、有限サイズスケーリングを利用します。量子臨界点では $\rho_s \sim L^{-(d+z-2)}$ となります。ここで $d=2$ は次元、 $z=1$ はこの普遍性クラス（3 次元 XY）の動的指数です。 $d=2$ 、 $z=1$ の場合これは $\rho_s \sim L^{-1}$ を与えるので、組み合わせ $\rho_s L$ は臨界点で無次元となり、異なる $L$ の曲線は $t_c$ で交わります。

予想される臨界点の周辺で細かいホッピングの格子点を取り、3 つのシステムサイズ $L = 4, 6, 8$ をシミュレートします。

コマンドライン

パラメータファイル parm5b：

LATTICE="square lattice"
MODEL="boson Hubbard"
NONLOCAL=0
U=1.0
mu=0.5
Nmax=2
T=0.05
SWEEPS=600000
THERMALIZATION=150000
{L=4; t=0.045;}
{L=4; t=0.05;}
{L=4; t=0.0525;}
{L=4; t=0.055;}
{L=4; t=0.0575;}
{L=4; t=0.06;}
{L=4; t=0.065;}
{L=6; t=0.045;}
{L=6; t=0.05;}
{L=6; t=0.0525;}
{L=6; t=0.055;}
{L=6; t=0.0575;}
{L=6; t=0.06;}
{L=6; t=0.065;}
{L=8; t=0.045;}
{L=8; t=0.05;}
{L=8; t=0.0525;}
{L=8; t=0.055;}
{L=8; t=0.0575;}
{L=8; t=0.06;}
{L=8; t=0.065;}

parm5a と比べて低い温度（ $T = 0.05$ ）と長い実行時間が必要なのは、交点をはっきりと分解するためです。

parameter2xml parm5b
worm --Tmin 10 --write-xml parm5b.in.xml

Python

スクリプト tutorial5b.py：

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

parms = []
for L in [4, 6, 8]:
    for t in [0.045, 0.05, 0.0525, 0.055, 0.0575, 0.06, 0.065]:
        parms.append(
            {
                'LATTICE'        : "square lattice",
                'L'              : L,
                'MODEL'          : "boson Hubbard",
                'NONLOCAL'       : 0,
                'U'              : 1.0,
                'mu'             : 0.5,
                'Nmax'           : 2,
                'T'              : 0.05,
                'SWEEPS'         : 600000,
                'THERMALIZATION' : 150000,
                't'              : t
            }
        )

input_file = pyalps.writeInputFiles('parm5b', parms)
pyalps.runApplication('worm', input_file, Tmin=10)

評価と作図

各システムサイズのスティフネスを読み込み、 $L$ を掛けてプロットします。

data = pyalps.loadMeasurements(pyalps.getResultFiles(prefix='parm5b'), 'Stiffness')
rhos = pyalps.collectXY(data, x='t', y='Stiffness', foreach=['L'])

for s in rhos:
    s.y = s.y * float(s.props['L'])

plt.figure()
pyalps.plot.plot(rhos)
plt.xlabel('Hopping $t/U$')
plt.ylabel('$\\rho_s L$')
plt.legend()
plt.title('Finite-size scaling: Bose-Hubbard model')
plt.show()

$L = 4, 6, 8$ の曲線は $(t/U)_c$ の近くで交わるはずです。単位充填の 2 次元ボース・ハバード模型に対する厳密な結果は $(t/U)_c = 0.05974\ldots$ です。

問題

粗い走査では、 $\rho_s$ はどのホッピングの値から明確にゼロでなくなりますか。これは $\rho_s L$ のプロットにおける交点と一致しますか。
$\rho_s L$ の曲線はどこで交わりますか。あなたの見積もりは厳密値 $(t/U)_c = 0.05974$ と比べてどうでしょうか。
有限温度のシミュレーションは $(t/U)_c$ を系統的に過大評価します。それはなぜでしょうか。 $T$ をさらに下げると交点はどうなるでしょうか。
Nmax を 3 に増やしてみてください。結果はどれくらい変わりますか。どのような充填率や相互作用強度で Nmax=2 は悪い近似になるでしょうか。
（発展）より大きなシステムサイズで有限サイズスケーリング解析を繰り返してください。交点の精度はどのように改善しますか。

MC-04 測定 MC-06 量子 Wang-Landau