Lack of convergence of neutron star initial data

[Pedro Espino]

I'd like to report an issue with the convergence of neutron star initial data produced by some thorns packaged with the toolkit. Firstly, I assess the converge of TOV solutions using polytropic equations of state by considering whether the Hamiltonian constraints converge to zero at the expected rate (using the fourth order finite differencing scheme in the ML_BSSN thorn). I consider Cactus grids with 64. 80, and 100 grid points per NS radius. I find that TOVSolver does not produce convergent ID at the default resolution, and only after the resolution has been turned up significantly (e.g, setting TOVSolver::TOV_dr[0] = 1e-6, TOVSolver::TOV_Num_Radial = 1e7 in the cases I considered) did I find the expected fourth order convergence of the Hamiltonian constraints within the star. The story is worse with RNSID, as there seems to be a problem interpolating the solutions onto the Cactus grids. With RNSID, I set the solution accuracy to be as low as possible (e.g, setting RNSID::accuracy ~ 1e-16), for the model in question but could not obtain the desired convergence rate. Due to its similarity to RNSID, I have also tested a private thorn which imports NS ID built by the Cook et al. code, and could not obtain the desired convergence there either. The Cook et al. code uses a compactified radial coordinate and polar coordinates similar to those of RNS, and this may result in poor interpolation onto the Cartesian grids of Cactus.

Besides NS ID built assuming a polytropic EOS, I have also tested stars built with RNSID using realistic, tabulated EOS with the goal of using cold beta equilibrium versions of the finite temperature tables found on stellarcollapse.org. To have a controlled environment, I generated piecewise polytropic tables, for which I know everything analytically. When building these stars with RNSID, the EOS interpolation appears incorrect and if I check if the EOS is satisfied, I find errors up to few percent. Increasing the number of points in the table to large values like 70,000 improves the situation but does not solve the problem. Errors of a few percent in the internal energy changes ID which is supposed to describe cold NSs to warm or hot profiles instead.

The convergence rate of different hydro variables has also been tested for the expected second order convergence using different MHD thorns (GRHydro and IllinoisGRMHD), and the same conclusions can be made: TOVSolver produces convergent ID at very high resolution, and RNSID cannot produce convergent ID. Below I attach a few plots which show the convergence (or lack thereof) of the constraints for a TOV star with central rest mass density rho_central=0.129285 , built with TOVSolver, RNSID, and the Cook et al. code assuming a Gamma=2, K=1 polytropic EOS.

I’ve spent some time trying to fix these problems on my own and have tried many things including changing the internal grid structure in RNSID to allow for more radial grid points, increasing the solution resolution as much as possible, using a high number of points in tabulated EOSs, using very high resolution Cactus grids (about twice the resolution of the results for the grids shown in the plots), and more, but could not achieve the desired convergence rate with RNSID or the similar Cook et al. code. I am currently working to expand the unreleased TOVSolverHot thorn to work with finite temperature EOSs (out of the box, it seems to be incompatible with the nuc_eos EOS type within EOSOmni) which will hopefully converge at the same rate as TOVSolver.