引力波数据分析第三章题目

本系列文章旨在充实博客内容记录引力波数据分析课程的作业成果,标记为引用的题目部分是直接复制自胡教授引力波数据处理讲义的内容,大部分公式也来自该讲义,不再一一标出。

如果真的有人搜索到了本系列文章并参考,最好能标记一下引用哦。

Exercise 3.1

考虑质量为 的双黑洞系统,假设其距离我们 ,倾角 。请编写程序作图,

时间内频率随时间变化的图。

时间内相位随时间变化的图。

时间域波形(假设 )。

对应的频率域波形,考虑 0PN(以 Bode 图形式给出)

考虑质量为 的双黑洞系统,距离为 ,倾角 ,画图:

频率随时间变化

在 Newtonian 近似下,我们可以得到旋近阶段引力波频率的演化公式:

其中 为啁啾质量,积分后易得频率与时间的关系为:

那么我们很容易画出在 时间内频率随时间变化的关系,如图所示,具体代码附在最后。

频率随时间变化图

相位随时间变化

根据相位随时间变化的公式:

那么我们很容易画出在 时间内相位随时间变化的关系,令初始相位 ,如图所示。

相位随时间变化图

时间域波形

假定 ,那么

时间域信号为:

其中

那么我们很容易画出在 时间内时间域的信号,令初始相位 ,如图所示。

时间域波形

频率域波形

对于正频率部分,相位项为

那么我们很容易画出在 时间内频率域波形,令初始相位 ,如图所示。

频率域波形

最后给出绘图使用的代码

Exercise 3.2

设中子星表面磁场强度 ,有 ,其中 是磁偶极矩, 是中子星半径, 是转动惯量。

请推导出 的表达式(请用 表示)

若该中子星有观测量 ,假定中子星辐射完全由磁偶极矩辐射主导,请计算出其年龄

若该中子星有观测量 ,假定中子星辐射完全由磁偶极矩辐射主导,请计算出其磁场强度

若该中子星有观测量 ,假定中子星辐射完全由磁偶极矩辐射主导,请计算出其辐射功率

假定中子星辐射完全由引力波辐射主导,请根据公式3.29,计算出上述三种情况下的

是中子星半径, 是中子星转动惯量

B的表达式

中子星系统动能的减慢速率

同时磁偶极矩辐射的功率为

其中磁偶极矩 ,可以得到

年龄

若中子星辐射完全由磁偶极辐射主导,有 ,那么年龄

磁场强度

若中子星辐射完全由磁偶极辐射主导,有 ,那么磁场强度

辐射功率

若中子星辐射完全由磁偶极辐射主导,有 ,那么辐射功率

引力波辐射

若中子星辐射完全由引力波辐射主导,有

可以得到

,那么

,那么

,那么

Exercise 3.3

请在 平面上作图, 的取值范围为 的取值范围为 ,以对数展示,画出

的等年龄线(假定中子星辐射完全由磁偶极矩辐射主导)

的等磁场强度线

的等辐射强度线

参考了psrqpy的画图配置 psrqpy: a python interface for querying the ATNF pulsar catalogue,如图所示,绘图代码在最后给出。

$P-\dot{P}$图

Exercise 3.4

考虑 ,假设磁场强度 、椭率 、转动惯量 等均为常数。记 ,试证明 。已知 PSR J1640-4631 的制动指数 ,计算相应的

已知

可以得到

已知 PSR J1640-4631 的制动指数 ,相应的

绘图代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as cons

# constants
# sun mass
sun = 1.98855e30
# chirp mass
chirp = np.power(5*sun,3/5)*np.power(20*sun,2/5)
distance = 1e8*cons.parsec

def time(t, n):
return np.linspace(-t,0,n)

def frequency(t):
return np.power(5,3/8)/(8*cons.pi)*np.power(np.power(cons.c,3)/(cons.G*chirp),5/8)*np.power((-t),-3/8)

plt.figure(figsize = (6, 6))
plt.plot(time(cons.year, 10000), frequency(time(cons.year, 10000)))
plt.xlabel('Time(s)')
plt.ylabel('Frequency(Hz)')
plt.title('Inspiral', fontsize = 14)
plt.savefig('frequency.jpg')
plt.show()

def phi(t):
return -2*np.power(np.power(cons.c,3)/(5*cons.G*chirp)*(-t),5/8)

plt.figure(figsize = (6, 6))
plt.plot(time(cons.year, 10000), phi(time(cons.year, 10000)))
plt.xlabel('Time(s)')
plt.ylabel(r'$\Phi$')
plt.title('Inspiral', fontsize = 14)
plt.savefig('phi.jpg')
plt.show()

def radius(t):
return np.power(256*np.power(cons.G,3)/(5*np.power(cons.c,5))*5*sun*np.power(10*sun,2),1/4)*np.power(-t,1/4)

def h(t):
omega = np.sqrt(cons.G*20*sun/(np.power(radius(t),3)))
return cons.G*5*sun*np.power(radius(t),2)*np.power(omega,2)/(np.power(cons.c,4)*distance)*np.cos(phi(t))

plt.figure(figsize = (6, 6))
plt.plot(time(1, 10000), h(time(1, 10000)))
plt.xlabel('Time(s)')
plt.ylabel('h')
plt.title('Inspiral', fontsize = 14)
plt.savefig('ht.jpg')
plt.show()

def hf(f):
phif = -cons.pi/4+3/4*np.power(8*cons.pi*cons.G*chirp*f/np.power(cons.c,3),-5/3)
return np.sqrt(5/384)/distance*np.power(cons.G,5/6)*np.power(cons.c,-3/2)*np.power(chirp,5/6)*np.power(np.pi,-2/3)*np.power(f,-7/6)*np.cos(phif)

freq = np.linspace(np.log(frequency(-cons.year)), np.log(frequency(-1e-3)), 5000)
freq = np.exp(freq)

plt.figure(figsize = (6, 6))
plt.loglog(freq, hf(freq))
plt.xlabel('Frequency(Hz)')
plt.ylabel('h')
plt.title('Inspiral', fontsize = 14)
plt.savefig('hf.jpg')
plt.show()

# G in CGS
G = 6.67428e-8

def freq_dot_tau(freq, tau):
return(freq/(2*tau*cons.year))

def freq_dot_gauss(freq, gauss):
return(np.power(1e18*gauss,2)*8*cons.pi*cons.pi/(3*np.power(cons.c*100,3)*1e45*freq))

def freq_dot_ergs(freq, ergs):
return(ergs*np.power(freq,3)/(4*cons.pi*cons.pi*1e45))

freq = np.power(10, np.linspace(-3,1,1000))

fig, ax = plt.subplots()
plt.ylim((1e-22,1e-8))
plt.xlim((1e-3,10))

for i in (1e3, 1e5, 1e7, 1e9):
plt.loglog(freq, freq_dot_tau(freq, i), 'k--', linewidth=1.5)
xytext = (0, 5)
text = ax.annotate(
f'$10^{{ {int(np.log10(i)):d} }}$ yr',
xy=(1e-3, np.min(freq_dot_tau(freq, i)*1.5)),
xytext=xytext,
textcoords="offset points",
size=15,
color='k',
zorder=1,
horizontalalignment="left",
verticalalignment="center_baseline",
)
text.set_rotation_mode("anchor")
text.set_rotation(10)

for i in (1e9, 1e10, 1e11, 1e12, 1e13, 1e14):
plt.loglog(freq, freq_dot_gauss(freq, i), 'r:', linewidth=1.5)
xytext = (0, 5)
text = ax.annotate(
f'$10^{{ {int(np.log10(i)):d} }}$ G',
xy=(6, np.min(freq_dot_gauss(freq, i)*2)),
xytext=xytext,
textcoords="offset points",
size=15,
color='r',
zorder=1,
horizontalalignment="left",
verticalalignment="center_baseline",
)
text.set_rotation_mode("anchor")
text.set_rotation(-10)

for i in (1e30, 1e33, 1e36):
plt.loglog(freq, freq_dot_ergs(freq, i), 'b-.', linewidth=1.5)
xytext = (0, 5)
text = ax.annotate(
f'$10^{{ {int(np.log10(i)):d} }}$ erg/s',
xy=(1, freq_dot_ergs(1.1, i)),
xytext=xytext,
textcoords="offset points",
size=15,
color='b',
zorder=1,
horizontalalignment="left",
verticalalignment="center_baseline",
)
text.set_rotation_mode("anchor")
text.set_rotation(25)

plt.xlabel('Period(s)')
plt.ylabel('Period Derivative')
plt.show()