Setting up stars and tidal disruption events

Setting up and relaxing a star

First, follow the usual procedure for initiating a new simulation with phantom. We’ll use the “star” setup, but you can also use the “polytrope” or “neutronstar” configurations (the first two use self-gravity for the star, the last one uses an external potential). For tidal disruption events in general relativity use“grtde”. That is:

make a new directory and write a local Makefile

$ mkdir star
$ cd star
$ ~/phantom/scripts/ star > Makefile

compile phantom and phantomsetup

$ make
$ make setup
$ ls
Makefile    phantom*    phantomsetup*

run phantomsetup

 $./phantomsetup star

 star.setup not found: using interactive setup

 1) Uniform density profile
 2) Polytrope
 3) Density vs r from ascii file
 4) KEPLER star from file
 5) MESA star from file
 6) Piecewise polytrope
 7) Evrard collapse
Enter which density profile to use ([1:7], default=1): 2
Setting up Polytrope
Enter mass unit (e.g. solarm,jupiterm,earthm) (blank="blank",default="solarm"):
Enter distance unit (e.g. au,pc,kpc,0.1pc) (blank="blank",default="solarr"):
Enter the approximate number of particles in the sphere ([0:], default=100000):
Enter the desired EoS for setup (default=2):
Enter gamma (adiabatic index) ([1.000:7.000], default=1.667):
Enter the mass of the star (code units) ([0.000:], default=1.000):
Enter the radius of the star (code units) ([0.000:], default=1.000):
Relax star automatically during setup? (default=no): y
 Writing star.setup
STOP please check and edit .setup file and rerun phantomsetup

once you have answered the questions once, they are written to a file called star.setup, which you can use for subsequent runs of phantomsetup without having to answer prompts (see below)

relax the star

Open the star.setup file and make sure the parameter “relax_star” to True:

relax_star =       T    ! artificial damping of velocities (if on, v=0 initially)

Then run phantomsetup:

./phantomsetup star.setup

Which will generate a sequence of snapshots of the relaxation process:

RELAX-A-STAR-O-MATIC: Etherm:  0.499     Epot: -0.854     R*:   1.00
    WILL stop WHEN: dens error <   1.00% AND Ekin/Epot <   1.000E-07 OR Iter=0

-------->   TIME =    0.000    : full dump written to file relax_00000   <--------

 Relaxing star: Iter   1/1000, dens error: 20.39%, R*:  0.978     Ekin/Epot:  3.157E-03
 Relaxing star: Iter   2/1000, dens error: 15.91%, R*:  0.978     Ekin/Epot:  4.858E-03
 Relaxing star: Iter   3/1000, dens error: 13.05%, R*:  0.978     Ekin/Epot:  3.303E-03
 Relaxing star: Iter   4/1000, dens error: 11.11%, R*:  0.977     Ekin/Epot:  2.276E-03
 Relaxing star: Iter   5/1000, dens error:  9.69%, R*:  0.975     Ekin/Epot:  1.655E-03

-------->   TIME =   0.3756E-01: full dump written to file relax_00001   <--------
 Relaxing star: Iter   6/ 500, dens error: 22.17%, R*:   1.03     Ekin/Epot:  1.280E-03

This process will stop when either the density error is less than 1% and the ratio of kinetic to potential energy is less than the specified tolerance, OR the number of iterations reaches the specified maximum.

Once complete, you should obtain a relaxed initial conditions snapshot with the correct density and pressure profiles:

-------->   TIME =    0.000    : full dump written to file star_00000.tmp   <--------

 input file written successfully.
 To start the calculation, use:


evolve the star

If you really want to, you can evolve the star in isolation with the regular code to double-check that the relaxation process worked:


check the output

splash star_0*

(if you have version 2 of splash, the relevant command is “ssplash”)

Putting the star on an orbit for a tidal disruption event

If you used the “tde” or “grtde” setup then simply compile moddump:

$ make moddump

otherwise you need to specify the tidal moddump file:

$ make moddump MODFILE=moddump_tidal.f90

Then run moddump on your relaxed star:

$ ./phantommoddump star_00000 tde 0.0
 writing moddump params file tde.tdeparams
  Edit tde.tdeparams and rerun phantommoddump

When you first run this, a “tde.tdeparams” file will be created. Edit this to set the star on your desired orbit, and then rerun phantommoddump:

# parameters file for a TDE phantommodump
                beta =       1.000    ! penetration factor
                  mh =   1.000E+06    ! mass of black hole (code units)
                  ms =       1.000    ! mass of star       (code units)
                  rs =       1.000    ! radius of star     (code units)
               theta =       0.000    ! stellar rotation with respect to x-axis (in degrees)
                 phi =       0.000    ! stellar rotation with respect to y-axis (in degrees)
                  r0 =        490.    ! starting distance

After this you can simply run phantom:

$ ./phantom

Adding a magnetic field to the star

compile phantommoddump

The module used to compile this utility is specified using MODFILE= in phantom/build/Makefile. The default for the “polytrope” setup is currently moddump_spheres.f90:


Change this to moddump_default.f90. You can do this temporarily on the command line by compiling phantommoddump as follows:

make moddump MODFILE=moddump_default.f90 MHD=yes

run phantommoddump

$ ./phantommoddump
PhantomSPH: (c) 2007-2023 The Authors

 Usage: moddump dumpfilein dumpfileout [time] [outformat]

in our case we want:

./phantommoddump star_00010 magstar_00000

which will give some errors:

ERROR! MHD arrays not found in Phantom dump file: got            0

but then prompt you to add magnetic fields:

add/reset magnetic fields? (default=no): yes

you can follow the prompts to add uniform magnetic fields using this routine.

now implement something decent in src/setup/set_Bfield.f90

you can either use the pre-cooked magnetic field setups in this routine, or you can just make a new moddump module that sets up the magnetic field in a custom way.