RKKY起源スピングラスの数値モデル

RKKY相互作用を起源とするスピングラスは、長距離・振動型の交換結合と乱れが結び付いた多体問題であり、その理解には数値シミュレーションが重要な役割を果たしている。ここでは、RKKYスピングラスを扱う数値計算の枠組みを、モデル化・アルゴリズム・物理量・実験との対応という観点から体系的に整理する。

参考ドキュメント

1. Hukushima, K., Nemoto, K., Exchange Monte Carlo Method and Application to Spin Glass Simulations, Journal of the Physical Society of Japan 65 (6), 1604-1608 (1996). https://www.jstage.jst.go.jp/article/jpsj/65/6/65_6_1604/_article
2. Walstedt, R. E., Walker, L. R., Numerical simulation of spin-glass transition phenomena, Journal of Applied Physics 53, 7985 (1982). https://pubs.aip.org/aip/jap/article-pdf/53/11/7985/18399433/7985_1_online.pdf
3. 大関 真之, スピングラス模型の臨界点と双対変換1 – 厳密解を求めて – (2010). https://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/bk2010.pdf

1. 位置づけ

RKKYスピングラスとは、金属中の局在磁気モーメント同士が伝導電子を介して結合し、その結合が長距離かつ符号振動（強磁性・反強磁性の混在）を示すことで生じるスピングラス状態である。希薄合金 CuMn, AuFe などのいわゆるカノニカル・スピングラスは、その代表例である。

数値計算の観点からは、RKKYスピングラスの特徴は次の4点に要約される。

• 局在スピンの位置が乱雑（希薄合金、準結晶、近似結晶など）
• 相互作用 $J_{ij}$ が長距離かつ距離とともに振動し符号が変わる
• スピンがイジング型・XY型・ハイゼンベルク型など連続対称性を持ち得る
• エネルギー地形が多極小構造を持つため緩和時間が非常に長い

これらを忠実に数値的に扱うことは計算コストの面で容易ではないため、実材に近い「RKKY型模型」と、臨界現象や普遍性の理解に焦点を当てた「簡約化模型（Edwards–Anderson模型など）」の両方が用いられている。

2. 有効ハミルトニアンとRKKY結合の生成

2.1 有効ハミルトニアン

古典スピングラス模型では、一般に次のHeisenberg型有効ハミルトニアンが用いられる。

$$H=-\frac{1}{2}\sum _{i\ne j}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,|{\mathbf{S}}_i|=1$$

イジング模型では ${\mathbf{S}}_i\to S_i=\pm 1$ とし、XY模型では2次元ベクトルスピンを用いる。RKKYスピングラスでは $J_{ij}$ の距離依存と符号がRKKY形式に従う。

2.2 連続RKKYから格子モデルへの写像

等方的3次元自由電子模型に基づくRKKY相互作用は、遠距離で

$$J(r)=J_0\frac{\sin (2k_{F}r)-2k_{F}r\cos (2k_{F}r)}{(2k_{F}r{)}^4}\simeq -J_0\frac{\cos (2k_{F}r)}{(2k_{F}r{)}^3}$$

と書ける:contentReference[oaicite:3]{index=3}。ここで $k_F$ はフェルミ波数、$J_0$ は局在スピンと伝導電子の微視的結合から決まる係数である。

数値模型では、格子上のサイト $i,j$ またはランダム配置された位置 ${\mathbf{r}}_i$ に対し、

$$J_{ij}=J_{0}f(2k_F{r}_{ij})\frac{1}{r_{ij}^3}D(r_{ij})$$

のように離散化することが多い。$f$ は振動因子（$\cos$ など）、$D(r)$ は不純物散乱や電子局在を反映した指数減衰

$$D(r)=\exp (-\frac{r}{\lambda })$$

を含む場合もある。希薄合金では局在スピンは母格子の一部サイトにランダムに配置されるため、$J_{ij}$ の強さと符号は距離と配置の両方に依存する。

2.3 Edwards–Anderson模型との関係

解析的・数値的な扱いやすさのため、スピングラスの臨界現象の研究には、最近接結合のみを持つEdwards–Anderson(EA)模型

$$H_{\mathrm{EA}}=-\sum _{⟨ij⟩}{J}_{ij}{S}_i{S}_j$$

が広く用いられる。ここで $J_{ij}$ は平均ゼロ・分散有限の分布（ガウス分布、$\pm J$ 分布など）から独立にサンプリングされる。

RKKYスピングラスとEA模型の関係は次のように整理される。

• RKKY模型：長距離・振動型 $J_{ij}$、実材の電子状態を反映
• EA模型：短距離・無相関乱数 $J_{ij}$、普遍性と臨界指数の抽出に適する

数値解析では、まずEA模型で普遍的性質（臨界指数、下限臨界次元など）を調べ、その後、RKKY特有の要素（長距離性、振動、希薄性、結晶構造）を導入した模型でより具体的な物質挙動を検討するという戦略が多い。

3. RKKYスピングラス関連模型の分類

RKKY相互作用とスピングラスを結びつける数値模型を整理すると、おおまかに以下のようになる。

模型	相互作用 $J_{ij}$	スピン	長距離性	用途
EA Ising/Heisenberg模型	最近接乱数 $J_{ij}$（ガウス, $\pm J$）	Ising, Heisenbergなど	なし（短距離）	臨界現象・普遍性の研究
長距離べき則模型	$J_{ij}\propto {ϵ}_{ij}/r_{ij}^{\alpha }$（${ϵ}_{ij}$は乱数）	Ising, Heisenberg	あり（べき減衰）	長距離スピングラスの臨界挙動、RKKY近似
格子RKKY模型	RKKY形式 $J_{ij}\propto [\cos (2k_F{r}_{ij})]/r_{ij}^3$ など	Ising, Heisenberg	あり（振動＋減衰）	カノニカル合金や近似結晶に近い記述
バンド構造由来 $J_{ij}$ 模型	DFTなどから得た $J_{ij}$	Heisenberg	実材に依存	物質固有の磁気相図の再現

ここで長距離べき則模型は、RKKYの長距離性を保持しつつ、位相の揺らぎや多バンド性などの詳細を平均化して扱う効果的な模型として用いられている。

4. 古典モンテカルロによる平衡状態の計算

4.1 基本更新則

古典スピングラス模型の平衡統計は、ボルツマン分布

$$P(\{{\mathbf{S}}_i\})\propto \exp (-\beta H)$$

に従う。これをサンプリングする標準的手法はメトロポリス法やヒートバス法である。

イジングスピングラスの場合、メトロポリス更新の1ステップは

1. サイト $i$ を一様ランダムに選ぶ
2. $S_i\to -S_i$ の試行スピン反転後のエネルギー差 $\mathrm{\Delta }E$ を計算
3. 受理確率 $p=min[1,\exp (-\beta \mathrm{\Delta }E)]$ に従って反転を受理・棄却

Heisenbergスピングラスでは、$|{\mathbf{S}}_i|=1$ の制約を保つように微小回転やヒートバス更新を用いる。

スピングラスではエネルギー地形が複雑で低温での緩和が極めて遅いため、単純なメトロポリス法のみでは十分な平衡化が困難になる。この困難を緩和するために発展的手法が導入されている。

4.2 レプリカ交換モンテカルロ法（Exchange Monte Carlo）

レプリカ交換モンテカルロ法（レプリカ交換法、並列テンパリングとも呼ばれる）は、福島・根本らによってスピングラスシミュレーションのために提案された手法であり、現在ではRKKYスピングラスを含む多くの乱雑スピン系に用いられている。

基本的な考え方は次の通りである。

• 温度 $\{T_1<T_2<\cdots <T_M\}$ の異なる $M$ 個のレプリカ系を同時に走らせる
• 各温度では通常のモンテカルロ更新を行う
• 一定ステップごとに隣接温度レプリカ $(T_k,T_{k+1})$ 間でスピン配位を入れ替える試行を行う
• 入れ替えは詳細釣り合いを満たす確率$$p_{\mathrm{ex}}=min[1,\exp (({\beta }_k-{\beta }_{k+1})(E_{k+1}-E_k))]$$で受理する

これにより、低温レプリカは高温レプリカを経由して高エネルギー状態を探索できるようになり、広い位相空間を効率よく探索できる。

4.3 スピングラス特有の物理量

スピングラスでは、通常の磁化 $m$ よりも、スピン配位の重なり（オーバーラップ）を指標とする。

• スピングラス秩序パラメータ（Edwards–Anderson order parameter）： レプリカ $a,b$ に対し

  $$q^{ab}=\frac{1}{N}\sum _i{S}_i^{(a)}{S}_i^{(b)},q_{\mathrm{EA}}^2=[⟨(q^{ab}{)}^2⟩{]}_{\mathrm{av}}$$

  ここで $⟨\cdots ⟩$ は熱平均、$[\cdots {]}_{\mathrm{av}}$ は結合乱れのサンプル平均である。

• スピングラス感受率：

  $${\chi }_{\mathrm{SG}}=N[⟨(q^{ab}{)}^2⟩{]}_{\mathrm{av}}$$

• Binder比、相関長： Binder比 $g$ やスピングラス相関長 ${\xi }_{\mathrm{SG}}(L)$ の有限サイズスケーリングから、転移温度 $T_g$ と臨界指数 $\nu ,\eta$ を決定する。

RKKYスピングラスでは、これらの量に加えて、特定波数の磁気構造因子 $S(\mathbf{q})$ を評価し、螺旋秩序やスピン密度波との競合を調べることも重要である。

5. 長距離RKKY相互作用の数値的扱い

RKKY相互作用は長距離かつ振動型であるため、相互作用の数値的扱いが計算効率を左右する。

5.1 全結合和とその計算量

スピン数を $N$ とすると、単純には全てのペア $(i,j)$ について

$$E=-\frac{1}{2}\sum _{i\ne j}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

を評価する必要があり、計算量は $O(N^2)$ となる。希薄系では有効スピン数が小さいためこの負担は軽くなるが、密な格子模型では大きな負担となる。

5.2 カットオフと有効レンジ

振動しつつ $1/r^3$ で減衰する3次元RKKYでは、ある距離 $r_c$ で相互作用を打ち切る近似がしばしば用いられる。$r_c$ の選び方は、系サイズ $L$ と物理的長さスケール（スピングラス相関長、平均自由行程など）とのバランスで決めることになる。

• 小さ過ぎる $r_c$：RKKYの長距離性が失われ、物理がEA模型に近づき過ぎる
• 大き過ぎる $r_c$：計算量が急増する

したがって、具体的な研究では、$r_c$ を変化させて有限サイズ依存性と合わせて収束性を検証することが多い。

5.3 Ewald和・FFTを用いた高速化

系が周期境界条件を持ち、$J_{ij}$ が平行移動に対して不変な場合、エネルギー評価を畳み込みとして書き換え、FFTを用いて $O(N\log N)$ に高速化することが可能である。RKKYのような振動＋減衰ポテンシャルに対しても、Ewald和に基づく分割（短距離部分は実空間で、長距離部分は波数空間で扱う）が検討されている。

しかし、希薄合金や準結晶のような不規則構造では平行移動対称性が弱く、FFTによる単純な高速化が難しい。この場合は、スパースな隣接リストやGPUによる並列化を用いて $O(N^2)$ 計算を現実的な時間に抑える方法がとられている。

6. RKKYスピングラス数値研究の具体例

6.1 カノニカル合金（CuMn, AuFeなど）の模型化

カノニカル合金では、母格子（fccなど）のサイト群の一部にランダムに磁性イオン（MnやFe）が配置され、これらがRKKY相互作用で結合するとみなされる。数値シミュレーションでは、以下のような戦略がある。

• 結晶格子を固定し、磁性サイトをランダムに占有
• 距離 $r_{ij}$ に基づくRKKY型 $J_{ij}$ を割り当てる
• レプリカ交換モンテカルロ法で平衡状態をサンプリングし、$T_g$ や臨界指数を有限サイズスケーリングから決定する

この種の研究では、実験から得られた凍結温度や非線形磁化率の振る舞いとの比較を通じて、模型の妥当性や普遍性クラスの確認が行われている。

6.2 Tsai型近似結晶・準結晶におけるRKKY磁性

Tsai型近似結晶 Au–Al–Gd や類縁系では、Tsaiクラスタ頂点にあるGd局在スピンがRKKY相互作用を介して多彩な磁気秩序（スピングラス、強磁性、反強磁性など）を示すことが報告されている。

これに対する数値解析では、Gdイオンのみを有効スピンとして取り出し、

• クラスタ構造を保持したグラフ上にRKKY型 $J_{ij}$ を導入
• フェルミ波数 $k_F$ をパラメータとして変化させ、Weiss温度や低温磁気構造の変化を古典モンテカルロで調べる

という手順が用いられている。準結晶特有の幾何学的フラストレーションとRKKY振動の組み合わせにより、多様な磁気状態が生じる点が重要である。

6.3 長距離べき則スピングラス

1次元や3次元で、距離のべき乗で減衰する長距離スピングラス模型

$$J_{ij}=\frac{{ϵ}_{ij}}{r_{ij}^{\alpha }}$$

（${ϵ}_{ij}$ は平均ゼロ・分散1の乱数）は、RKKYの長距離性を抽象化した有効模型として精力的に研究されている。この模型では、指数 $\alpha$ の値を変化させることで、短距離EA模型から無限次元SK模型への連続的な接続を実現でき、理論的な比較が容易である。

6.4 スピングラス薄膜・次元クロスオーバー

RKKY起源のカノニカルスピングラス薄膜では、膜厚を薄くすると凍結温度が低下し、3次元から2次元への次元クロスオーバーが起こることが実験・理論の両面から検討されている。数値的には、有限厚みの系でスピングラス相関長と凍結温度のスケーリングを調べることで、下限臨界次元や有限サイズ効果を解析する枠組みが構築されつつある。

7. 動的性質と非平衡現象のシミュレーション

RKKYスピングラスは、平衡状態だけでなく、遅い緩和・エイジング・履歴依存性などの動的現象でも特徴的な挙動を示す。数値シミュレーションでは、モンテカルロ時間を疑似時間として扱う「ダイナミカルモンテカルロ」枠組みが用いられる。

• コロナレーション関数：

  $$C(t,t_w)=\frac{1}{N}\sum _i[⟨S_i(t_w)S_i(t_w+t)⟩{]}_{\mathrm{av}}$$

  からエイジングや時間スケールを評価する。

• 線形応答： 微小磁場を印加したときの応答と $C(t,t_w)$ を組み合わせ、フラクチュエーション・ディシペーション関係(FDR)の破れを解析する。

RKKY特有の長距離相互作用は、スピン回転が空間的に広がる速度や、局所的なガラスドメインの成長に影響を与えると考えられており、これを反映した時間発展則の検討が進められている。

8. PythonによるIsing RKKY模型の疑似コード

以下に、単純化した2次元正方格子上のIsing RKKY模型をメトロポリス法で扱う疑似コード例を示す。実践的な大規模計算にはMPIやレプリカ交換法などが必要になるが、ここでは枠組みの雰囲気を把握することを目的とする。

import numpy as np

# 格子サイズと物理パラメータ
L = 16
N = L * L
kF = 1.0             # フェルミ波数（無次元化）
J0 = 1.0             # 結合のスケール
T = 0.5              # 温度
beta = 1.0 / T
steps = 100000       # MCステップ数

# 2次元正方格子上のサイトインデックス
coords = [(x, y) for x in range(L) for y in range(L)]

def distance(i, j):
    xi, yi = coords[i]
    xj, yj = coords[j]
    dx = min(abs(xi - xj), L - abs(xi - xj))  # 周期境界
    dy = min(abs(yi - yj), L - abs(yi - yj))
    r = np.sqrt(dx * dx + dy * dy)
    return r

# RKKY型結合の構築（簡略版：2次元で r^{-2} * cos(2kF r) を仮定）
J = np.zeros((N, N), dtype=float)
for i in range(N):
    for j in range(i + 1, N):
        r = distance(i, j)
        if r < 1e-8:
            continue
        val = np.cos(2.0 * kF * r) / (r ** 2)
        J[i, j] = J[j, i] = J0 * val

# スピンの初期状態（ランダム）
spins = np.random.choice([-1, 1], size=N)

def local_field(i, spins):
    # 有効場 h_i = sum_j J_ij S_j
    return np.dot(J[i, :], spins)

def energy(spins):
    # 全エネルギー E = -1/2 sum_{ij} J_ij S_i S_j
    return -0.5 * np.sum(J * np.outer(spins, spins))

# メトロポリス更新
E = energy(spins)
energies = []
for step in range(steps):
    i = np.random.randint(N)
    h_i = local_field(i, spins)
    dE = 2.0 * spins[i] * h_i   # S_i -> -S_i のエネルギー変化
    if dE <= 0.0 or np.random.rand() < np.exp(-beta * dE):
        spins[i] *= -1
        E += dE

    if step % N == 0:
        energies.append(E / N)

# スピングラス重なり q の例（レプリカ2本で測る）
spins2 = np.copy(spins)
# 2本目を別途熱化したと仮定して重なりを計算
q = np.dot(spins, spins2) / N
print("q =", q)

このような単純なコードでも、結合が長距離かつ符号振動を持つことで、多数の準安定状態が生じ、低温での緩和が著しく遅くなる様子を確認できる。ただし、RKKYスピングラスの物理を真に反映するためには、3次元系やHeisenbergスピン、レプリカ交換法、大規模系サイズなどへの拡張が不可欠である。

9. 大規模並列計算とアルゴリズムの発展

スピングラス系は臨界近傍で相関長が大きくなり、有限サイズ効果を抑えた解析には大きな系サイズが必要となる。このため、専用計算機（Janusなど）やスーパーコンピュータを用いた大規模モンテカルロシミュレーションが行われてきた。

アルゴリズム面では、

• レプリカ交換モンテカルロ法
• 集団アニーリング（population annealing）
• クラスタ更新の工夫（低次元スピングラスへの特殊アルゴリズム）

などが提案されている。RKKYスピングラスのような長距離相互作用系では、これらのアルゴリズムと高速なエネルギー評価（FFT、GPU、ベクトル化など）を組み合わせることが鍵となる。

さらに近年は、RKKYスピングラス薄膜や準結晶系など複雑な幾何構造を持つ系に対しても、大規模レプリカ交換法を適用する取り組みが国内外で進んでいる。

まとめと展望

RKKY相互作用を起源とするスピングラスは、長距離・振動型相互作用と乱れが組み合わさった多体問題であり、その数値的理解には、適切に簡約化された模型（EA模型、長距離べき則模型など）と、より実材に即したRKKY型模型の両方を使い分ける枠組みが有効である。古典モンテカルロ法、とりわけレプリカ交換モンテカルロ法は、これらの模型の平衡・臨界性を解析するための基盤として確立されている。

今後の展望としては、第一原理計算や強相関電子計算から得られる実在物質の交換定数 $J_{ij}$ を RKKYスピングラス模型に直接組み込むことで、実験と一体化した「物質別」数値解析がさらに精緻化されると期待される。また、機械学習を用いた有効模型の同定や、動的現象（エイジング、メモリー効果、温度カオスなど）の長時間スケールシミュレーションも重要なテーマであり、RKKYスピングラスの数値計算枠組みは、今後も計算物性・統計物理の交差点として発展していくと考えられる。

参考文献

• Hukushima, K., Nemoto, K., Exchange Monte Carlo Method and Application to Spin Glass Simulations (preprint, arXiv:cond-mat/9512035). https://arxiv.org/abs/cond-mat/9512035

• Reger, J. D., Bhatt, R. N., Young, A. P., Computer simulation of the Heisenberg spin glass with long-range interactions, Physical Review B 37, 5493 (1988). https://link.aps.org/doi/10.1103/PhysRevB.37.5493

• Okamoto, Y., 拡張アンサンブル法：状態空間の効率的なサンプリング（講演スライド）. https://www2.ccs.tsukuba.ac.jp/workshop/algorithm09/slides/tsukuba-09-okamoto.pdf

• 小和口 昌愛, 大規模並列計算を活用したレプリカ交換モンテカルロ法の最適化：一次相転移を示す液晶モデルへの展開（学位論文要旨, 2024, 慶應義塾大学KOARA）. https://koara.lib.keio.ac.jp/xoonips/modules/xoonips/detail.php?koara_id=KO50002002-20246420-0002&ml_lang=ja

• 岡本 祐幸, 「拡張アンサンブル法」の特集にあたって, 分子シミュレーション研究会会誌 “アンサンブル” 12(2) (2010). https://www.jstage.jst.go.jp/article/mssj/12/2/12_2_2_6/_pdf

• 西森 秀稔, スピングラス理論と情報統計力学【新物理学選書】（書誌情報）. https://www.hanmoto.com/bd/isbn/9784007303654

• 東京大学新聞, 【研究室散歩】@統計力学 福島孝治教授（レプリカ交換法とスピングラスの周辺紹介, 2022）. https://www.todaishimbun.org/prof_hukushima_20220819/