Large shift of the center of mass of the star

[cradesto]

During the relaxation process for long times ($\sim 10^{4}$ [unit]), the center of mass of the neutron star is shifted by a distance comparable to the radius of the star. To reproduce this issue one need to execute:

$PHANTOM_DIR/scripts/writemake.sh nsmerger > Makefile
make
make setup
./phantomsetup star

During the 'setup' action one should choose the following options:

cat star.setup

# Phantom v2022.0.1 (c) 2007-2020 The Authors
# input file for Phantom star setup
            iprofile =           2    ! 1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Piec, 7=Evra

# units
           mass_unit =      solarm    ! mass unit (e.g. solarm)
           dist_unit =          km    ! distance unit (e.g. au)

# resolution
                  np =       50000    ! approx number of particles (in box of size 2R)
          use_exactN =           T    ! find closest particle number to np

# equation of state
        use_var_comp =           F    ! Use variable composition (X, Z, mu)
                ieos =           2    ! 1=isothermal,2=adiabatic,10=MESA,12=idealplusrad
               gamma =       2.000    ! Adiabatic index
                  mu =       1.000    ! mean molecular weight

# star properties
               Rstar =         10.    ! radius of star
               Mstar =       1.000    ! mass of star

# relaxation options
          relax_star =           F    ! relax star automatically during setup
   write_rho_to_file =           F    ! write density profile to file

These parameters are selected to reproduce a typical neutron star.

Here we choose not automated relaxation of the star. So one needs to run ./phantomsetup star again. This will create star.in In the last file, we change the tmax value from 100, which was set by default, to 3e+04.

After that run command

./phantom star.in

Results are pictured on the following figure:

The rms velocity of sph-particals is smaller than $5 \cdot 10^{-3}$

The coordinates of center mass of the star ${x{com},y{com},z_{com}}$ is greater than 30 km at the end of the calculation.

The star seems to shift as a whole, although the velocities of all particles are small (see following figure).

Where can there be a mistake?

Thank you for your work.

[jameswurster]

Hello Sorry for the delay in my response.
I have been able to reproduce your issue -- thanks for clearly providing your setup files & information.

I first ran the default case, and I found similar drift, however, that drift was on order 1e-3 compared to the spatial scale on order 1; the drift for the default case is unimportant.

The smaller radius in your case is amplifying the issue. Moreover, I have seen similar drifts in my own star formation simulations. I have tested several patches using your setup without any success. I suspect that any slight asymmetries in the sphere -- either caused by the initial lattice, particle groupings (for gravity or openmp), or computational roundoff, ultimately leads to the shift. At the moment, I don't have any viable patch, however, it seems reasonable that if you 'simply' remove the shift during post-processing, then you should have reasonable data to work with.

I know that this isn't that this isn't the answer that you are looking for, but you can be assured that it is not a mistake on your part.

Cheers

[danieljprice]

I'm not understanding why the automatic relaxation is not being used here.

The best way to fix this is to use a smaller value for the tree_tolerance parameter. The tree code does not exactly conserve momentum and this is the main source of the drift