LoginSignup
22
21

More than 5 years have passed since last update.

@popondeli

イジングモデル一次相転移シミュレーション①

Posted at

反強磁性-強磁性磁気一次相転移のシミュレーション

室温では磁石としての性質を持っている(くっつく)が、温度下げていくとある温度で
磁石としての性質を持たなくなる(くっつかなくなる)物質が存在する。
この現象を磁気相転移といい、今回はこれを数値計算によりシミュレートを試みる。

磁性体のミクロなイメージ

各原子はスピン(下記図の矢印)というものをもっており、スピンが一方に偏って揃っている状態(aやc)は磁石状態(くっつく）となり、
ばらばらになってたり(d)逆向きにそろって打ち消し合っている状態(b)では非磁石状態（くっつかない)となる。

image.png

引用：http://www.nims.go.jp/mmu/tutorials/magnetism_j.html

スピン間の相互作用

スピン間の相互作用は交換相互作用($J$)と呼ばれ、$J>0$のときは同じ向きにスピンは結合しようし、
$J<0$のときは逆向きに揃おうとする。
$|J|$が結合の強さを表している。

スピンベクトルを$\vec{s}$とすると、スピン間の交換相互作用によるエネルギーは

E=-2J\vec{s_i}・\vec{s_j}

と表され、$J$が正の時は同じ向きに向いていたほうがエネルギーが低く、$J$が負のときは逆向きに向いた
ほうがエネルギーが低くなり、磁気エネルギー的に安定となる。

一方、外部から与えられるエネルギー(熱エネルギー等)は、スピン同士をバラバラにしようとする働きをし、
例えば、交換相互作用エネルギー＜熱エネルギー程度になると上図の(d)常磁性のようになる。

一般の磁石は温度を上げていくと熱エネルギー増加により(a)強磁性⇒(d)常磁性へ磁気二次相転移を起こす。
しかし、昇温に伴い(b)反強磁性⇒(a)強磁性へ磁気一次相転移、(a)強磁性⇒(d)常磁性へ磁気二次相転移
と二段階相転移を起こす物質があり、今回はそれをシミュレーション対象とする。

■縦軸を磁石としての強さ(磁化)、横軸を温度とした時のイメージ

  • (a)強磁性⇒(d)常磁性

image.png
引用：https://astamuse.com/ja/published/JP/No/2016532778

自由エネルギー

どちらの状態が実現するか(例えば強磁性or反強磁性)を定量的に求めるために、ヘルムホルツの自由エネルギー$F$を扱う。
各状態(強磁性、反強磁性等)毎に自由エネルギーが決まり、自由エネルギーが低い状態のほうが確率的に実現しやすいため、
時間経過に従い自由エネルギーが低い状態へ相転移していく。

ヘルムホルツの自由エネルギーは

F=U-TS

で与えられ、$U$は内部エネルギー、$T$は温度、$S$はエントロピーである。

一次相転移

自由エネルギーを温度等の物理量で一階偏微分した物理量(温度微分の場合エントロピー)が、
相転移時に不連続変化する場合一次相転移という。
一次相転移の場合、安定状態と準安定状態が存在し得、下図でいうと②が安定状態、①が準安定状態である。
準安定状態から安定状態へ移るにはエネルギー壁$E_{12}$(活性化エネルギー)程度の熱エネルギーが
与えられなければならない。

image.png
引用：http://www.c.u-tokyo.ac.jp/info/news/topics/20150203150801.html

反応速度論では、①⇒②状態へ移る速度$v$は$e^{-\frac{E_{12}}{k_BT}}$に比例する。
$k_B$はボルツマン定数、$T$は温度であり、$k_BT$は外部の熱エネルギーに相当する。
$k_BT(熱エネルギー)<<E_{12}(活性化エネルギー)$の時、$v⇒0$となり相転移はほとんど進まない。
(ポテンシャルを乗り越えるための外部エネルギーが十分でないため)

Kittelのエクスチェンジインバージョンモデル

反強磁性-強磁性磁気一次相転移のモデルとして、Kittelのエクスチェンジインバージョンモデルがある。

エクスチェンジインバージョンモデルでは、スピン間の相互作用$J$が原子間距離に対して線形的に変化し、
ある臨界原子間距離で$J$の符号が変わる。

J(c)=ρ(c-c_C)

ρは比例係数、cは原子間距離、$c_C$は臨界原子間距離である。

降温に伴い、物質が熱収縮し原子間距離が縮まることで$J$の符号が変わり相転移が起こることをモデルしている。

原子i,jの自由エネルギーは

F=\frac{1}{2}R(c-c_T(T))^2-ρ(c-c_C)\vec{s_i}・\vec{s_j}-\vec{H}・(\vec{s_i}+\vec{s_j})…①

で与えられる。
$R$は弾性係数、$H$は磁場、$c_T(T)$は磁気的な相互作用がない時の平衡原子間距離である。
第一項は弾性エネルギー、第二項は交換相互作用エネルギー、第三項はゼーマンエネルギーを表している。
ゼーマンエネルギーはスピンを磁場と同じ方向に向かせようとする働きをする。
今回はゼーマンエネルギー項は省略する。($H=0$)

自由エネルギーを原子間距離で偏微分し0とおき、平衡時の原子間距離を求めると

\frac{∂F}{∂c}=R(c_{eq}-c_T(T))-ρs^2cosφ=0\\
c_{eq}=c_T(T)+\frac{ρs^2}{R}cosφ…②

φはスピンiとjのなす角である。
第二項が磁気的な相互作用による膨張/縮小を表しており、強磁性状態(φ=0)と、反強磁性状態(φ=π)では$2\frac{ρs^2}{R}$の差がある。
(=相転移時不連続な格子定数変化が生じる=一次相転移)

②を①に代入すると

F=-\frac{ρ^2}{2R}s^4cos^2φ-ρ(c_T(T)-c_C)s^2cosφ…③

$cos^2φ$の項が入ってくることで、φが0～πの間で上図のようなポテンシャル壁が導入される。(=準安定状態出現)
強磁性状態(φ=0)での自由エネルギー$F_{FR}$と、反強磁性状態(φ=π)での自由エネルギー$F_{AFM}$は

F_{FR}(T)=-\frac{ρ^2}{2R}s^4-ρ(c_T(T)-c_C)s^2…④\\
F_{AFM}(T)=-\frac{ρ^2}{2R}s^4+ρ(c_T(T)-c_C)s^2…⑤

Kittelのモデルでは2つのスピン間のエネルギーを考慮しているが、今回のシミュレーションでは
対象とするスピンと、隣接スピン4個間でのエネルギーを考慮するよう拡張する。
対象とするスピンを$i$とすると、スピン$i$と隣接4個のスピン間での自由エネルギーは
①を拡張して、

F_i=\frac{1}{2}R(c-c_T(T))^2-ρ(c-c_C)\sum_j\vec{s_i}・\vec{s_j}\\
=\frac{1}{2}R(c-c_T(T))^2-ρ(c-c_C)s^2\sum_jcos(φ_i-φ_j)…⑥

$φ_i$はスピンiの向き(↑なら0、↓ならπ)
$j$は隣接4個に対してとる。

②、③と同様に

c_{eq}=c_T(T)+\frac{ρs^2}{R}\sum_jcos(φ_i-φ_j)…⑦\\
F_i=-\frac{ρ^2}{2R}s^4(\sum_jcos(φ_i-φ_j))^2-ρ(c_T(T)-c_C)s^2\sum_jcos(φ_i-φ_j)…⑧

磁気的な相互作用がない時の平衡格子定数$c_T(T)$は温度に対して比例するとする。

c_T(T)=αT+β…⑨

また、上記の交換相互作用線形変化の仮定では、原子間距離が離れて行くにつれて交換相互作用の
絶対値が大きくなってしまい、高温側での強磁性⇒常磁性相転移の振る舞いが正しくシミュレート
できないことが予想されるので

J(c)=ρ(c-c_C)e^{-γ(c-(c_C+ζ))}…⑩

上記のような減衰項を独自に追加したいが、$\frac{∂F}{∂c}=0$を満たす$c$を解析的に求めるのが難しいため、
高温側($T>T_T+ζ$)では交換相互作用$J$は正の定数として計算する。

F_i=-Js^2\sum_jcos(φ_i-φ_j), (T≧T_T+ζ)…⑪\\
F_i=-\frac{ρ^2}{2R}s^4(\sum_jcos(φ_i-φ_j))^2-ρ(c_T(T)-c_C)s^2\sum_jcos(φ_i-φ_j), (T<T_T+ζ)…⑫

・$J(c)∝(c-1)e^{(-(c-(1+2)))}$のグラフ(交換相互作用として設定したかった形)
image.png

しかし、交換相互作用を定数項にするとうまく相転移を再現できなかったため、高温側($T>T_T+ζ$)では自由エネルギー全体に減衰項をかけ、$T>>T_T+ζ$$F_i⇒0$(磁気的相互作用がない状態)となるように設定することにした。

F_i=e^{-γ(ρ(c_T(T)+\frac{ρs^2}{R}\sum_jcos(φ_i-φ_j)-(c_C+ζ)))}(-\frac{ρ^2}{2R}s^4(\sum_jcos(φ_i-φ_j))^2-ρ(c_T(T)-c_C)s^2\sum_jcos(φ_i-φ_j)), (T≧T_T+ζ)…⑬\\
F_i=-\frac{ρ^2}{2R}s^4(\sum_jcos(φ_i-φ_j))^2-ρ(c_T(T)-c_C)s^2\sum_jcos(φ_i-φ_j), (T<T_T+ζ)…⑭

相転移シミュレーションアルゴリズム

以下の流れでスピン状態を計算する。

  1. 温度$T$を設定する。(=原子間距離が変わる=交換相互作用が変わる)
  2. ランダムに一つスピンを選び、⑬、⑭式に従い上向き、下向き状態それぞれ自由エネルギー$F_↑$、$F_↓$を計算する。
  3. 確率$\frac{exp(-\frac{F_↑}{k_BT})}{exp(-\frac{F_↑}{k_BT})+exp(-\frac{F_↓}{k_BT})}$で↑状態、$\frac{exp(-\frac{F_↓}{k_BT})}{exp(-\frac{F_↑}{k_BT})+exp(-\frac{F_↓}{k_BT})}$で↓状態を次の状態候補とする。
  4. 一次相転移の効果を取り入れるため、3で次の状態の候補が、今の状態(自由エネルギー$F_c$とする)と異なる場合、$\frac{exp(-\frac{F_{max}-F_c}{k_BT})}{\int_0^∞exp(-\frac{F_{max}-F_c}{k_BT})d(F_{max}-F_c)}$の確率で次の状態への遷移を確定する。ここで$F_{max}$は、上の自由エネルギーの図でいう極大値($E_{12}$の場所)である。例えば、今の状態が①、次の状態が②だと$exp(-\frac{F_{max}-F_c}{k_BT})=exp(-\frac{E_{12}}{k_BT})$である。$E_{12}>>k_BT$の場合、高確率で遷移はキャンセルされる。横軸はスピンの向き$φ$に対応させる。$F(φ)>F_c$となるような$φ$がない場合は、3の遷移を確定する。
  5. 2～4のステップをある程度磁化が収束するまで繰り返す。
  6. 1～5を降温/昇温させながら各温度で実施。

$F_{max}$の計算法としては、現在のスピンの向きをπ/N刻みでπ回転させながら⑦を計算し、$F(φ)>F_c$をみたし
最大となる$F(φ)$を$F_{max}$として求める。

また、今回は厳密な定量性にはこだわらず、各種パラメータは相転移が確認できる程度に適当に調整して選ぶ。

実装

import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime as dt
import os

class SpinSystem:
    def __init__(self, spin_num_col, spin_num_row, temperature, magnetic_field, alpha, beta, rho, r, s, gamma, zeta, kb, cc):
        self.spin_num_col = spin_num_col
        self.spin_num_row = spin_num_row
        self.temperature = temperature
        self.magnetic_field = magnetic_field
        self.alpha = alpha
        self.beta = beta
        self.rho = rho
        self.r = r
        self.s = s
        self.gamma = gamma
        self.zeta = zeta
        self.spin_map = self.create_spinmap(spin_num_row, spin_num_col)
        self.kb = kb
        self.cc = cc

    # 温度の設定
    def set_temperature(self, temperature):
        self.temperature = temperature

    # 温度の取得
    def get_temperature(self):
        return self.temperature

    # 磁場の設定
    def set_magnetic_field(self, magnetic_field):
        self.magnetic_field = magnetic_field

    # 温度Tにおけるc_Tを取得
    def get_elastic_lattice_parameter(self, temperature):
        return max(0.0, self.alpha * self.temperature + self.beta)

    # 温度Tにおけるc_eqを取得
    def get_equilibrium_lattice_parameter(self, i, j):
        try:
            cos_sum = np.cos(np.pi*(self.spin_map[i][j] - self.spin_map[i-1][j])) + np.cos(np.pi*(self.spin_map[i][j] - self.spin_map[i+1][j])) + np.cos(np.pi*(self.spin_map[i][j] - self.spin_map[i][j+1])) + np.cos(np.pi*(self.spin_map[i][j] - self.spin_map[i][j-1]))
        except Exception as e:
            print("Out of index. {}, {}".format(i, j))
            cos_sum = 0.0
        return self.get_elastic_lattice_parameter(self.temperature) + ((self.rho * self.s**2) / self.r) * cos_sum

    # スピン(i, j)がφ方向を向いているときの自由エネルギー計算
    def calc_free_energy(self, i, j, phi):
        try:
            cos_sum = np.cos(phi - np.pi*self.spin_map[i-1][j]) + np.cos(phi - np.pi*self.spin_map[i+1][j]) + np.cos(phi - np.pi*self.spin_map[i][j+1]) + np.cos(phi - np.pi*self.spin_map[i][j-1])
        except Exception as e:
            print("Out of index. {}, {}".format(i, j))
            cos_sum = 0.0
        # 高温側では減衰項を掛ける
        if self.temperature < 400:
            energy = - ((self.rho**2 * s**4) / (2.0 * self.r)) * cos_sum**2 - self.rho * (self.get_elastic_lattice_parameter(self.temperature) - self.cc) * self.s**2 * cos_sum
        else:
            J = self.rho * (self.get_elastic_lattice_parameter(400) - self.cc) * np.exp(-0.01 * (self.temperature - 400))
            energy = - J * self.s**2 * cos_sum
        return energy

    # (i, j)のスピンを更新
    def update_spin(self, i, j):
        up_prob = np.exp(-self.calc_free_energy(i, j, 0)/(self.kb*self.temperature)) / (np.exp(-self.calc_free_energy(i, j, 0)/(self.kb*self.temperature)) + np.exp(-self.calc_free_energy(i, j, np.pi)/(self.kb*self.temperature)))
        rand = np.random.rand()
        # 現在↑向き&次の状態も↑、もしくは現在↓向き&次の状態も↓のため何もしない
        if (rand <= up_prob and self.spin_map[i][j] == 0) or (rand > up_prob and self.spin_map[i][j] == 1):
            return

        # 低温側(400K以下)では一次相転移効果を取り入れる
        if self.temperature < 400:
            # スピンをπ回転させて自由エネルギーの極大値を計算
            cur_free_energy = max_free_energy = self.calc_free_energy(i, j, np.pi*self.spin_map[i][j])
            N = 100
            dphi = np.pi / N
            for n in range(N):
                f = self.calc_free_energy(i, j, np.pi*self.spin_map[i][j] + (n+1)*dphi)
                if f > max_free_energy:
                    max_free_energy = f
            # 自由エネルギー極大値が存在する場合
            if max_free_energy > cur_free_energy:
                # ポテンシャル壁を乗り越えられなかった場合は何もしない
                if np.random.rand() >= np.exp(-(max_free_energy - cur_free_energy)/(self.kb*self.temperature)) / (self.kb * self.temperature):
                    return

        # 自由エネルギー極大値が存在しない　or ポテンシャル壁を乗り越えた場合、スピンをフリップする。
        self.spin_map[i][j] = 1 if self.spin_map[i][j] == 0 else 0

    # スピンをランダムにnum_spin個更新
    def update_spins(self, num_spin):
        for i, j in zip(np.random.randint(1, self.spin_num_row-2 ,num_spin), np.random.randint(1, self.spin_num_col-2 ,num_spin)):
            self.update_spin(i, j)

    # 規格化した磁化を取得
    def get_normalized_magnetization(self):
        return np.abs(sum([np.cos(np.pi * spin) for spin in self.spin_map.reshape(-1,)])) / len(self.spin_map.reshape(-1,))

    # スピンマップを取得
    def get_spinmap(self):
        return self.spin_map

    # spin_num_row×spin_num_colのスピンマップを作成
    def create_spinmap(self, spin_num_row, spin_num_col):
        spin_map = np.zeros([spin_num_row, spin_num_col])
        for i in range(spin_num_row):
            for j in range(spin_num_col):
                spin_map[i][j] = 0 if np.random.rand() <= 0.5 else 1
        return spin_map

if __name__ == '__main__':
    spin_num_col = 50
    spin_num_row = 50
    temperature = 10
    magnetic_field = 0
    alpha = 0.1
    beta = -1.0
    rho = 10.0
    r = 1.0
    s = 1.0
    gamma = 0.1
    zeta = 15.0
    kb = 0.5
    cc = 15.0 # T=150Kくらいで相転移が起きるように設定

    spin_system = SpinSystem(spin_num_col = spin_num_col, spin_num_row = spin_num_row, temperature = temperature, magnetic_field = magnetic_field, alpha = alpha, beta = beta, rho = rho, r = r, s = s, gamma = gamma, zeta = zeta, kb = kb, cc = cc)
    spin_system.update_spins(1000000)

    now = dt.now().strftime('%Y%m%d%H%M%S')
    os.mkdir(dt.now().strftime('%Y%m%d%H%M%S'))

    spin_map = spin_system.get_spinmap()
    plt.imshow(spin_map)
    plt.text(5, 5, 'T : {}K up, M : {}'.format(str(temperature), str(spin_system.get_normalized_magnetization())) , bbox={'facecolor': 'white', 'pad': 10})
    plt.savefig(now + '\\' + str(spin_system.get_temperature()) + 'K_up.png')

    delta_t = 5

    temperature_hist2 = []
    magnetization_hist2 = []
    temperature = 10.0
    # 10K⇒800Kまで5K刻みで温度を上げる
    while temperature <= 800:
        spin_system.set_temperature(temperature)
        spin_system.update_spins(3000000)
        temperature_hist2.append(temperature)
        magnetization = spin_system.get_normalized_magnetization()
        magnetization_hist2.append(magnetization)
        # スピンマップの書き出し
        if temperature % 5 == 0:
            spin_map = spin_system.get_spinmap()
            plt.clf()
            plt.imshow(spin_map)
            plt.text(5, 5, 'T : {}K up, M : {}'.format(str(temperature), str(spin_system.get_normalized_magnetization())) , bbox={'facecolor': 'white', 'pad': 10})
            plt.savefig(now + '\\' + str(int(temperature))+ 'K_up.png')
            print("Temperature : {}, Magnetization : {}".format(temperature, magnetization))
        temperature += delta_t

    temperature_hist = []
    magnetization_hist = []
    # 500K⇒10Kまで5K刻みで温度を下げる
    while temperature >= 10:
        spin_system.set_temperature(temperature)
        spin_system.update_spins(3000000)
        temperature_hist.append(temperature)
        magnetization = spin_system.get_normalized_magnetization()
        magnetization_hist.append(magnetization)
        # スピンマップの書き出し
        if temperature % 5 == 0:
            spin_map = spin_system.get_spinmap()
            plt.clf()
            plt.imshow(spin_map)
            plt.text(5, 5, 'T : {}K up, M : {}'.format(str(temperature), str(spin_system.get_normalized_magnetization())) , bbox={'facecolor': 'white', 'pad': 10})
            plt.savefig(now + '\\' + str(int(temperature)) + 'K_down.png')
            print("Temperature : {}, Magnetization : {}".format(temperature, magnetization))
        temperature -= delta_t

    plt.clf()
    plt.plot(temperature_hist, magnetization_hist, label="Decrease Temperature")
    plt.plot(temperature_hist2, magnetization_hist2, label="Increase Temperature")
    plt.legend() 
    plt.ylabel("Normalized Magnetization")
    plt.xlabel("Temperature")
    plt.show()

結果

磁化-温度曲線

  • 温度10K⇒800K⇒10K

更新数を変えて何度か実施。
緑が昇温過程、青が降温過程。

MagneticTransition_rho-10_kb-1_decay-500-3.png

奇麗でないが、昇温/降温で一次相転移特融のヒステリシスらしきものがみられる。

MagneticTransition_rho-10_kb-1_decay-500-4.png

MagneticTransition_rho-10_kb-1_decay-500-5.png

MagneticTransition_rho-10_kb-1_decay-500-5000000.png

更新数を多くしないと、反強磁性⇒強磁性の相転移がスムーズに進まないことが多い。
400Kで一気に相転移が進むのは、一次相転移効果は400K以下しか計算していないためと思われる。

  • 一次相転移効果計算省略(上記手順4スキップ)

MagneticTransition_rho-10_kb-1_decay-500-no-FOMT.png

ヒステリシスは消えた。

スピンマップ

  • 10K (反強磁性状態)

spin_map_10K.png

↑、↓のスピンがほぼ交互に並んでおり秩序状態であることが確認できる。

  • 400K (強磁性状態)

更新回数少な目バージョン
規格化磁化：0.2392

spin_map_400K_magnetization-0.2392.png

クラスタができて局所安定状態に落ち込んでほとんど相転移がとまってしまっている？

更新回数増やしたバージョン
規格化磁化：0.8592

spin_map_400K_magnetization-0.8592_iter-10000000.png

ほぼ強磁性への相転移は完了している。端に赤があるのは端の更新処理を行っていないため。

  • 800K (常磁性状態)

規格化磁化：0.02

spin_map_800K_magnetization-0.02_iter-10000000.png

反強磁性状態みたいに交互に並んでおらず無秩序状態であることが確認できる。

  • 昇温時変化過程（一次相転移効果あり）

y4T8s52aa.gif

反強磁性⇒強磁性相転移時にクラスタができ、Mが0.1くらいでほぼ進まなくなってしまっている。

  • 降温時変化過程（一次相転移効果あり）

rsXX0r0b1.gif

  • 昇温時変化過程（一次相転移効果なし）

VPWsoZ93e.gif

クラスタは多少できているがスムーズに強磁性状態へ完全転移している。

  • 降温時変化過程（一次相転移効果なし）

udxmDo1d6.gif

まとめ

  • 温度変化する交換相互作用を導入し、反強磁性-強磁性磁気一次相転移らしきものは再現させることができた。
  • 一次相転移効果を取り入れると、反強磁性-強磁性相転移がスムーズに進まないことがあるが、クラスタができて局所安定状態に落ち込んでからが特に進み方が遅いようにみられる。一度クラスタができるとスピンフリップはほとんどクラスタ境界でしか起きないため、完全相転移までにかなりの回数更新をかけないといけないと思われる。
  • 強磁性-常磁性相転移に関しては、高温側は自由エネルギーを指数関数的に減少させて両状態の自由エネルギーを0に近づけることで転移確率を等しくなるようにして計算したが、交換相互作用エネルギーと熱エネルギーの競合の意味での相転移になっていないため修正したい。
  • 高温側・低温側で自由エネルギーを別々に設けているため一つに統一したい。

磁気エントロピー

磁気エントロピーは、ミクロな視点ではスピンの揃っている度合いを表す指標と考えられる。
全部そろっている場合エントロピーは最小となり、バラバラになるほど大きくなっていく。
式で書くと、エントロピー$S$は

S(W)=k_BlogW

で表される。$k_B$はボルツマン定数、$W$は状態数である。

例えば、上下のみを向くスピンが横に4個並んでいる場合を考えると、全部上向きにそろっている
状態だと状態数$W$は1であるため$S(1)=0$となり、一つだけ上を向く状態は、一番左が上、左から二番目が上…
と4通りあるため状態数$W$は4となり$S(4)=k_Blog4>S(1)$となり、後者のほうがエントロピーが大きい。
自由エネルギーに$-TS$項があるが、高温になるとスピンがそろって交換相互作用項で稼ぐより、
スピンがバラバラになりエントロピー項で稼いだ方が利得が大きくなるため、スピンがバラバラ状態(常磁性)が実現する。

一方、マクロな視点では、磁気エントロピーは

S=-\frac{∂F_{mag}(T)}{∂T}

と表される。$F_{mag}$は自由エネルギーの中の磁気エネルギーである。

相転移時における各状態間のエントロピー差$|S_{FRI}(T_T)-S_{AFM}(T_T)|$は相転移の駆動力に影響し、
エントロピー差が大きくなるほど相転移は進みやすい。

各磁場下での磁化-温度曲線が与えられている場合、特定の温度、磁場における磁気エントロピー差$S(T,H)$は、
以下のマクスウェル方程式から計算できる。

ΔS(H,T)=\int_0^H\frac{dM}{dT}dH

$H$は磁場、$M$は磁化、$T$は温度である。
離散データに対しては以下で近似できる。

ΔS(H,T)≃\sum_H\frac{dM}{dT}ΔH

上記の手順で各磁場での$M-T$曲線を求めれば、エントロピー差はこの離散式から計算できる。

今後の課題

・元素置換していない状態では相転移は起きず、元素置換すると相転移が起きる物質が存在する。置換効果をシミュレーションに取り入れないか。
・磁場が印加された状態だと、低温で2つの状態が共存したまま相転移がほとんど進まなくなる現象がある。この現象をシミュレーションできないか。
・エクスチェンジインバージョンモデルとは別に、交換相互作用が原子間距離に対して符号を変えながら減衰振動する
伝導電子を介したRKKY相互作用があるが、この交換相互作用を導入してシミュレーションするとどうなるか。
https://ja.wikipedia.org/wiki/RKKY%E7%9B%B8%E4%BA%92%E4%BD%9C%E7%94%A8
・スピンの向きを上下2方向だけでなく$n$方向許容し、ヘリカルスピン状態等シミュレートできないか。
・単純格子ではなく、より複雑な結晶構造で、関数形が異なる複数種類の交換相互作用を導入してシミュレート
・フラストレーションが生じるように交換相互作用の関数形を設定してシミュレート
・3次元へ拡張

強磁性-反強磁性一次相転移関連論文

・Model of Exchange-Inversion Magnetization
https://journals.aps.org/pr/abstract/10.1103/PhysRev.120.335
・Theory of Antiferromagnetic-Ferromagnetic Transitions in Dilute Magnetic Alloys and in the Rare Earths
https://journals.aps.org/pr/abstract/10.1103/PhysRev.108.1233
・Interpretation of Magnetic Properties of Dysprosium
https://journals.aps.org/pr/abstract/10.1103/PhysRev.116.1464
・Magnetic Structure Transitions
https://journals.aps.org/pr/abstract/10.1103/PhysRev.90.55

22
21
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
22
21