[Python] What's making these two neutron stars not get closer to each other? (Kilonova...

import numpy as np
import matplotlib.pyplot as plt

G = 6.6743e-11
c = 299792458
R_ns = 12000
m1 = m2 = 3.5e31
init_dist = 1e4
dt = 1e-8

pos1 = np.array([-init_dist / 2, 0], dtype=float)
pos2 = np.array([init_dist / 2, 0], dtype=float)

# Orbital velocity magnitude for each star
v = np.sqrt(G * (m1 + m2) / init_dist)

# Initial velocities perpendicular to separation vector (x-axis)
v1 = np.array([0, v], dtype=float) # Star 1 moves +y
v2 = np.array([0, -v], dtype=float) # Star 2 moves -y


def compute_accelerations(pos1, pos2):
r_vec = pos2 - pos1
r_mag = np.linalg.norm(r_vec)
r_hat = r_vec / r_mag
force_mag = G * m1 * m2 / r_mag ** 2
a1 = force_mag / m1 * r_hat
a2 = -force_mag / m2 * r_hat
return a1, a2


a1, a2 = compute_accelerations(pos1, pos2)


def verlet_integration(pos1, pos2, a1, a2, v1, v2):
pos1_new = pos1 + v1 * dt + 0.5 * a1 * dt ** 2
pos2_new = pos2 + v2 * dt + 0.5 * a2 * dt ** 2

a1_new, a2_new = compute_accelerations(pos1_new, pos2_new)

v1_new = v1 + 0.5 * (a1 + a1_new) * dt
v2_new = v2 + 0.5 * (a2 + a2_new) * dt

return pos1_new, pos2_new, v1_new, v2_new, a1_new, a2_new


def peters_mathews(pos1, pos2, v1, v2):
r_vec = pos2 - pos1
r = np.linalg.norm(r_vec)
dr_dt = - (64 / 5) * (G ** 3 * m1 * m2 * (m1 + m2)) / (c ** 5 * r ** 3)
delta_r = dr_dt * dt
r_new = r + delta_r # dr_dt negative, so r_new < r

if r_new < R_ns / 2:
print("Merger condition reached. Exiting.")
exit()

r_hat = r_vec / r
pos1_new = pos1 + 0.5 * (r_new - r) * r_hat
pos2_new = pos2 - 0.5 * (r_new - r) * r_hat

# Calculate new orbital speed magnitude
v_mag = np.sqrt(G * (m1 + m2) / r_new)

# Calculate tangential direction by rotating r_hat by 90 degrees CCW
tangential_dir = np.array([-r_hat[1], r_hat[0]])

v1_new = v_mag * tangential_dir
v2_new = -v1_new # opposite direction for star 2

return pos1_new, pos2_new, v1_new, v2_new


# For plotting
pos1_list = []
pos2_list = []

for i in range(100000):
pos1, pos2, v1, v2, a1, a2 = verlet_integration(pos1, pos2, a1, a2, v1, v2)
if i % 100 == 0:
pos1, pos2, v1, v2 = peters_mathews(pos1, pos2, v1, v2)

pos1_list.append(pos1.copy())
pos2_list.append(pos2.copy())

pos1_arr = np.array(pos1_list)
pos2_arr = np.array(pos2_list)

plt.figure(figsize=(8, 8))
plt.plot(pos1_arr[:, 0], pos1_arr[:, 1], label="Star 1")
plt.plot(pos2_arr[:, 0], pos2_arr[:, 1], label="Star 2")
plt.scatter([0], [0], color='black', marker='x', label="Center of Mass")
plt.xlabel("x position (m)")
plt.ylabel("y position (m)")
plt.title("Binary system orbit simulation with Peters-Mathews decay")
plt.legend()
plt.axis('equal')
plt.grid()
plt.show()


Here's my code that I've made to simulate two neutron stars resulting in a kilonova. I want these two neutron stars to grow closer to each other as a result of the peters_mathews function which accounts for orbital decay. Additionally, using velocity verlet I calculated their orbital position. However, instead of the orbit decaying and them getting closer, they grow further apart. How should I fix this?

I've tried changing m1, m2, and initial distance by becoming larger or smaller. Also, the acceleration and velocity have been modified to bring them closer together. I have also implemented artificial decay which doesn't seem to work either. I want them to grow closer, not further apart.

I'm really frustrated with this program as I've been working on it for a couple of days but I can't get this problem to change. Tried debugging as well but can't identify the problem.

I'm a high school student trying to simulate a kilonova, so please keep that in mind in your responses.

Continue reading...