反磁界（長距離相互作用）の数値解法

反磁界（demagnetizing field, stray field）は磁化分布全体に依存するため、マイクロ磁化計算において計算量と精度を最も支配しやすい項である。本稿では、方程式の定式化から、代表的な高速計算法（FFT畳み込み、FEM/BEM、FMMなど）と境界条件の扱いを体系的に整理する。

参考ドキュメント

• A. Vansteenkiste et al., The design and verification of MuMax3, AIP Advances 4, 107133 (2014) https://pubs.aip.org/aip/adv/article-pdf/doi/10.1063/1.4899186/12878560/107133_1_online.pdf
• K. M. Lebecki et al., Periodic boundary conditions for demagnetization interactions in micromagnetic simulations (OOMMF related) https://math.nist.gov/~MDonahue/pubs/pbc-in-oommf_4.6.pdf
• 日本磁気学会スクール資料：マイクロマグネティックス（講義資料） https://www.magnetics.jp/archive/seminar/school/ss_034/ss_034_matsuyama.pdf

1. 反磁界とは何か

1.1 静磁場の基本式

時間変化の遅い（準静的）マイクロマグネティクスでは、磁束密度 $\mathbf{B}$ と磁場 $\mathbf{H}$、磁化 $\mathbf{M}$ の関係は

$$\mathbf{B}={\mu }_0(\mathbf{H}+\mathbf{M})$$

である。真空透磁率は ${\mu }_0=4\pi \times {10}^{-7}\mathrm{H}/\mathrm{m}$ である。

磁化が作る静磁場は、いわゆる反磁界（demagnetizing field, stray field）であり、LLGで用いる有効磁場の一部として

$${\mathbf{H}}_{\mathrm{eff}}={\mathbf{H}}_{\mathrm{ex}}+{\mathbf{H}}_{\mathrm{ani}}+{\mathbf{H}}_{\mathrm{ext}}+{\mathbf{H}}_d+\cdots$$

に加算される。

1.2 反磁界の方程式

磁性体を含む全空間で、静磁場は

$$\mathrm{\nabla }\times {\mathbf{H}}_d=\mathbf{0},\mathrm{\nabla }\cdot ({\mathbf{H}}_d+\mathbf{M})=0$$

を満たす。回転がゼロなのでスカラーポテンシャル $\varphi$ を導入でき、

$${\mathbf{H}}_d=-\mathrm{\nabla }\varphi$$

とおける。これを代入すると（磁性体領域で）

$${\mathrm{\nabla }}^2\varphi =\mathrm{\nabla }\cdot \mathbf{M}$$

という Poisson 形式が得られる。ここが Poisson 法の出発点である。

1.3 反磁界エネルギー

反磁界エネルギー（静磁エネルギー）は、一般に

$$E_d=-\frac{{\mu }_0}{2}\int \mathbf{M}\cdot {\mathbf{H}}_{d}dV$$

で与えられる。多くの場合 ${\mathbf{H}}_d$ は $\mathbf{M}$ と反平行成分を持つため、$\mathbf{M}\cdot {\mathbf{H}}_d<0$ となり、結果として $E_d>0$ になる。

1.4 単位の注意

• $\mathbf{H}$ は A/m である
• $\mathbf{B}$ は T である（$B={\mu }_{0}H$ は真空でのみ成立）
• 実装内で「外場を A/m で与える」「磁化は $\mathbf{M}=M_s\mathbf{m}$（$\mathbf{m}$ は無次元）」を徹底すると混乱が減る

目安として、$H={10}^6\mathrm{A}/\mathrm{m}$ は $B={\mu }_{0}H\approx 1.26\mathrm{T}$ に相当する。

1.5 磁荷

磁化の発散が体積磁荷、表面法線成分が表面磁荷として振る舞うと見なすと、

• 体積磁荷密度 $\rho m=-\nabla ·M$
• 表面磁荷密度 $\sigma m=M·n$

を用いて磁位を積分表示できる。これにより、反磁界が本質的に長距離（$1/r^2$）で結合する理由が明確になる。

2. 離散化の基本

有限差分（直交格子）では、セル $i$ の平均磁化 $M_i$ からセル中心の反磁界 $Hdi$ を求める離散化がよく用いられる。

$${\mathbf{H}}_{d,i}=-\sum _j{\mathbf{N}}_{i-j}{\mathbf{M}}_j$$

ここで $N$ は 3×3 の離散反磁界テンソル（demagnetizing tensor kernel）であり、相対位置 $(i-j)$ のみに依存する（平行移動不変性）。この構造が高速化の鍵となる。

重要な設計点は次である。

• セル内磁化を一様とみなす近似（piecewise constant）か、平均場として扱うか
• セル中心で評価するか、セル平均として評価するか
• 形状境界が階段状（staircase）になる誤差をどう扱うか

3. FFTによる高速畳み込み

3.1 FFT化の狙い

上式をそのまま計算すると、$N$セルに対して $O(N^2)$ となり、大規模系で破綻する。畳み込み定理を使えば

• 反磁界テンソル $N$ のフーリエ変換 $𝓕[N]$
• 磁化 $M$ のフーリエ変換 $𝓕[M]$

から、周波数空間で点ごとの積を取り、逆変換することで

$${\mathbf{H}}_d=-{\mathcal{F}}^{-1}\{\mathcal{F}[\mathbf{N}]\cdot \mathcal{F}[\mathbf{M}]\}$$

を $O(NlogN)$ で評価できる。

3.2 畳み込みカーネル（離散反磁界テンソル）

直交格子上で「一様に磁化した直方体セル同士の相互作用」を厳密積分した式を用いて $N$ を構成する方法が広く使われる。 このとき N は対称性を持ち、自己項（同一セル）と近距離項の扱いが精度を左右する。

運用上の要点は次である。

• $N$ は最初に一度だけ計算し、以後は再利用する
• 形状が固定なら $N$ は固定でよい（外場掃引や時間発展でも不変）
• 複数材料（Msが空間変化）では、$M=Msm$ として $M$ を更新しつつ同じ $N$ を使うのが基本になる

3.3 無限空間（開境界）をFFTで扱う

FFTは本質的に周期畳み込みを計算するため、そのままだと“折り返し相互作用”が混入する。無限空間に近い開境界を近似するには、計算領域外をゼロ磁化で埋めるパディングを行い、循環畳み込みが実質的に線形畳み込みになるようにする。

典型的には

• 各方向でサイズを2倍程度に拡張（$2Nx,2Ny,2Nz$）
• 域外は $M=0$ として埋める
• 逆FFT後に元の領域だけ切り出す

という操作を行う。パディング量が少ないと、境界近傍で反磁界が歪む。

3.4 周期境界条件（PBC）とEwald的扱い

薄膜の2次元周期、バルクの3次元周期など、周期境界条件を課したい場合は、カーネル $N$ 自体を周期和した形に置き換える必要がある。単純なゼロパディングではなく、

• PBCに対応する反磁界テンソル $N_{PBC}$ を構成する
• その上でFFT畳み込みを行う

という流れになる。2D周期（面内周期）では特に、真空層の取り方や実効的な遠方境界の扱いが結果に影響する。

3.5 多層膜・スピントロニクス構造

多層膜を一様格子で全空間離散化すると、非磁性スペーサや真空層のために無駄なセルが増えやすい。層方向を不等間隔に扱いつつ、面内はFFT畳み込みで高速化する「多層畳み込み」は、薄膜多層系で計算効率を改善する発想である。

4. Poisson ポテンシャル法

4.1 目的

• ${\mathrm{\nabla }}^2\varphi =\mathrm{\nabla }\cdot \mathbf{M}$ を数値的に解き、${\mathbf{H}}_d=-\mathrm{\nabla }\varphi$ を得る
• 空気箱（airbox）を明示し、境界条件を明確にした検証系を作る
• 将来的にFEM（曲面形状、磁気弾性との統合）へ自然につなぐ

4.2 離散化の流れ

1. 磁化を無次元 $\mathbf{m}$ から $\mathbf{M}=M_s\mathbf{m}$ に戻す
2. 計算箱（磁性体＋空気）に拡張し、右辺 $\mathrm{\nabla }\cdot \mathbf{M}$ を作る
3. Poisson 方程式を離散化して線形方程式 $A\varphi =b$ を作る
4. 反復解法（CG など）で $\varphi$ を解く
5. その勾配から ${\mathbf{H}}_d=-\mathrm{\nabla }\varphi$ を得る（中心差分など）

4.3 境界条件と airbox の意味

Poisson を有限領域で解く以上、外周境界条件を選ぶ必要がある。簡単な選択肢が Dirichlet 境界（外周で $\varphi =0$）である。

• airbox が小さいと、外周条件が試料の場に直接影響し、反磁界が系統的に歪む
• airbox を大きくすると、試料から外周までの距離が増え、境界条件の影響が弱まる

したがって air_cells を増やして結果が頭打ちになるかを見れば、「境界の影響が十分小さい」ことを確認できる。この sweep はその確認でもある。

4.4 長所と注意点

長所

• 境界条件を明示するため、差の原因が追いやすい
• 空気箱や材料の領域分割を入れやすい
• FEM化・異方的メッシュ・他物理（弾性）との統合に親和的である

注意点

• FFTに比べて重い（3Dで未知数が増えるため）
• 収束性が前処理に依存する（pyamg 等の前処理が効く）
• Dirichlet/Neumann の選び方で系統誤差が出るので、airboxで抑える設計が必要である

5. 2手法の比較指標

あなたの demag_compare では、同一の $\mathbf{m}(x)$ を固定し、FFT と Poisson で ${\mathbf{H}}_d$ と $E_d$ を計算し、以下で比較している。

5.1 ベクトル場の差：相対L2誤差

$$\mathrm{relL}2(\mathbf{H})=\frac{\parallel {\mathbf{H}}_{fft}-{\mathbf{H}}_{poi}{\parallel }_2}{\parallel {\mathbf{H}}_{fft}{\parallel }_2}$$

• 0 に近いほど一致である
• 局所的な強度差にも敏感である

5.2 パターン一致：相関係数

ベクトル成分を並べた相関係数 $\mathrm{corr}(\mathbf{H})$ を取ると、スケール差があっても空間パターンが一致しているかを見やすい。

5.3 方向一致：$\cos \theta$ の平均

各セルで

$$\cos \theta =\frac{{\mathbf{H}}_{fft}\cdot {\mathbf{H}}_{poi}}{|{\mathbf{H}}_{fft}||{\mathbf{H}}_{poi}|}$$

を計算し、その平均を見る。これは「方向のバグ」を強く検出する指標である。$\cos \theta \simeq 1$ なら方向は揃っている。

5.4 エネルギー一致：$E_d$ の相対差

$$E_{\mathrm{rel}}=\frac{|E_{fft}-E_{poi}|}{|E_{poi}|}$$

LLG緩和や磁区形成ではエネルギーが重要であるため、$E_d$ の一致は特に意味が大きい。

5.5 sweep の読み方

• air_cells を増やすと Poisson が外周の影響から解放され、エネルギー差が縮む傾向が出る
• fft_padding を増やすと FFT の開境界近似が改善し、同様に差が縮む傾向が出る
• 方向一致（$\cos \theta$）が先に 1 に近づき、強度差（relL2）が後から詰まることがある
  • 境界条件が違うと、同じパターンでも強度が系統的にずれるためである
  • その場合、airbox/padding を増やしたときに差が頭打ちになる点が「この2手法の差の下限」である

6. FEM/BEMハイブリッド

反磁界問題の難点は外部の無限領域である。境界要素法（BEM）は外部領域を境界積分に落とし込めるため、無限空間に自然に対応できる。

ただしBEMは密行列になりやすく計算量が重い（$O(N^2)$ になりやすい）。そこで実装としては

• 内部をFEM、外部をBEMで扱うハイブリッド
• BEM行列の高速化（FMMやH行列） などが採用されることがある。

形状が複雑で、反磁界が支配的（永久磁石、磁性粒子集合、微細加工構造）な問題では、有力な選択肢となる。

7. FMM・階層法によるO(N)近似

磁荷（または双極子）相互作用を遠方では多重極展開でまとめて近似し、近傍のみを直接計算することで $O(N)$ あるいは $O(NlogN)$ を目指すのがFMM系の発想である。

概略は次である。

• 空間を階層的に分割（木構造）
• 遠方セル群の寄与を多重極展開で集約
• 近傍セルは直接和で計算
• 精度パラメータ（展開次数、近傍範囲）で誤差と計算量を制御

高速多重極（FFT）法が直交格子に強いのに対し、FMMは点群・不規則配置・適応メッシュとも相性がよい。ただし実装の複雑さと定数項が大きい点が課題となりやすい。

8. 精度検証

反磁界は長距離相互作用ゆえに、わずかな境界・カーネル誤差が全体に波及する。代表的な検証観点は次である。

8.1 解析解との比較

• 一様磁化楕円体の反磁界（反磁界係数が一様）
• 薄膜の近似式（面内・面外磁化での極限挙動）
• 標準問題（μMAGなど）でのエネルギー比較

8.2 エネルギーの整合性

反磁界エネルギーは

$$E_d=-\frac{{\mu }_0}{2}{\int }_{\mathrm{\Omega }}\mathbf{M}\cdot {\mathbf{H}}_{d}dV$$

の形で評価されることが多い。時間発展でエネルギーが不自然に増減する場合、時間積分の問題だけでなく反磁界計算のラップアラウンドや単位系混在も疑うべきである。

8.3 メッシュ依存性と境界依存性

• セルサイズを変えたときに、反磁界エネルギーが収束するか
• パディング量を増やしたときに、境界近傍の反磁界が安定するか
• PBCの有無で、期待する極限（薄膜、バルク）へ向かうか

9. 手法の比較

手法	典型の計算量	強い状況	注意点
直接和	$O(N^2)$	ごく小規模、検証用	大規模では不可
FFT畳み込み（FD）	$O(NlogN)$	直交格子、大規模、GPU	開境界はパディング、形状は階段誤差
FEM（Poisson/Laplace）	反復解法に依存	複雑形状、曲面、局所改良	外部無限領域処理が要点
BEM / FEM-BEM	問題依存（重め）	無限領域を厳密に扱いたい	密行列の高速化が必要
FMM / 階層法	$O(N)〜O(NlogN)$	不規則配置、適応メッシュ	実装が難しく定数項が大きいことがある
多層畳み込み	$O(NlogN)$系	薄膜多層、スペーサが多い系	近似条件と実装依存がある

10. 実装と再現性

• 単位系：$M（A/m）$、$H（A/m）$、$B（T）$、エネルギー密度（$J/m^3$）を混在させない
• カーネル定義：$N$の自己項、対称性、セル中心/セル平均の定義を明示する
• 境界条件：開境界（パディング量）とPBC（周期次元）を明示する
• 検証：一様磁化形状、μMAG標準問題、格子収束で最低限の整合性を確認する
• 形状誤差：階段近似が支配的な場合、FEMや補正の導入を検討する

まとめ

反磁界は磁化分布の全セル間結合として現れ、計算コストと境界条件の扱いが結果の信頼性を左右する。直交格子では離散反磁界テンソルのFFT畳み込みが標準であり、開境界ならパディング、周期系ならPBC対応カーネルが要点となる。複雑形状や無限領域の厳密性が重要な場合はFEM/BEMや階層法が有効であり、解析解・標準問題・格子収束による検証を通じて反磁界計算の妥当性を担保する必要がある。