-
Notifications
You must be signed in to change notification settings - Fork 13
White Dwarf Binary
The merging of compact objects, such as neutron stars (NSs) and white dwarfs (WDs), is an interesting phenomenon for study in astrophysics. Here, we present a binary and example of a binary white dwarf (BWD) simulation using FleCSPH.
To set up the system, we begin by initializing two individual stars. For this example, we consider 0.45 and 0.75 solar mass WDs. Provide radial profiles for the WDs with the following normalized columns: r/Rtot, rho*Rtot^3/Mtot, m/Mtot, and (drho/dr)*Rtot^4/Mtot. Radial profile can be found in here (You can download this file or create WD_M0.49_prof.dat by yourself).
Create input_star1.par using below parameters:
#
# White Dwarf
# Initialization
# Use ./sedov_3d_generator <this_file>.par
#
initial_data_prefix = "wd_0.49_initial"
# geometry:
domain_type = 1 # 0:box, 1:sphere
sphere_radius = 9.824616359799e+08
density_profile = "from file"
input_density_file = "WD_M0.49_prof.dat"
# icosahedra lattice with small perturbations
lattice_nx = 59 # particle lattice dimension
lattice_type = 4 # 0:rectangular, 1:hcp, 2:fcc, 3:icosahedral
lattice_perturbation_amplitude = 0.10 # in units of sm. length
# equation of state type and parameters
eos_type = "white dwarf"
# density and pressure for relaxation stage
rho_initial = 1.822425000000e+06
pressure_initial = 6.665335965511e+22
# since we only need spherical distribution of particles,
# set Sedov energy to zero
sedov_blast_energy = 0.0
# use a good kernel
sph_kernel = "Wendland C6"
sph_eta = 1.6
Using this, we can generate initial data via below command
mpirun -np 1 ./sedov_3d_generator input_star1.par
Repeat for the second star configuration.
Run the WVT relaxation. Select either Diehl or Arth method in the input file. The number of iterations is currently set to 4000.
Create relaxation_wd_wvt.par using below parameters:
#
# White Dwarf Relaxation
# via WVT
# Use ./wvt_3d <this file>
#
# initial data
initial_data_prefix = "wd_0.49_initial"
initial_iteration = 0
domain_type = 1 # 0:box, 1:sphere
sphere_radius = 9.824616359799e+08
rho_initial = 1.822425000000e+06
density_profile = "from file"
input_density_file = "WD_M0.49_prof.dat"
# evolution
final_iteration = 4000
adaptive_timestep = false
# output
out_screen_every = 1
out_scalar_every = 10
out_h5data_every = 20
output_h5data_prefix = "wvt_wd_relaxed_m0.49"
# eos
eos_type = "no eos"
# sph
sph_kernel = "Wendland C6"
sph_eta = 1.6
sph_variable_h = true
# wvt
wvt_method = "diehl" # "arth" or "diehl"
wvt_boundary = "reflective" # "reflective" or "frozen"
wvt_mu = 1.e-3 # small fraction of smoothing length
wvt_ngb = 60 # number of desired neighbors
wvt_convergence_check = true
wvt_convergence_point = 0.1
wvt_h_ngb = true
wvt_cool_down = 0
wvt_radius = 9.824616359799e+08. # if problems occur at the edge, reduce to less than total star radius
Using this, we can generate initial data with the following command
mpirun -np <#> ./wvt_3d relaxation_wd_wvt.par
You can check the evolution of the particles with any data analysis and visualization application, e.g., ParaView. However, beware that the density that you will see is not accurate. It is only provided for rough cross-check. For an accurate density.
(OPTIONAL) extract last iteration from wvt_wd_relaxed_m0.45.h5part, with the following requirements: extract_h5part_iteration.py (in flecsph/tools) and wvt_wd_relaxed_m0.45.h5part by running
python extract_h5part_iteration.py -f wvt_star_relaxed.h5part -l
Finally, check that the star is stable under it's own gravity by running
mpirun -np <#> ./newtonian_3d evolution_wd.par
with the following parameter file, evolution.par:
#
# White Dwarf
# Evolution with FMM
# run by ./newtonian_3d <this_file>
#
# initial data
initial_data_prefix = "wd_relaxed_m0.49_<########>"
initial_iteration = 0
# equation of state
eos_type = "wd"
# sph kernel
sph_kernel = "Wendland C6"
sph_eta = 1.2
sph_variable_h = yes
# evolution parameters
final_iteration = 200
initial_dt = 2.e-14
timestep_cfl_factor = 0.5
adaptive_timestep = yes
thermokinetic_formulation = FALSE
# output
out_screen_every = 10
out_scalar_every = 10
out_h5data_every = 20
output_h5data_prefix = "wd_relaxed_m0.49_ev"
# gravity related parameters:
enable_fmm = yes
fmm_macangle = 0.8
gravitational_constant = 6.67408e-8
If stable, you can optionally extract the last
In order to accurately reflect orbital effects on the star's configuration, run the single star under gravity with the following parameter file, eforce_evolution.par:
#
# White Dwarf
# Evolution with FMM
# run by ./newtonian_3d <this_file>
#
# initial data
initial_data_prefix = "wd_relaxed_m0.49_ev"
initial_iteration = 0
# equation of state
eos_type = "wd"
# sph kernel
sph_kernel = "Wendland C6"
sph_eta = 1.2
sph_variable_h = yes
# evolution parameters
final_iteration = 10000
initial_dt = 2.e-14
timestep_cfl_factor = 0.75
adaptive_timestep = yes
thermokinetic_formulation = FALSE
# output
out_screen_every = 20
out_scalar_every = 20
out_h5data_every = 100
output_h5data_prefix = "wd_relaxed_m0.49_roche"
# gravity related parameters:
enable_fmm = yes
fmm_macangle = 0.8
gravitational_constant = 6.67408e-8
# eforce parameters:
external_force_type = "orbit"
mass_neutron_star = 1.4913525e33 # companion star
mass_white_dwarf = 9.916419542156e32 # primary star
orbital_separation = 3.064070749e9 # choose orbital separation
Run this with the following command:
mpirun -np <#> ./newtonian_3d eforce_evolution.par
Once you confirm that the star has relaxed into its Roche potential, you can extract the final iteration as before.
Run the following python script to place these stars into a binary. Requirements: make_binary_system.py (in flecsph/tools), star_1.h5part, star_2.h5part.
python make_binary_system.py -f <star1.h5part> <star2.h5part> -a <orbital separation in cm>
make_binary_system.py can be found in tools directory. After successful generation, you will have
binary.h5aprt