第二章 理論架構
2.2 質量公式推演
我們可以更進一步的定義函數H(r,r’)： 然後把delta function表示成平面波的形式來展開：
∫
−∞∞∫
∞∂經過演算後，我們可以改寫成：
2.3.1 Famaey and Binney Model
在FB模型裡[3]，牛頓速度 和觀測到的速度 v 由(1)和(3)可得知有以
之後我們就可以推出速度的公式：
第三章 程式及數值方法
(Maxwellian)熱攪動情況下，譜線的輪廓主要將由熱攪動所引起的都普勒效應來 生的效應。這一輪廓的半寬高稱為「都普勒寬度」(Doppler width)：
ν τ
觀測在星系中央的H-alpha和[NII] line的放射物得出CCD的光譜。最後利用 envelop-tracing方法來從PV圖得到旋轉曲線[13]。所以我們只要觀察星系分式所給出： 絕對星等(Absolute Magnitude)：
「視星等」是對在地球上所接收到的流量作的測量。它與我們與天體的距離有
(43)
VN without gas
V Vgas
Figure 1 rotation curve and Newtonian velocity and gas velocity
在這裡橫軸使用的單位為 KPC，而縱軸為 2
2
s
Km
。(本論文後有詳細的單位轉換) 粗黑線為旋轉曲線，虛線逗點線為 Newtonian 速度曲線，而短虛線為不含氣體的Newtonian 速度，最後長虛線為氣體的速度。而旋轉曲線的誤差範圍大概在 +4.1%~ -5.1%之間。旋轉曲線的取得是來自於：
http://www.ioa.s.u-tokyo.ac.jp/~sofue/RC99/2403.dat
VelocitySquare:
Figure 2 Square of velocity and fitting curve and delta velocity
虛線為擬合的曲線，擬合的方程式為(32)，參數為 C0= a1=7.2 KPC b1= 0.1 KPC 實線為數據資料由(42)所算出的曲線，而細黑線為兩者之相 差的數據曲線縱軸單位為公里的平方除以秒平方，橫軸單位為 1000 乘以秒差距。
10 27
7 . 7 × −
2
)2
(
pc s
Mass density
Figure3 mass density and curve fitting
橫軸依然是 KPC，為了等下我們處理資料方便，我們把縱軸不使用 M.K.S 制而 使用了太陽質量除以 PC 平方。
而相關的參數為：C0=7.7×10−27 a1=7.2 KPC b1= 0.1 KPC 2 )2
(
s pc
實線為數據資料所算出的質量密度(40)，而虛線為擬合質量密度(31)。
M/L
Figure 4 mass density over luminosity 縱軸為太陽質量除以太陽光度，橫軸為 1000 秒差距
3.3.2 NGC2903
DIAM1 : 21.9 DIAM2 : 12.3 D=3.25(MPC)
Figure 5 rotation curve and Newtonian velocity 在這裡橫軸使用的單位為 KPC，而縱軸為 2
2
s
Km
。粗黑線為旋轉曲線，虛線逗點 線為 Newtonian 速度曲線，而短虛線為不含氣體的 Newtonian 速度，最後長虛線 為氣體的速度。誤差範圍在+8%~ -18%之間。旋轉曲線的取得是來自於：http://www.ioa.s.u-tokyo.ac.jp/~sofue/RC99/2903.dat
VelocitySquare:
Figure 6 Square of velocity and fitting curve and delta velocity
虛線為擬合的曲線，擬合的方程式為(32)，參數為 C0= a1=5 KPC b1= 0.14 KPC 實線為數據資料由(42)所算出的曲線，而細黑線為兩者之 相差數據曲線縱軸單位為公里的平方除以秒平方，橫軸單位為 1000 乘以秒差距。
10 26
4 . 4 × −
2
)2
(
pc
s
Mass density
Figure 7 mass density and curve fitting
橫軸依然是 KPC，為了等下我們處理資料方便，我們把縱軸不使用 M.K.S 制而 使用了太陽質量除以 PC 平方。實線為數據資料所算出的質量密度(40)，而虛線 為擬合質量密度(31)。相關的參數為：C0= a1=5 KPC b1=
0.14 KPC
10 26
4 .
4 × − 2
)
( 2
s pc
M/L
Figure 8 mass density over luminosity
縱軸為太陽質量除以太陽光度，橫軸為
1000秒差距
3.3.3 NGC 6503
DIAM1 : 7.1 DIAM2 : 2.4 D=5.9(MPC)
0 5 10 15 20 25
0 20 40 60 80 100 120
kpc
km/s
NGC6503
VN
VN without gas
V Vgas
Figure 9 rotation curve and Newtonian velocity
在這裡橫軸使用的單位為
KPC，而縱軸為
22
s
Km 。粗黑線為旋轉曲線，虛線
逗點線為 Newtonian 速度曲線，而短虛線為不含氣體的 Newtonian 速度，最後長
虛線為氣體的速度。縱軸單位為公里的平方除以秒平方，橫軸單位為
1000乘以
秒差距。旋轉曲線的取得是來自於[12]。
VelocitySquare
Figure 10 Square of velocity and fitting curve and delta velocity
虛線為擬合的曲線，擬合的方程式為(32)，參數為 C0=
C01=
10 26
24 .
4 × − 10 27
76 .
8 × −
2
)2
(
pc s
a1=1.2 KPC ， a2= 4.6 KPC
2
)2
(
s pc
b1=0.9 KPC， b2=1 KPC
實線為數據資料由(42)所算出的曲線，而細黑線為兩者之相差數據曲線，縱軸單 位為公里的平方除以秒平方，橫軸單位為 1000 乘以秒差距。
Mass density
Figure 11 mass density and curve fitting 相關的參數為：C0=4.24×10−26 C01=8.76×10−27
a1=1.2 KPC ， a2= 4.6 KPC b1=0.9 KPC， b2=1 KPC
2
)
( 2
s pc
2
)2
( s pc
橫軸依然是 KPC，為了等下我們處理資料方便，我們把縱軸不使用 M.K.S 制而 使用了太陽質量除以 PC 平方。實線為數據資料所算出的質量密度(40)，而虛線 為擬合質量密度(31)。
M/L
Figure 12 mass density over luminosity 縱軸為太陽質量除以太陽光度，橫軸為 1000 秒差距
第四章 結論
在前面的部分我們討論了螺旋星系的動力學，得到了旋轉曲線的解析解，
也利用數值運算擬合出了解析解的可靠性。現在我們想更進一步的討論亮度與質 量之間的關係。無蹤質量一直是天文學的一個問題，暗物質在目前看來是提供這 個問題的解答之一，但是經由 MOND 的理論我們可以得知，修改重力定律或許 也是這個問題的解答。在前面的章節，我們經由 MOND 理論算出了發光物質所 佔星系的質量，接下來我們將要比較發光物質的質量與亮度的比值。
亮度：
NGC 2403：
Figure 13 Luminosity of NGC 2403 縱軸為亮度橫軸為 1000 秒差距
NGC 2903
Figure 14 Luminosity of NGC 2903 縱軸為亮度橫軸為 1000 秒差距
NGC 6503
Figure 15 Luminosity of NGC 6503 縱軸為亮度橫軸為 1000 秒差距，以上的亮度資料來自[6],[7]。
我們將這些亮度資料分別與其對應的發光物質質量相比較，得到了下圖
Figure 16 M/L of NGC2403 NGC2903 NGC6503 一般天文學裡認定的 M/L 值如下：
NGC 2403 1.4 NGC 2903 3.5 NGC 6503 0.041
橫軸的單位為 KPC，縱軸的單位為「太陽質量除以 L 的太陽光度」，在天文學理 我們定義太陽的 M/L=1。在一般的天文學裡，我們知道質量是與 r 成反比，而光
度是與 exp(-r)成比例，所以當 r 很大時，M/L 必然增加得很強烈。這是推斷無蹤 質量的一個關鍵。但當我們使用 MOND 理論來計算 M/L 時，可以得到較平緩的 數值，這顯示了可能並非是有其他看不見的物質在星系邊緣，而是有可能在大尺 度下重力定律必須被修正，就如同本論文的方法。到目前為止，還沒有一個方法 能夠有效的說明無蹤質量這個問題，但 MOND 理論似乎提供了我們一條新的 路，可以去嘗試解決這個問題。而且透過了這個理論，我們更進一步瞭解了螺旋 星系的動力學，這樣對於我們要去探討天文學的問題時，我想是提供了一些小小 的幫助。而這個理論正在發展中，也在一些星系理得到了成功的吻合，真理是要 棄而不捨得努力去追尋，我想 MOND 理論提供了我們更多的線索，讓我們可以 一步步逼向最後的真理。
參考資料：
[1]Arfken, G. B. and Weber, H. J., 2005, Mathematical method for physicists, [2]B. Famaey and J. Binney, “Modified Newtonian dynamics in the Milky Way”,
astro-ph/0506723
[3] Kao,W.F.,’’ Modified Newtonian Dynamics In Dimensionless Form’’, astro-ph/0504009;
[4] Kao, W.F., “Modified Newtonian dynamics and induced gravity”, 2005, gr-qc/0512062;
[5] Kao, W.F. “MOND as an extra dimensional effect”, astro-ph/0603311;
[6] Kent, S. M., “Dark matter in spiral galaxies. I. galaxies with optical rotation curves”, 1987, Astron. J. 93,816k
[7] Kent, S. M., “Dark matter in spiral galaxies. I. galaxies with optical rotation curves”, 1986, Astron. J. 91, 1301
[8] Mak,M. K. and Harko, T., “Can the galactic rotation curves be explained in brane world models?” Phys.Rev. D70 (2004)024010;
[9] Milgrom, M., “A modification of the newtonian dynamics: Implications for galaxies”, 1983, ApJ, 270, 371;
[10]”The physical universe-an introduction to Astronomy” Frank Shu [11]”Galaxies and Cosmology” F.Combes P.Boissé A.Mazure A.Blanchard [12]J. R. Brownstein1 and J. W. Moffat2
,”
Galaxy Rotation Curves WithoutNon-Baryonic Dark Matter”,arXiv:astro-ph/0506370 v3 28 Jun 2005[13] Y.
Sofue, Y. Tutui, M. Honma, A. Tomita, T. Takamiya, J. Koda and Y.
Takeda ,”CENTRAL ROTATION CURVES OF SPIRAL
GALAXIES”(Astrophysical Journal Vol. 523, pp136-146, 1999)
[14] “觀測天文物理學＂ 孫維新，胡景耀 譯
APPENDIX Ⅰ
物理和天文常數單位換算
光速 c 299792458m/s 重力常數G 6.673 x 10-11 m3/kgs 秒差距pc 3.0856775807 x 10-16 m ＝3.261633 light yr 太陽質量 1.989×1030
kg
太陽亮度 3.826×1026W
名稱
λ
0[ ] µ m
∆λ
0[ ] µ m
0[
2 1]
−
−
m
Wm
e µ e
0[ ] Jy
U 0.36 0.068 4.35×10−8 1,880 紫外 B 0.44 0.098 7.20×10−8 4,650 藍 V 0.55 0.089 3.92×10−8 3,950 可見 R 0.70 0.22 1.76×10−8 2,870 紅
I 0.90 0.24 8.3×10−9 2,240 紅外 J 1.25 0.30 3.4×10−9 1,770 紅外 H 1.65 0.35 7×10−10 636 紅外 K 2.20 0.40 3.9×10−10 629 紅外 L 3.40 0.55 8.1×10−11 312 紅外 M 5.0 0.3 2.2×10−11 183 紅外 N 10.2 5 1.23×10−12 43 紅外 Q 21.0 8 6.8×10−14 10 紅外
Table 1 標準測光系統 1Jy=10−26
Wm
−2Hz
−1APPENDIX Ⅱ 數值計算程式
1. 計算星系質量密度(MATLAB)：
clear; clc; close all; tic;
D=9.1; band=4.83; %variable
KPC=3.09.*(10.^19); g0=1.2.*(10.^(-10)); G=6.67.*(10.^(-11));
PC=3.09.*(10.^16); SOLARMASS=1.99.*(10.^30); %constant
fid=fopen( '2903.txt','rt');
data1=fscanf(fid,'%f',[2 inf]);
fclose(fid);
data1=data1'; % Rotation curve total
fid=fopen( '2903gas.txt','rt');
data2=fscanf(fid,'%f',[2 inf]);
fclose(fid);
data2=data2';
data2=[0 0;data2]; %Rotation curve gas
fid=fopen( '2903L1.txt','rt');
data3=fscanf(fid,'%f',[2 inf]);
fclose(fid);
data3=data3'; %bright
fid=fopen( '2903L2.txt','rt');
data4=fscanf(fid,'%f',[2 inf]);
fclose(fid);
data4=data4'; %bright
R=20; VelocitySquare=[0];
while VelocitySquare(end)>=0 R=R+20;
r=0:0.05:R; r=r';
newdata1=zeros(length(r),1);
for n=1:length(r)
newdata2=zeros(length(r),1);
for n=1:length(r)
VelocitySquare=[0 0;r(2:n) VelocitySquare(1:n-1);zeropoint 0];
dyx=diff(VelocitySquare(:,2))./diff(VelocitySquare(:,1).*KPC);
midpoint=zeros(length(VelocitySquare)-1,1);
for n=1:length(midpoint)
midpoint(n)=(VelocitySquare(n,1)+VelocitySquare(n+1,1))./2;
end
diffVelocitySquare=[0 dyx(1);midpoint dyx;zeropoint 0];
M=0.0005:0.001:0.9995; W=1./M;
K=ellipke(M.^2);
i1=diffVelocitySquare(:,1)*M; i1=i1';
i2=diffVelocitySquare(:,1)*W; i2=i2';
M=M'; W=W'; K=K';
I1=zeros(length(M),length(diffVelocitySquare));
for n=1:length(M)
for m=1:length(diffVelocitySquare)
if i1(n,m)<=diffVelocitySquare(end,1)
I1(n,m)=interp1(diffVelocitySquare(:,1),diffVelocitySquare(:,2),i1(n,m));
end end end
K1=zeros(length(diffVelocitySquare),1);
for m=1:length(diffVelocitySquare) K1(m)=trapz(M,K.*I1(:,m));
end
I2=zeros(length(M),length(diffVelocitySquare));
for n=1:length(M)
for m=1:length(diffVelocitySquare)
if i2(n,m)<=diffVelocitySquare(end,1)
I2(n,m)=interp1(diffVelocitySquare(:,1),diffVelocitySquare(:,2),i2(n,m));
end end end
K2=zeros(length(diffVelocitySquare),1);
for m=1:length(diffVelocitySquare) K2(m)=trapz(M,K.*I2(:,m).*W);
end
density=(K1+K2)./(G.*(pi.^2));
n=1;
while density(n)>=0 n=n+1;
end
density=(PC.^2).*density(1:n-1)./SOLARMASS;
totalmass=zeros(length(density)-1,1);
for n=1:length(totalmass)
totalmass(n)=trapz(diffVelocitySquare(1:n+1,1),2.*pi.*diffVelocitySquare(1:n+1,1).*
density(1:n+1));
end
totalmass=[0;totalmass];
bright=[data3;data4(:,1).*60 data4(:,2)];
radius=bright(:,1).*pi.*D.*1000./(180.*3600);
Mag=bright(:,2)-(5.*log10(D.*100000));
luminosity=(2.512.^(band-Mag))./((D.*1000000.*pi./(180.*3600)).^2);
totalluminosity=zeros(length(luminosity)-1,1);
for n=1:length(totalluminosity)
totalluminosity(n)=trapz(radius(1:n+1),2.*pi.*radius(1:n+1).*luminosity(1:n+1));
end
totalluminosity=[0;totalluminosity];
if diffVelocitySquare(length(density),1)<=radius(end) last=diffVelocitySquare(length(density),1);
else
last=radius(end);
end
r=0:0.1:last; r=r';
newdensity=zeros(length(r),1);
for n=1:length(r)
newdensity(n)=interp1(diffVelocitySquare(1:length(density),1),density,r(n),'spline' );
end
newtotalmass=zeros(length(r),1);
for n=1:length(r)
newtotalmass(n)=interp1(diffVelocitySquare(1:length(density),1),totalmass,r(n),'spl ine');
end
newluminosity=zeros(length(r),1);
for n=1:length(r)
newluminosity(n)=interp1(radius,luminosity,r(n),'spline');
end
newtotalluminosity=zeros(length(r),1);
for n=1:length(r)
newtotalluminosity(n)=interp1(radius,totalluminosity,r(n),'spline');
end
ml=newdensity./newluminosity;
ML=newtotalmass(2:end)./newtotalluminosity(2:end);
figure(1);
plot(r,ml);
figure(2);
plot(r(2:end),ML);
toc;
2 擬合曲線程式(MATHEMATICA)：
s=ReadList["2903_12.txt",{Number,Number}]
<<Statistics`NonlinearFit`
NonlinearFit[s,theta1*((x^2+theta2^2)^0.5-(x^2+theta3^2)^0.5)/((x^2+theta2^2)^0.5
*(x^2+theta3^2)^0.5)+theta4*((x^2+theta5^2)^0.5-(x^2+theta6^2)^0.5)/((x^2+theta5
^2)^0.5*(x^2+theta6^2)^0.5),x,{theta1,theta2,theta3,theta4,theta5,theta6}]
BestFitParameters/.NonlinearRegress[s,theta1*((x^2+theta2^2)^0.5-(x^2+theta3^2)^0 .5)/((x^2+theta2^2)^0.5*(x^2+theta3^2)^0.5)+theta4*((x^2+theta5^2)^0.5-(x^2+theta 6^2)^0.5)/((x^2+theta5^2)^0.5*(x^2+theta6^2)^0.5),x,{theta1,theta2,theta3,theta4,the ta5,theta6},RegressionReport BestFitParameters]
f3=Block[{theta1=1303.4055606525544,theta2=-5.6503077814271405,theta3=-0.148 39929878772273,theta4=0.006716632178861634,theta5=-0.677464557717977,theta6
=-2.506588798706922*10^(-6)},Plot[theta1*((x^2+theta2^2)^0.5-(x^2+theta3^2)^0.5 )/((x^2+theta2^2)^0.5*(x^2+theta3^2)^0.5)+theta4*((x^2+theta5^2)^0.5-(x^2+theta6
^2)^0.5)/((x^2+theta5^2)^0.5*(x^2+theta6^2)^0.5),{x,0,15}]]
f4=ListPlot[s,PlotRange {0,20000},PlotJoined True,PlotStyle Hue[.6] ] Show[f3,f4]
3.擬合參數回推程式(MATLAB)：
newdata1=zeros(length(r),1);
for n=1:length(r)
newdata2=zeros(length(r),1);
for n=1:length(r)
VelocitySquare=[0 0;r(2:n) VelocitySquare(1:n-1);zeropoint 0];
theta1=0; theta2=28153.2; theta3=0; theta4=275; theta5=7; theta6=0.1121;
theta7=theta1.*2.*pi.*G2;theta8=theta4.*2.*pi.*G2
a2=theta2.*KPC; a3=theta3.*KPC; a5=theta5.*KPC; a6=theta6.*KPC
vn3=(theta7.*((theta2./((p.^2+theta2.^2).^0.5))-(theta3./((p.^2+theta3.^2).^0.5)))+
legend('VN','VN without gas','','V','Vgas',1)
title('NGC2403');
print -depsc 2403Velocity.eps;
figure(4);
plot(p,newvn2,'k--','LineWidth',2);
hold on
plot(p,vsq,'k','LineWidth',3);
hold on
plot(p,dv,'k');
legend('VN2f','VN2','delta VN2');
axis([0 25 -inf inf]);
xlabel('kpc');
ylabel('km^2/s^2');
print -depsc 2403Vn2.eps;
load dif24 figure(2);
plot(diffVelocitySquare(1:last,1),totaldensity(1:last),'b',diffVelocitySquare(1:las t,1),fitdensity(1:last),'r',diffVelocitySquare(1:last,1),density(1:last),'k');
legend('mu','mu_f') axis([0 10 -inf inf]) xlabel('kpc')
ylabel('solarmasss/(pc)^2');
print -depsc 2403mu.eps;
figure(3)
plot(r,ml,'k','LineWidth',3);
xlabel('kpc');
ylabel('solarmass/solarluminosity');
axis([0 8 0 inf]);
print -depsc 2403ml.eps;