絶縁体中の電子ポーラロンの局在状態計算：ルチルTiO₂におけるDFT+Uの較正

金属酸化物半導体・絶縁体中では、余剰の電子や正孔が周囲の格子を局所的に歪ませ、その歪みによって自らを捕捉した「小ポーラロン (small polaron)」を形成します。キャリアはバンド伝導ではなく、隣接サイトへの熱活性化ホッピングによって輸送されるため、移動度は同じ材料のバンド計算から予想される値より数桁小さくなります。ルチル型 TiO₂ はその代表例で、余剰電子は特定の Ti サイトに局在して Ti³⁺ (3d¹) を形成し、光触媒活性や n 型導電性を支配します。しかし、この局在状態を第一原理計算で得ること自体が容易ではありません。GGA/LDA の非局在化誤差 (自己相互作用誤差) により、素の計算では余剰電子がセル全体に一様に広がった非物理的な解が安定になってしまうためです。本解析では、第一原理計算ソフトウェア Advance/PHASE の DFT+U 機能を用い、ルチル TiO₂ 中の電子ポーラロンの局在解を実際に得るまでの手順と、その鍵となる Hubbard $U$ の較正手続きを示します。あわせて、局在の成否を定量的に判定する指標を提示し、得られた局在解の構造・電子状態・自己束縛エネルギーを検証します。

Keywords: 第一原理計算, DFTシミュレーション, 小ポーラロン, DFT+U, 電子局在, ルチルTiO₂, スピン分極, 荷電欠陥

計算方法：ポーラロン局在解の求め方

自己束縛と DFT+U

伝導帯に加えられた余剰電子は、Ti 3d 準位に入って周囲の酸素を静電反発で押し広げます。この格子緩和が電子をさらに深く束縛し、ギャップ中に Ti³⁺ 由来の占有準位を生じます。この相乗効果 (電子の局在 → 格子の歪み → さらなる局在) が成立して初めて小ポーラロンが得られます。一方、GGA では余剰電子の自己相互作用誤差により電子が非局在化する方向に強いバイアスがかかるため、Ti 3d に Hubbard $U$ 補正を加えて局在解を安定化させる必要があります。

ポーラロンの束縛の強さは、自己束縛エネルギー $E_{\text{ST}}$ で定量化されます。これは、理想構造上に広がった非局在解と、緩和した局在解のエネルギー差として定義されます。

$$E_{\text{ST}} = E_{\text{deloc}} - E_{\text{loc}}$$

ここで $E_{\text{deloc}}$ は理想構造上の非局在解、$E_{\text{loc}}$ は緩和した局在解の全エネルギーです。$E_{\text{ST}} > 0$ であればポーラロンが束縛されていることを意味します。

対称性の破り方と 2 段階手順

理想結晶は完全な対称性を持つため、そのまま計算すると余剰電子はすべての等価な Ti サイトに均等配分されます。したがって、局在させたい Ti (以下 Ti0) の周囲の TiO₆ 八面体を意図的に膨張させた初期構造から出発し、対称性を明示的に破ります。本解析では文献 [1] の報告する Ti³⁺ の緩和パターン (面内 O を +1%、面外 O を +4%) を初期歪みとして与えました。

手順上の要点として、構造緩和の前に「局在しているか」を確認します。電子状態が非局在のまま構造緩和を始めると、原子に働く力が初期歪みを元に戻す方向に働くため、与えた歪みが単調に潰れて計算資源を浪費します。そこで、まず全原子を固定した一点計算で局在の成否を判定し、局在が確認できてから緩和に進む 2 段階の手順を採りました。

局在度の判定指標

局在の判定には、DFT+U が内部で用いる Ti 3d の占有行列 (occmat) を利用します。サイト $i$ の射影磁気モーメント $m_i = n_{i\uparrow} - n_{i\downarrow}$ を全 Ti について求め、その総和に対する比として局在度 $f$ を定義します。

$$f_{\text{top}} = \dfrac{m_{\text{Ti0}}}{\sum_i m_i}, \qquad f_{\text{2nd}} = \dfrac{m_{\text{2nd}}}{\sum_i m_i}$$

ここで比を用いる点が重要です。射影球の内側に入るモーメントの絶対量は射影子半径に依存するため、$m$ の絶対値そのものは物理的な判定基準になりません。実際、文献 [2] は本系について「電荷は +3/+4 という酸化数が示唆するも互いの値が近く、Ti³⁺ が +2.2、Ti⁴⁺ が +2.6」と報告しており、その差はわずか 0.4 電子です。一方、分布の形 (どのサイトに何割が乗っているか) は射影子半径に依存しにくく、文献の報告値と直接比較できます。

計算モデルと計算条件

計算モデルとして、ルチル TiO₂ の慣用単位胞 (正方晶 P4₂/mnm、6 原子) を 2×2×3 に拡張した 72 原子のスーパーセルを用いました (図1a)。ルチルの TiO₆ 八面体は c 軸方向に稜共有して直線鎖を作り、この鎖に沿った Ti–Ti 距離は 2.959 Å と近接しているため、電子ポーラロンの主要な輸送経路となります (図1b)。格子定数は実験値に固定しました。GGA は体積を過大評価し、ポーラロンの相対安定性を歪めることが指摘されているためです [4]。

図1. (a) 実際に計算した 2×2×3 スーパーセル (72 原子)。濃灰がポーラロンを局在させた Ti0 サイトです。(b) Ti0 サイトを含む a–c 断面における Ti 原子の配置を一部拡大して示したもの。橙が Ti0、青線が [001] 方向の稜共有鎖 (Ti–Ti = 2.959 Å) で、これが電子ポーラロンの主要なホッピング経路にあたります。

局在解の構造を構築する際に、余剰電子は、Ti0 サイトに 1 つの電子を付与する形で導入しました (系全体の電荷は −1)。スピン分極は必須であり、母体は非磁性 (d⁰) であるため「スピン密度が立ったサイト＝ポーラロンの位置」が一対一で対応します。余剰電子 1 個に対応して $n_\uparrow - n_\downarrow = 1$ を全 SCF 反復で固定しました。主な計算条件を表1にまとめます。

表1. 計算モデル・計算条件

原子構造	2×2×3スーパーセル / 72 原子
擬ポテンシャル	PAW
交換相関汎関数	GGA (PBE) ＋ Hubbard \(U\) (Ti 3d)
射影子 (DFT+U)	\(l\) = 2、半径 2.30 bohr
カットオフ (波動関数 / 電荷密度)	36 / 324 Rydberg
k点サンプリング	\(U\) スキャン：Γ点、緩和：2×2×2 (Monkhorst–Pack)
スピン	スピン分極、全スピン \(n_\uparrow-n_\downarrow = 1\) を固定
系の電荷	−1 (余剰電子 1 個)
構造最適化	全原子緩和、最大力 \(5\times10^{-5}\) hartree/bohr

計算コストを削減しつつ、定性的な局在転移の傾向を掴むための $U$ スキャンでΓ点のみを設定しました。実際の構造緩和は2x2x2 k点メッシュを設定しました。

計算結果と考察

Hubbard $U$ の較正：局在転移

まず、初期歪みを与えた構造を固定したまま $U_{\text{eff}}$(Ti 3d) を 4, 6, 8, 10 eV と変化させ、局在度を評価しました (図2)。

図2. 局在度の $U_{\text{eff}}$(Ti 3d) 依存性。青＝ポーラロンを局在させた Ti0 サイト、橙＝[001] 方向の隣接 Ti サイト。緑帯は文献 [3] が報告する局在解の範囲です。$U$ = 8 → 10 eV で 22% の不連続な増加が現れ、局在転移を越えたことを示しています。

$U$ = 4 eV では $f_{\text{top}}$ = 28.1% にとどまり、余剰電子は多数の Ti サイトに分散していました。$U$ の増加とともに局在度は上昇しますが、4 → 6 → 8 eV では 9% 程度の緩やかな増加であるのに対し、8 → 10 eV では 22% と不連続に急増し、$f_{\text{top}}$ = 68.5% に達します。この不連続な変化が局在転移に対応します。

$U$ = 10 eV での分布は、Ti0 に 68.5%、[001] 方向の第 1 近接 Ti に各 5.0% であり、文献の報告値 (対象 Ti に約 68%、軸方向隣接に各約 8% [3]) とよく一致します。2 番目に大きいサイトが [001] 鎖の隣接 Ti であるという空間パターンが再現されている点も、単なる数値の一致ではなく物理的に正しい解が得られていることを裏付けます。

この $U$ = 10 eV という値は一見大きく見えますが、バルクルチルに対する先行研究 [2] と整合します。同文献は $U$ を系統的に走査した結果、「交差点エネルギーは $U_{\text{eff}}$ ≈ 10 eV 付近で飽和し、これは (部分的な非局在ではなく) 完全な局在に達したことを示す」と述べ、以降の全計算を $U_{\text{eff}}$ = 10 eV で行っています。

局在解の構造：Ti³⁺ 小ポーラロン

$U$ = 10 eV で構造緩和を行った結果を図3に示します。36 イオンステップで収束しました (最大力 $3.3\times10^{-5}$ hartree/bohr)。

図3. 緩和後のスーパーセルにおける Ti サイト別の平均 Ti–O 結合長。ポーラロンを局在させた Ti0 (橙、サイト番号 3) のみが理想値から +3.5% 膨張し、他の 23 サイトはすべて理想ルチルの値 (1.959 Å、青破線) に一致しています。

Ti0 の平均 Ti–O 結合長は 2.0269 Å となり、他の 23 サイト (1.9572〜1.9591 Å) に対して +3.5% 膨張しています。結合の内訳は、面外 4 本が 1.948 → 2.042 Å (+4.8%)、面内 2 本が 1.982 → 1.997 Å (+0.8%) であり、ハイブリッド汎関数による文献値 (面外 +4%、面内 +1% [1]) とよく一致します。母体側の 23 サイトが理想値を保っていることは、ポーラロンが単一サイトに収まり、有限サイズ効果が制御できていることを示します。

スピン密度と状態密度

スピン密度分布を図4に示します。局在解では余剰電子が Ti0 の 1 サイトに集中しているのに対し、非局在参照ではすべての Ti にほぼ均等に分散しており、ピーク値には約 15 倍の差があります。母体が d⁰ 非磁性であるルチルでは、スピン密度がそのままポーラロンの位置を可視化します。

図4. スピン密度 $n_\uparrow - n_\downarrow$ の分布 (Ti0 を通る a–c 面のスラブ平均、両パネル共通スケール)。(a) $U$ = 10 eV の局在解では単一サイトに集中。(b) 理想構造上の非局在参照では全 Ti に分散し、ピーク値は約 1/15 です。橙の丸は Ti 原子位置を示します。

図5. 状態密度。上向きが多数スピン、下向きが少数スピンです。(a) 局在解では伝導帯下端から 0.87 eV 下に孤立した占有準位が現れ、少数スピン成分はギャップ領域でゼロです。(b) 非局在参照ではこの準位が消え、2.67 eV のクリーンなギャップとなります。

状態密度 (図5) はさらに明確です。局在解ではギャップ中に孤立した占有準位が現れ、伝導帯下端からの深さは 0.87 eV で、実験・計算の報告値 (伝導帯下端から 0.7〜1.0 eV) と整合します。少数スピン成分が全ギャップ領域でゼロであるため、この準位はTi³⁺ (3d¹) の単一電子に対応する 100% スピン分極した状態であることが確認できます。一方、非局在参照ではこの準位が完全に消失し、2.67 eV のクリーンなギャップとなります。

自己束縛エネルギー

局在解と非局在参照のエネルギー差から、自己束縛エネルギーを評価しました (表2)。

表2. 自己束縛エネルギーの評価 ($U_{\text{eff}}$ = 10 eV)

状態	相対エネルギー \(\Delta E\)
局在解 (緩和構造、余剰電子がTi0 に局在)	0 (基準)
非局在参照 (理想構造)	+0.399 eV (= +0.01465 hartree)
自己束縛エネルギー \(E_{\text{ST}}\)	+0.399 eV

$E_{\text{ST}}$ は正であり、局在解が非局在解より 0.40 eV 安定である、すなわちポーラロンが確かに束縛されていることを示します。この値は、同じ DFT+U の枠組みにおける文献 [2] の再配置エネルギー $\lambda$ ≈ 1.2 eV から期待される範囲 (2 サイト極限の $\lambda/4$ = 0.30 eV とバンド極限の $\lambda/2$ = 0.60 eV の間) に収まります。

ただし、ハイブリッド汎関数 (HSE) による報告値 0.15 eV [1] と比べると大きい点には注意が必要です。これは、バルクでの局在を実現するために大きな $U$ を課したことによる過剰束縛であり、DFT+U とハイブリッド汎関数の系統差として理解されます。

まとめ

第一原理計算ソフトウェア Advance/PHASE の DFT+U 機能を用い、ルチル TiO₂ 中の電子ポーラロンの局在状態を求めました。Hubbard $U$ の系統的な走査により $U$ = 8 → 10 eV に局在転移が存在することを見出し、$U_{\text{eff}}$ = 10 eV で局在度 68.5% (文献値 68〜76%) の Ti³⁺ 小ポーラロンを得ました。得られた局在解は、TiO₆ 八面体の膨張 (面外 +4.8%、面内 +0.8%、文献値 +4%/+1%)、伝導帯下端から 0.87 eV 下の 100% スピン分極したギャップ内準位、正の自己束縛エネルギー +0.399 eV という、小ポーラロンの特徴を再現しています。あわせて、局在度を射影占有行列から定量評価する指標と、構造緩和の前に局在を確認する 2 段階手順を提示しました。本解析で構築した局在解は、Marcus / Holstein 理論に基づくホッピング解析の始状態となります。さらに局在させた終状態を計算すれば、両者を結ぶ経路上のエネルギー曲線から再配置エネルギー、電子カップリング、活性化エネルギーを評価することで、ホッピング速度と移動度が算出可能です。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献

1. A. Janotti, C. Franchini, J. B. Varley, G. Kresse, and C. G. Van de Walle, "Dual behavior of excess electrons in rutile TiO₂", Phys. Status Solidi RRL 7, 199 (2013).
2. N. A. Deskins and M. Dupuis, "Electron transport via polaron hopping in bulk TiO₂: A density functional theory characterization", Phys. Rev. B 75, 195212 (2007).
3. C. S. Ahart, D. Li, J. Blumberger, and S. Liu, "Polaron Transport in TiO₂ from Machine Learning Molecular Dynamics", arXiv:2606.01763 (2026).
4. M. Reticcioli, U. Diebold, and C. Franchini, "Modeling polarons in density functional theory: lessons learned from TiO₂", J. Phys.: Condens. Matter 34, 204006 (2022).

関連ページ

• 第一原理計算ソフトウェア Advance/PHASE
• 解析分野：ナノ・バイオ
• 産業分野：材料・化学