{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 核磁気共鳴（NMR）と核四重極共鳴（NQR） \n", "\n", "（初版：2020年3月、更新：2023年2月21日） \n", "\n", "原子核は電気的にも磁気的にも中性ではない。外場と原子核の間に生じる相互作用として、下記の3つがある。 \n", "\n", "　・ゼーマン相互作用 \n", "　・核四重極相互作用 \n", "　・超微細相互作用 \n", "\n", "### ゼーマン相互作用 \n", "\n", "　核磁気モーメント $\\boldsymbol{\\mu}_n$ と外部磁場 $\\boldsymbol{H}_0$ による、磁気的な相互作用。 \n", "\n", "### 核四重極相互作用 \n", "\n", "　核四重極モーメント $eQ$ と核周りの電場勾配 $V_{\\alpha \\alpha}$ $(\\alpha = x, y, z)$ による、電気的な相互作用。\n", " \n", "\n", "### 超微細相互作用 \n", "\n", "　対象原子核とその周り（他の原子核や電子など）の間に生じる相互作用。 \n", "　物質の性質を調べる上で非常に重要で、実験的にはこれが知りたい。ここでは触れません。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "　\n", "# ゼーマン相互作用\n", "ゼーマン相互作用のハミルトニアンは以下のように書ける。 \n", "\n", "　　$\\mathcal{H}_Z = -\\boldsymbol{\\mu}_n \\cdot \\boldsymbol{H}_0 = -\\gamma_n h \\boldsymbol{I} \\cdot \\boldsymbol{H}_0$　　\n", "\n", "ここで、$\\boldsymbol{I}$ は核スピンである。 \n", "\n", "磁場方向を $z$ 軸とすると、ハミルトニアンは \n", "\n", "　　$\\mathcal{H}_Z = -\\gamma_n h H_0 \\boldsymbol{I}_z$　　\n", "\n", "と書け、固有エネルギーは \n", "\n", "　　$E_m = - \\gamma_n h H_0 m$ \n", "\n", "となる。$m$ は磁気量子数である。 \n", "$I = \\frac{3}{2}$ （ $m = 3/2 \\sim -3/2$ ）の場合の行列表示は、 \n", "\n", "　　$ \\mathcal{H}_Z = -\\gamma_n h H_0\n", " \\left[\\begin{array}{cccc}\n", " \\frac{3}{2} & 0 & 0 & 0 \\\\\n", " 0 & \\frac{1}{2} & 0 & 0 \\\\\n", " 0 & 0 & -\\frac{1}{2} & 0 \\\\\n", " 0 & 0 & 0 & -\\frac{3}{2} \\\\\n", " \\end{array}\\right]\n", "$\n", "\n", "と書け、対角行列なのでそれぞれの値が固有値（固有エネルギー）となる。 \n", "隣接準位間のエネルギー差は \n", "\n", "　　$\\Delta E = \\gamma_n h H_0$ \n", "\n", "となり、等間隔である。 \n", "よって、共鳴周波数は以下のようになる。 \n", "\n", "　　$\\nu =\\Delta E/ h = \\gamma_n H_0$ \n", "\n", "これを、julia上で書くと以下のようになる。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 10.0\n", " 10.0\n", " 10.0" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "γn = 10.0; #核磁気回転比（MHz/T）\n", "h = 1.0; #プランク定数。共鳴周波数にするときにhで割るので、１にしておく。\n", "H0 = 1.0; #外部磁場（T）\n", "IZ = [+3/2 0 0 0;\n", " 0 +1/2 0 0;\n", " 0 0 -1/2 0;\n", " 0 0 0 -3/2]\n", "Ene = -γn*h*H0*IZ # 固有エネルギーの行列要素\n", "ΔE = [Ene[1,1], Ene[2,2], Ene[3,3], Ene[4,4]] #エネルギー準位（MHz単位）\n", "ν = [ΔE[4]-ΔE[3], ΔE[3]-ΔE[2], ΔE[2]-ΔE[1]] #共鳴周波数（MHz）" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "　\n", "# 核四重極相互作用\n", "核四重極相互作用のハミルトニアンは以下のように書ける。 \n", "\n", "　　$\\mathcal{H}_Q = \\frac{h\\nu_Q}{6}[ 3I_z^2 - I^2 +\\frac{\\eta}{2}(I_+^2 + I_-^2) ]$ \n", "\n", "ここで、$\\nu_Q \\equiv 3eQV_{zz}/2I(2I-1)h$ は核四重極モーメント、$\\eta \\equiv (V_{xx}-V_{yy})/V_{zz}$ は非対称パラメータである。 \n", "$\\eta \\neq 0$ の場合、ハミルトニアンの行列要素に非対角成分があるため、対角化して固有エネルギーを求める必要がある。 \n", "\n", "$I = \\frac{3}{2}$の時のハミルトニアンの行列要素は、 \n", "\n", "　　$ \\mathcal{H}_Q = h\\nu_Q\n", " \\left[\\begin{array}{cccc}\n", " \\frac{1}{2} & 0 & \\frac{\\eta}{2\\sqrt{3}} & 0 \\\\\n", " 0 & -\\frac{1}{2} & 0 & \\frac{\\eta}{2\\sqrt{3}} \\\\\n", " \\frac{\\eta}{2\\sqrt{3}} & 0 & -\\frac{1}{2} & 0 \\\\\n", " 0 & \\frac{\\eta}{2\\sqrt{3}} & 0 & \\frac{1}{2} \\\\\n", " \\end{array}\\right]\n", "$\n", "\n", "となる。この場合は解析的に対角化できて、固有エネルギーは \n", "\n", "　　$E_{\\pm\\frac{3}{2}} = \\frac{h\\nu_Q}{2}\\sqrt{1+\\frac{\\eta^2}{3}}$、$E_{\\pm\\frac{1}{2}} = -\\frac{h\\nu_Q}{2}\\sqrt{1+\\frac{\\eta^2}{3}}$ \n", "\n", "となる。\n", "隣接準位間のエネルギー差は \n", "\n", "　　$\\Delta E = h\\nu_Q\\sqrt{1+\\frac{\\eta^2}{3}}$ \n", "\n", "であるので、共鳴周波数 $\\nu$ は \n", "\n", "　　$\\nu =\\Delta E/ h = \\nu_Q\\sqrt{1+\\frac{\\eta^2}{3}}$ \n", "\n", "となる。 \n", "\n", "これを、julia上で書くと以下のようになる。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×4 Matrix{Float64}:\n", " 5.0 0.0 0.0 0.0\n", " 0.0 -5.0 0.0 0.0\n", " 0.0 0.0 -5.0 0.0\n", " 0.0 0.0 0.0 5.0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "νQ = 10.0; #核四重極周波数（MHz）\n", "η = 0.0; #非対称パラメータ（0≦|η|≦1）\n", "h = 1.0; #プランク定数。共鳴周波数にするときにhで割るので、１にしておく。\n", "Hami = h*νQ*[1/2 0 η/(2*sqrt(3)) 0;\n", " 0 -1/2 0 η/(2*sqrt(3));\n", " η/(2*sqrt(3)) 0 -1/2 0;\n", " 0 η/(2*sqrt(3)) 0 1/2] #ハミルトニアンの行列要素" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n", "values:\n", "4-element Vector{Float64}:\n", " -5.0\n", " -5.0\n", " 5.0\n", " 5.0\n", "vectors:\n", "4×4 Matrix{Float64}:\n", " 0.0 0.0 1.0 0.0\n", " 1.0 0.0 0.0 0.0\n", " 0.0 1.0 0.0 0.0\n", " 0.0 0.0 0.0 1.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "Eval, Evec = eigen(Hami) #固有値と固有行列を求める" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 10.0\n", " 10.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ν = [Eval[3]-Eval[1], Eval[4]-Eval[2]] #共鳴周波数（MHz）" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$I = \\frac{5}{2}$ などについても、$\\mathcal{H}_Q$ の行列要素を求め対角化することにより、共鳴周波数を求めることができる。 \n", "$|\\pm m\\rangle$ のエネルギーが縮退するため、$I$ が半整数の場合、共鳴線の数は $(I-1/2)$ 本となる。 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "　\n", "# ゼーマン相互作用＋核四重極相互作用\n", "前述の2つの相互作用がある場合の共鳴周波数を求める。 \n", "「 $I = \\frac{3}{2}$ で $H_0 \\parallel V_{zz}$ 」の場合を考える。 \n", "ハミルトニアンは \n", "\n", "　　$\\mathcal{H} = \\mathcal{H}_Z + \\mathcal{H}_Q = -\\gamma_n h H_0 \\boldsymbol{I}_z + \\frac{h\\nu_Q}{6}[ 3I_z^2 - I^2 +\\frac{\\eta}{2}(I_+^2 + I_-^2) ]$ \n", "\n", "であり、行列要素で書くと \n", "\n", "　　$ \\mathcal{H} = \n", " \\left[\\begin{array}{cccc}\n", " -\\frac{3}{2}\\gamma_n h H_0 + \\frac{1}{2}h\\nu_Q & 0 & \\frac{\\eta}{2\\sqrt{3}} & 0 \\\\\n", " 0 & -\\frac{1}{2}\\gamma_n h H_0 -\\frac{1}{2}h\\nu_Q & 0 & \\frac{\\eta}{2\\sqrt{3}} \\\\\n", " \\frac{\\eta}{2\\sqrt{3}} & 0 & \\frac{1}{2}\\gamma_n h H_0 -\\frac{1}{2}h\\nu_Q & 0 \\\\\n", " 0 & \\frac{\\eta}{2\\sqrt{3}} & 0 & \\frac{3}{2}\\gamma_n h H_0 + \\frac{1}{2}h\\nu_Q\\\\\n", " \\end{array}\\right]\n", "$ \n", "\n", "となる。対角化して固有エネルギーを求めて、その差から共鳴周波数が計算できる。 \n", "これを、julia上で書くと以下のようになる。" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×4 Matrix{Float64}:\n", " -145.0 0.0 0.0 0.0\n", " 0.0 -55.0 0.0 0.0\n", " 0.0 0.0 45.0 0.0\n", " 0.0 0.0 0.0 155.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "γn = 10.0; #核磁気回転比（MHz/T）\n", "H0 = 10.0; #外部磁場（T）\n", "νQ = 10.0; #核四重極周波数（MHz）\n", "η = 0.0; #非対称パラメータ（0≦|η|≦1）\n", "h = 1.0; #プランク定数。共鳴周波数にするときにhで割るので、１にしておく。\n", "\n", "HamiZ = -γn*h*H0*[+3/2 0 0 0;\n", " 0 +1/2 0 0;\n", " 0 0 -1/2 0;\n", " 0 0 0 -3/2] #ゼーマン相互作用\n", "\n", "HamiQ = [h*νQ/2 0 h*η*νQ/(2*sqrt(3)) 0;\n", " 0 -h*νQ/2 0 h*η*νQ/(2*sqrt(3));\n", " h*η*νQ/(2*sqrt(3)) 0 -h*νQ/2 0;\n", " 0 h*η*νQ/(2*sqrt(3)) 0 h*νQ/2] #核四重極相互作用\n", "\n", "Hami = HamiZ + HamiQ #全体のハミルトニアン" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n", "values:\n", "4-element Vector{Float64}:\n", " -145.0\n", " -55.0\n", " 45.0\n", " 155.0\n", "vectors:\n", "4×4 Matrix{Float64}:\n", " 1.0 0.0 0.0 0.0\n", " 0.0 1.0 0.0 0.0\n", " 0.0 0.0 1.0 0.0\n", " 0.0 0.0 0.0 1.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra #前に実行してあれば、もう一度読み込む必要は無い\n", "Eval, Evec = eigen(Hami) #固有値と固有行列を求める" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "これまでは適当に固有値を引き算して共鳴周波数を求めていたが、色々と複雑になると真面目に遷移確率を計算する必要が出てくる。 \n", "$|m\\rangle \\leftrightarrow |m'\\rangle$ の遷移確率 $I_{m,m'}$ は \n", "\n", "　　$I_{m,m'} \\propto |\\langle m|I_x |m'\\rangle|^2$ \n", "\n", "であるので、固有ベクトルと演算子 $I_x$ から求めることができる。 \n", "$I = \\frac{3}{2}$ の時の $I_x$ の行列要素は \n", "\n", "　　$ I_x = \n", " \\left[\\begin{array}{cccc}\n", " 0 & \\frac{\\sqrt{3}}{2} & 0 & 0 \\\\\n", " \\frac{\\sqrt{3}}{2} & 0 & 1 & 0 \\\\\n", " 0 & 1 & 0 & \\frac{\\sqrt{3}}{2} \\\\\n", " 0 & 0 & \\frac{\\sqrt{3}}{2} & 0 \\\\\n", " \\end{array}\\right]\n", "$ \n", "\n", "である。\n", "先に計算した固有値と固有ベクトルについて共鳴周波数と遷移確率は以下のようになる。" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "freq. = 90.0 MHz, prob. = 0.7499999999999999\n", "freq. = 190.0 MHz, prob. = 0.0\n", "freq. = 300.0 MHz, prob. = 0.0\n", "freq. = 100.0 MHz, prob. = 1.0\n", "freq. = 210.0 MHz, prob. = 0.0\n", "freq. = 110.0 MHz, prob. = 0.7499999999999999\n" ] } ], "source": [ "Ix = [0 sqrt(3)/2 0 0;\n", " sqrt(3)/2 0 1 0;\n", " 0 1 0 sqrt(3)/2;\n", " 0 0 sqrt(3)/2 0]\n", "\n", "#共鳴周波数の計算\n", "f1 = abs(Eval[1] - Eval[2])\n", "f2 = abs(Eval[1] - Eval[3])\n", "f3 = abs(Eval[1] - Eval[4])\n", "f4 = abs(Eval[2] - Eval[3])\n", "f5 = abs(Eval[2] - Eval[4])\n", "f6 = abs(Eval[3] - Eval[4])\n", "#遷移確率の計算\n", "p1 = (Evec[:,1]'*Ix*Evec[:,2])^2\n", "p2 = (Evec[:,1]'*Ix*Evec[:,3])^2\n", "p3 = (Evec[:,1]'*Ix*Evec[:,4])^2\n", "p4 = (Evec[:,2]'*Ix*Evec[:,3])^2\n", "p5 = (Evec[:,2]'*Ix*Evec[:,4])^2\n", "p6 = (Evec[:,3]'*Ix*Evec[:,4])^2\n", "\n", "println(\"freq. = \", f1, \" MHz, prob. = \", p1)\n", "println(\"freq. = \", f2, \" MHz, prob. = \", p2)\n", "println(\"freq. = \", f3, \" MHz, prob. = \", p3)\n", "println(\"freq. = \", f4, \" MHz, prob. = \", p4)\n", "println(\"freq. = \", f5, \" MHz, prob. = \", p5)\n", "println(\"freq. = \", f6, \" MHz, prob. = \", p6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ガウス関数を畳み込んだ共鳴スペクトルは以下のようになる。" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "gr()\n", "\n", "freq0 = [f1, f2, f3, f4, f5, f6];\n", "prob0 = [p1, p2, p3, p4, p5, p6];\n", "freq1 = freq0[prob0 .> 0] #遷移確率ゼロのところは消しとく\n", "prob1 = prob0[prob0 .> 0] #遷移確率ゼロのところは消しとく\n", "\n", "n = length(freq1) #ピークの数\n", "w0 = 1 #線幅に関する定数\n", "divf = 500 #周波数の分割数\n", "fmin = minimum(freq1)*0.9 #計算する周波数範囲（下限）\n", "fmax = maximum(freq1)*1.1 #計算する周波数範囲（上限）\n", "\n", "calcfreq = (fmin:(fmax-fmin)/divf:fmax)\n", "calcint = calcfreq*0\n", "for x in 1:n\n", " calcint = calcint + prob1[x]*exp.(-(calcfreq .- freq1[x]).^2 ./ w0)\n", "end\n", "\n", "plot(calcfreq, calcint, size = (400, 200))\n", "plot!(xlabel = \"Frequency (MHz)\", ylabel = \"Intensity\", legend = false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "上記はなるべく何を計算しているのかを分かりやすく書いてあります。 \n", "ガチで計算するときには巧く書いて、利便性をあげた方が良いです。" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.9.0-beta3", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.0" } }, "nbformat": 4, "nbformat_minor": 4 }