MicrOMEGAsを用いて暗黒物質の残存量、散乱断面積を計算する #Python

　こんにちは、私は大学院の２年間で暗黒物質を研究をしていました。今回は(主に後輩向けに)暗黒物質のモデルから物理量を計算するツールを紹介していきたいと思います。暗黒物質の現象論の研究ではモデルから物理量を計算し、実験結果と比較します。その中には膨大な計算が必要なものもあり、プログラムによる数値計算が必須のものもあります。今回はその計算プログラムとしてFeynRulesとMicrOMEGAsを紹介します。
　この記事の1、2節では背景として、暗黒物質について得られている実験結果を紹介します。3節では本記事で扱う暗黒物質のモデル(実スカラーシングレットモデル)の紹介。4節では、FeynRulesを用いてそのモデルをコンピューターが理解できる形に直す方法(モデルファイルの作成法)を紹介します。5節ではMicrOMEGAsを用いてモデルファイルから物理量を計算し、実験結果との比較を行います。

#目次
1.ダークマターの量
 2.直接探索実験
 3.実スカラーシングレットモデル
 4.FeynRules
5.MicrOMEGAs

#1. ダークマターの量
　宇宙マイクロ波背景放射(宇宙初期に生じたおよそ３Kの黒体輻射)の温度揺らぎからダークマターの量を求めることが出来ています。この測定から下の円グラフのように、ダークマターのエネルギー密度が宇宙全体に対してどれくらいあるか分かりました。観測の結果ダークマターはおよそ26.8％、既知の物質は4.9％、宇宙を膨張させるダークエネルギーは68.3％存在することが分かります。これを宇宙論パラメータという物理量で表すとダークマターの存在比26.8％は$\Omega h^2=0.12$と表されます。後で詳しくお話ししますが、この結果からダークマターの理論に制限を与えることが出来ます。

[National Radio Astronomy Observatoryより引用]

#2. 直接探索実験
　他にもダークマターを直接的に検出しようという取り組みもあります。具体的にはXenon,　LZなどがすでに行われており、将来的にはより精度が高いDARWIN実験も行われます。この実験ではダークマターと原子核の散乱を検知します。現在、この実験ではダークマターの有力な信号は検知されていません。そのため、暗黒物質の理論を作る際は原子核とDMの反応率が十分小さくなるようにする必要があります。
　下のグラフがLZの実験結果です。縦軸が暗黒物質と核子との散乱断面積, 横軸が暗黒物質の質量です。黒線は直接探索実験LZにより与えられる暗黒物質-核子散乱断面積の上限値を表しています。この実験結果から、理論に対して「散乱断面積の理論計算がこの黒い線よりも下側にこなくてはならない」という制限が与えられます。

[LZ(Phys. Rev. Lett. 131, 041002 (2023))より引用]

#3. 実スカラーシングレットモデル
この記事では簡単な例として実スカラーシングレットモデル(J. M. Cline, K. Kainulainen, P. Scott and C. Weniger, Phys. Rev. D 88 (2013))を用います。このモデルは他もモデルよりも比較的シンプル(少ないインプットパラメータ)となっており、例としては最適です。
モデルについて具体的に紹介します。このモデルでは、標準模型のゲージ変換の下で変換しない実スカラー場$S$を考えます。また、この$S$については変換$S \rightarrow -S$の下で不変であるとしています。この$Z_2$対称性により$S$は安定な粒子となります($Z_2$対称性により$S$の相互作用は$Z_2$パリティーを保つものに限定されるため)。つまり、$S$はこのモデルにおける暗黒物質となります。また、真空期待値$\langle S \rangle=0$であると仮定します。この仮定の下でポテンシャルは以下のように表されます。
$$
V = \frac{\mu_S^2}{2}S^2+\frac{\lambda_{H S}}{2}S^2|H|^2+\frac{\lambda_S}{4!}S^4.
$$
ここで、$H$は標準模型のHiggs場です。次にこのポテンシャルのスカラー場を以下のようにパラメータ化します。
$$
H =
\begin{pmatrix}
\pi_{W^+}
\\1/\sqrt{2}(v+h+i\pi_{Z})
\end{pmatrix}.
$$
ここで、$\pi_{W^+}$、$\pi_{Z}$はそれぞれ、$W^+$、$Z$に対する南部-ゴールドストーン(NG)ボゾンです。$v$はHiggs場の真空期待値です。この下でスカラー場$S$と$h$に対する相互作用項は以下の様に与えられます。
$$
V\supset{\frac{m_S^2}{2} S^2+\frac{\lambda_{H S}}{2} v S^2 h+\frac{\lambda_{H S}}{4} S^2 h^2+\frac{\lambda_S}{4!} S^4}.
$$
ここで、$S$の質量$m_S=\sqrt{\mu_S^2+\lambda_{H S}\ v^2/2}$です。以上より、$\mu_S$は$m_S$と関係づいているので、このモデルのインプットパラメータを($m_S, \lambda_{H S}, \lambda_S$)とします。また、ポテンシャルにおける$\lambda_S$の相互作用項は、ダークマターの量、直接探索実験の現象論に影響を与えないため、以下の解析では$\lambda_S$に制限は与えられません。

#4. モデルファイルの作成
次にFeynRulesによりモデルファイルを作成します。これはモデルのラグランジアンを記入したファイル(frファイル)から、Calchepといった計算ソフトに読み込ませる用のファイル(Calchepファイル)を生成するプログラムです。FeynRulesはFeynRuleのページよりダウンロードして下さい。ダウンロードしたファイルを解凍すると、"feynrules-current"というフォルダが得られると思います。
　次にFeynRulesに読み込ませる、モデルのラグランジアンを記入したfrファイルを作成します。3節で述べた実スカラーシングレットモデルのラグランジアンのモデルファイル("singletDM.fr")を以下に示します。

M$ModelName = "singletDM";
(****************  Parameters *************)

M$Parameters = {
  \[Lambda]S == {
       ParameterType -> External,
       InteractionOrder -> {QED, 2},
       ParameterName -> lS,
       Value -> 0.5,
       BlockName -> INVSCALAR},
  \[Lambda]hS == {
       ParameterType -> External,
       InteractionOrder -> {QED, 2},
       ParameterName -> lhs,
       Value -> 0.5,
       BlockName -> INVSCALAR},
   muS=={
     ParameterType -> Internal,
     Value     -> Sqrt[2 MDM^2-vev^2 \[Lambda]hS]/Sqrt[2],
     TeX       -> \[Mu]S,
     Description -> "Coefficient of the quadratic piece of the singlet potential"
} 
}

(********* Particle Classes **********)
M$ClassesDescription = {
   (* The new scalar sector *)
  S[4] == {
        ClassName -> SDM,
        SelfConjugate -> True,
        Mass -> {MDM, 400},
        Width -> 0,
        PropagatorLabel -> "~S",
        PropagatorType -> D,
        PropagatorArrow -> None,
        PDG             -> 990025,
        ParticleName    -> "~S",
        FullName        -> "~S"
}
}

(*****************************************************************************************)

(* The lagrangian for the new sector *)

LScalar := Block[{ii, mu, feynmangaugerules},
feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];
ExpandIndices[1/2 del[SDM, mu]^2 - 1/2 muS^2 SDM^2 - \[Lambda]S/(24) SDM^4 - \[Lambda]hS/(2) SDM^2 Phibar[ii] Phi[ii], FlavorExpand -> {SU2D, SU2W}/.feynmangaugerules]
];

　このモデルファイルについて解説します。まず1行目の"ModelName="singletDM";"はモデルの名前を記入します。自分がわかりやすいように好きな名前をつけて下さい。その次の"Parameters"ではモデルのパラメータを指定します。このモデルのインプットパラメータは($m_S, \lambda_{H S}, \lambda_S$)としたため$\lambda_S$、$\lambda_{HS}$は"ParameterType　-> External"となっています。また、$\mu_S$については"ParameterType　-> Internal"となっており、その値"Value"は$m_S$と$\lambda_{HS}$で表されています。"Description"では自分が将来わかるようにパラメータの説明を書くことができます。次に"ClassesDescription"ではモデルで登場する粒子を定義します。$S$は自身が反粒子となるため"SelfConjugate -> True"となります。また$S$は$Z_2$により安定な粒子となるため、崩壊せず" Width -> 0"となります(もし粒子が崩壊する場合でも崩壊幅は0として大丈夫です。MicrOMEGAsの方でWidthを計算してくれます)。粒子の名前は"~S"として前に~をつけて下さい。MicrOMEGAsの方では~がついている粒子を暗黒物質として認識します。最後にラグランジアンを定義します。ラグランジアンでは"1/2 del[SDM, mu]^2"は運動項、"1/2 muS^2 SDM^2"は質量項、"\ [Lambda]S/(24) SDM^4"は$S^4$の項、"\ [Lambda]hS/(2) SDM^2 Phibar[ii] Phi[ii]"は$S^2 |H|^2$の項を表しています。ここで"ii"はSU(2)の添え字、barはダガーの意味を表しています。"Block[{ii, mu, feynmangaugerules},,,]"はダミーの添え字を指定しています。"ExpandIndices"ではSU(2)の場をcomponent field で書くように指示します。 "feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];"ではUnitaryゲージを選べるようにしています。全体を通して、コンマや{}を忘れないように注意して下さい。また半角スペースは掛け算を表します("*"で掛け算を表すこともできるのでこちらの方が打ち間違いは少ないかもしれません)。作成した"singletDM.fr"ファイルはfeynrule-currentディレクトリ内にある"Models"の中に置きましょう。

　標準模型のラグランジアンファイル"SM.fr"も必要ですが、これはダウンロードしたfeynrules-currentフォルダの中に用意されてあります(FeynRuleのモデルデータベースにてダウンロードすることもできます。)。
　次にこのモデルファイルからCalchepファイルをmathematicaにて作成します。そして、feynrules-currentのディレクトリに以下のように"FR.nb"を作成します。このノートの"LoadModel"では標準模型のモデルファイルと実スカラーシングレットモデルのラグランジアンを読み込みます。"FeynmanGauge=False"とすることで、Unitaryゲージを指定することができます。"WriteCHOutput"ではCalchepファイルを生成します。その中のLSM+LScalarが標準模型のラグランジアンと実スカラーシングレットモデルのラグランジアンを足したもの、"Output"では出力されるフォルダの名前を指定するものです。

これを全て実行すると、feynrules-currentディレクトリ内に"singletDM"という名前のディレクトリが生成されます。この中にはCalchepファイルが格納されています。もしエラーが発生したらfrファイルを修正してQuit[]から全てやり直して下さい。これでMimrOMEGAsに読み込ませる用のモデルファイルを生成することができました。これを用いて、MicrOMEGAsでこのモデルから様々な物理量を計算することができます。
　ちなみ、他のモデルのファイルを作成する際は、FeynRuleのモデルデータベースからなるべく近いモデルのファイルをダウンロードして、マニュアルを読みつつ書き換えて行くようにして下さい。
#4. MicrOMEGAs
　まず、MicrOMEGAsのサイトからDownload and Installへ行きダウンロードします。解凍し、以下のようにmicromegasのディレクトリーでmakeします。

 micromegas_6.0 % make

micromegas_6.0のディレクトリにはモデルの計算を行うためのディレクトリを作成することができます。コマンドは以下のようになります。

micromegas_6.0 %　./newProject realscalar

今回はディレクトリの名前をrealscalarとしました。新しく作られたディレクトリにはモデルファイルを入れる必要があります。modelsのディレクトリ(パスは"micromegas_6.0/hogehoge/work/models")には4節で作成したモデルファイルを入れてください。またMicrOMEGAsで物理量を計算する際に逐一、崩壊率を計算させるようにしたいので、以下に示すように、モデルファイルのうち"prtcls1.mdl"のwidthの前に!をつけて下さい。

　MicrOMEGAsで計算するには、物理量を計算する関数が書き込まれたmain.cファイル、モデルのパラメータの数値が書き込まれたテキストファイル(data.txt)が必要です。これらファイルの書き方について以下で解説します。
main.cには物理量を計算する関数を書き込みます。今回は暗黒物質の残存量を計算するdarkOmega(&Xf,fast,Beps,&err)という関数を用います。これはモデルのパラメータから暗黒物質の残存量$\Omega h^2$を計算して、返す関数です。これを持いてOMEGAのモジュールの部分(#ifdef OMEGA)を以下のように書きます。

#ifdef OMEGA
{ int fast=-1;
  double Beps=1.E-4, cut=0.01;
  double Omega;  
  int i,err; 

 double Xf;
     Omega=darkOmega(&Xf,fast,Beps,&err);
     printf("%.2e",Omega);
}
#endif

その他のモジュールは使わないのでコメントアウトしましょう(main.cの冒頭#define hogehogeをコメントアウト)。また余分なprint分もコメントアウトしておきましょう。また、今回はUnitaryゲージを用いるので、main.cの最初の方にある"ForceUG"を"ForceUG=1"として下さい。
　次にdata.txtにはモデルのパラメータを数値として書き込みます。具体的には以下のように書き込みます。

MDM 50            
lhs 0.0740999999999603          
lS 0

MDMは暗黒物質の質量、lhs、lSは結合定数$\lambda_{HS}$、$\lambda_{S}$を表しています。パラメータとその数値との間は半角スペースを開けましょう。

最後にdata.txtをmain.cファイルに読み込ませることで計算結果がコマンドプロンプトに出力されます。具体的には、以下のようなコマンドをコマンドプロンプトに入力します。

realscalar %　./main data.txt

出力結果は以下のようになります。

realscalar %　1.20e-01

となりました。そのため、data.txtにて設定したパラメータの値で、残存量の観測値$\Omega h^2=0.12$が説明できることがわかります。
次に、$\Omega h^2=0.12$となるMDM、lhsの値を求めてグラフに表します。MDMを50から10000 GeVまで変化させ、それぞれのMDMに対して$\Omega h^2=0.12$となるlhsを決定します。これを手作業でやると大変なのでpythonスクリプト(名称をrelic.py)を用いて計算します。その一例として私が用いたコードを以下に示します。

import linecache
import subprocess
import pandas as pd
from subprocess import PIPE

#omeghの初期値を設定(どんな値でもよい)これは後にdarkOmegaで得た結果を判定するために用いる
omegh=1

path = "data.txt"
path1="result.txt"

#data.txtからパラメータの値を読み込む
colname=range(1,8,1)
df = pd.read_csv(path,sep=" ",names=colname)
MX=df.iloc[0,1].astype(int)
lhs=df.iloc[1,1]

#MXを50から10000まで変化させる
for MX in range(50,10001,1):
 while(omegh>0.1201 or omegh<0.1199):
  df = pd.read_csv(path,sep=" ",names=colname)
  df.iloc[0,1]=MX
  lhs=df.iloc[1,1]
  #print(MX)
  #print(lhs)

  #subprocessを用いてmain.cにdata.txtを読み込ませた際の結果を読み込む
  rest=subprocess.run(["./main data.txt"],shell=True,stdout=PIPE,stderr=PIPE, text=True)
  
  #出力結果の整形
  res=rest.stdout
  res=res.split(" ")
  res=pd.DataFrame(res)
  res=res.astype(float)
  #print(res)

  #計算の結果(パラメータ、Omegah2の値)をまとめる
  res=pd.DataFrame(res)
  res[1]=MX
  res[2]=lhs
  print(res)
  resultcolname=range(0,3,1)


  
  #omeghが0.12より大きいか小さいかによって、lhsの値を更新
  omegh=res.iloc[0,0]
  #print(omegh)
  if(omegh<0.12):
   df.iloc[1,1]=lhs-0.0001
   df.to_csv("data.txt",sep=" ",index=False,header=False)
  else:
   df.iloc[1,1]=lhs+0.0001
   df.to_csv("data.txt",sep=" ",index=False,header=False)
  
 #result.txtという結果をまとめるtxtファイルを読み込む
 result=pd.read_csv(path1,sep=" ",names=resultcolname)
 #Omegah2=0.12となるパラメータの値をresultにまとめる
 result=pd.concat([result,res],ignore_index=True,axis=0)
 result.to_csv("result.txt",sep=" ",index=False,header=False)
 print(result)
 #omeghの初期化
 omegh=1

このプログラムにおいて、result.txtは$\Omega h^2=0.12$となるパラメータの値を書き込むファイルであり、このプログラムを動かす前にrealscalarのディレクトリ内に作成しておきます。このコードはrelic.pyとして保存しておきます。このスクリプトファイルを以下のコマンドで実行します。

realscalar %　python3 relic.py

そうすることで、result.txtに$\Omega h^2=0.12$となるMDM、lhsの値が書き込まれます。$\Omega h^2$の値が同じところを行ったり来たりしてループしてしまう場合はlhsの変化のはばを小さくするか、break文を追加してループを抜け出すように変更すると良いです。この結果をグラフにすると以下のようになります。

　次に求めたパラメータに対して、直接探索実験の散乱断面積を計算します。まず、main.cのCDM_NUCLEONのモジュールを以下のように書き換えます。

#ifdef CDM_NUCLEON
{ double pA0[2],pA5[2],nA0[2],nA5[2];
 double Nmass=0.939; /*nucleon mass*/
 double SCcoeff;
 double csSIp;
 int sI,sD;

 nucleonAmplitudes(CDM[1], pA0,pA5,nA0,nA5);
 SCcoeff=4/M_PI*3.8937966E8*pow(Nmass*Mcdm/(Nmass+ Mcdm),2.);
 csSIp=  SCcoeff*pA0[0]*pA0[0];
 printf("%.3e",csSIp);
}
#endif

その他のモジュールは使わないのでコメントアウトしましょう(main.cの冒頭#define hogehogeをコメントアウト)。また余分なprint分もコメントアウトしておきましょう。また、このモジュールを用いた散乱断面積の計算結果の単位はpbとなっていることに注意してください。次に、以下のようなプログラムを考えます。まずresult.txtから1行分を取り出して、data.txtを書き換えます。次に、そのdata.txtをmain.cに読み込ませることで散乱断面積を計算します。そしてその計算結果をresult_cs.txtに書き込みます。これをresult.txtの全ての行に対して行い、残存量を説明できるパラメータで直接探索実験での散乱断面積を計算します。これを手作業でやると大変なのでpythonスクリプト(名称をcs.py)を用いて計算します。その一例として私が用いたコードを以下に示します。

import linecache
import subprocess
import pandas as pd
from subprocess import PIPE

#data.txt, result.txtを読み込ませる
path = "data.txt"
path1="result.txt"
colname=range(1,15,1)
df = pd.read_csv(path,sep=" ",names=colname)
rs = pd.read_csv(path1,sep=" ",names=colname)

#result.txtの1行分を取り出し、data.txtに書き込む
for i in range(0,10001,1):
 vs=rs.iloc[i,2]
 MX=rs.iloc[i,1]
 df.iloc[1,1]=vs
 df.iloc[0,1]=MX
 df.to_csv("data.txt",sep=" ",index=False,header=False)

 #data.txtの内容をmain.cに読み込ませる
 rest=subprocess.run(["./main data.txt"],shell=True,stdout=PIPE,stderr=PIPE,  text=True)

 #出力結果の整形
 res=rest.stdout
 res=res.split(" ")
 print(res)
 res=pd.DataFrame(res)
 res=res.astype(float)
 res[1]=MX
 
 #計算結果pbをcm^2に変換
 res.iloc[0,0]=res.iloc[0,0]*pow(10, -36)
 print(res)

 #出力結果の書き込み
 resultcolname=range(0,2,1)
 path2="result_cs.txt"
 result=pd.read_csv(path2,sep=" ",names=resultcolname)
 result=pd.concat([result,res],ignore_index=True,axis=0)
 print(result)
 result.to_csv("result_cs.txt",sep=" ",index=False,header=False)

このスクリプトファイルを以下のコマンドで実行します。

realscalar %　python3 cs.py

そうすることで、result_cs.txtに$\Omega h^2=0.12$となるMDM、lhsの値で計算した散乱断面積が書き込まれます。この結果をグラフにすると以下のようになります。グラフの黒線が計算した散乱断面積です。他の実線は直接探索実験により与えられた断面積の上限値、点線は将来観測される範囲を表しています。グラフより, $m_S ≃ 62.5$ GeV(= $m_h /2$)の領域と, 4000 GeV $≲$ $m_S$の領域以外は理論計算が実験で与えられている散乱断面積の上限を上回ってしまうことが分かります。そのため, 直接探索実験の結果からそのパラメータ領域は除外されることがわかります。このように, 実スカラーシングレットモデルの様なシンプルなモデルではダークマターの量を説明し, かつ直接探索実験の結果を説明することが困難となってきています。そのため、この2つの実験結果を自然に説明できるモデルを作る研究が盛んに行われています。

　以上、MicrOMEGAsを用いて暗黒物質のモデルから様々な物理量を計算し実験結果と比較できる方法を紹介しました。その他の計算機能はぜひ、MicrOMEGAsのマニュアルを参照してみてください。また、"micromegas_6.0/STFM"のディレクトリにはdemo.cファイルがあり、main.cファイルを書くときの参考になると思います。他にもmicromegas_6.0のディレクトリ内にある他のモデルのディレクトリにもmain.cがあるので、そちらも活用して自分が計算したい物理量に合わせてmain.cを書き換えてみてください。