암흑물질 단면적 계산-micromegas

 

MicrOMEGAs 프로그램이 설치되어 있어야 합니다. 계산은 wsl 환경에서 수행, 시각화는 python 으로 했습니다.

 

 

import os
import subprocess
import re
import numpy as np
import matplotlib.pyplot as plt

# MicrOMEGAs 실행파일 및 .par 경로
MAIN = os.path.expanduser("~/micromegas_6.2.3/SingletDM/main")
PAR  = os.path.expanduser("~/micromegas_6.2.3/SingletDM/SingletDM.par")

# -----------------------------------
# .par 파일 작성
# -----------------------------------
def write_singletdm_par(path, Q, Mh, laS, laSH, Mdm1):
    with open(path, "w") as f:
        f.write(f"Q {Q}\n")
        f.write(f"Mh {Mh}\n")
        f.write(f"laS {laS}\n")
        f.write(f"laSH {laSH}\n")
        f.write(f"Mdm1 {Mdm1}\n")

# -----------------------------------
# MicrOMEGAs 실행
# -----------------------------------
def run_micromegas(main_path, parfile_path):
    result = subprocess.run([main_path, parfile_path], capture_output=True, text=True)
    return result.stdout

# -----------------------------------
# σv 파싱 (relic density + 채널별 기여)
# -----------------------------------
def parse_channels_and_sigmaV(text):
    channels = {}
    lines = text.splitlines()
    capture = False
    for line in lines:
        line = line.strip()
        if line.startswith("# Channels which contribute"):
            capture = True
            continue
        if capture:
            if line == "" or line.startswith("===="):
                break
            m = re.match(r'(\d+)% .*->\s*(\S+)', line)
            if m:
                frac = float(m.group(1))/100
                ch = m.group(2)
                channels[ch] = frac
    m = re.search(r'Omega=([0-9.Ee+-]+)', text)
    if m:
        Omega = float(m.group(1))
        sigmaV_tot = 3e-26 * 0.12 / Omega  # cm^3/s 근사
    else:
        sigmaV_tot = None
    sigmaV_channels = {k: sigmaV_tot*v for k,v in channels.items()} if sigmaV_tot else None
    return sigmaV_tot, sigmaV_channels

# -----------------------------------
# DM mass scan
# -----------------------------------
masses = np.linspace(10, 1000, 100)  # 10~1000 GeV
sigmaV_totals = []



channel_dict = {}

for Mdm in masses:
    write_singletdm_par(PAR, 14, 125, 0.1, 0.15, Mdm)
    output = run_micromegas(MAIN, PAR)
    sigmaV_tot, sigmaV_channels = parse_channels_and_sigmaV(output)
    sigmaV_totals.append(sigmaV_tot)
    
    if sigmaV_channels:
        # 새로운 채널이 있으면 리스트 초기화
        for ch in sigmaV_channels:
            if ch not in channel_dict:
                # 이전 mass들에 None으로 채워 길이 맞춤
                channel_dict[ch] = [None]*(len(sigmaV_totals)-1)
        # 각 채널별 값 추가
        for ch in channel_dict:
            channel_dict[ch].append(sigmaV_channels.get(ch, None))
    else:
        # sigmaV_channels가 None이면 모든 채널 None 추가
        for ch in channel_dict:
            channel_dict[ch].append(None)

# -----------------------------------
# 그래프 그리기
# -----------------------------------
plt.figure(figsize=(8,6))
plt.plot(masses, sigmaV_totals, marker='o', label='Total σv', color='black')

# 주요 채널별 σv (선택적으로 표시)
colors = ['red', 'blue', 'green', 'orange', 'purple']
for i, (ch, vals) in enumerate(channel_dict.items()):
    plt.plot(masses, vals, marker='x', linestyle='--', color=colors[i%len(colors)], label=f'{ch} channel')

plt.yscale('log')
plt.xlabel("DM Mass Mdm [GeV]")
plt.ylabel("Annihilation cross section σv [cm³/s]")
plt.title("SingletDM: σv vs DM Mass")
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.legend()
plt.show()