N-Body solver

N-Body solver#

Introduction#

One feature that FARGO3D has inherited from its ancestor FARGO is the possibility to include a number of point-like masses which can interact between themselves and which can also gravitationally interact with the gas. The set of such masses is defined in a file. This file must be defined in the parameter file with the parameter PlanetConfig. For instance, in the default parameter file of the fargo setup (that is, setups/fargo/fargo.par), we find the line:

PlanetConfig    planets/jupiter.cfg

Throughout this chapter, we call this file the config file, or the planetary configuration file. The content of this file is strictly similar to that defined in the previous FARGO code:

###########################################################
#   Planetary system initial configuration
###########################################################

# Planet Name         Distance        Mass     Accretion      Feels Disk      Feels Others
Jupiter                1.0            0.001    0.0            NO              NO

Lines beginning with a ‘#’ are ignored. A line containing a valid planet definition must begin with a letter (‘a-z’ or ‘A-Z’). So far this name is not reflected in the code output, but it gives a sense of what kind of planet or planetary system is intended. The second column gives the semi-major axis in units of R0, defined in src/fondam.h. The third column gives the planetary mass in units of MSTAR, defined also in src/fondam.h. The fourth column, relative to accretion, is not used at the present time. So far its only purpose is to allow backward compatibility with FARGO’s planetary config files. Finally, the content of the last two columns, which is self-explanatory, indicates whether the planet feels the gravity of the other planets and whether it feels the gravity of the disk.

Note

Our N-Body solver is a simple fifth order Runge-Kutta integrator. This may seem pretty simplistic compared to way more sophisticated N-body integrators, which are tailored to follow the dynamics of a collection of point-like masses in vacuum over very long timescales. It should be remembered that our primary purpose is the study of protoplanets embedded in disks. The lifetime of the disk is 3 to 4 orders of magnitude shorter than the age of the Solar System, and the variation of total energy of planet under a disk tide is usually many orders of magnitude larger than that inherent to our scheme. The latter is therefore perfectly fine for our purpose.

The planets are initialized at \(t=0\) such that they lie on the \(x\) axis, are prograde with the disk, have all same eccentricity (defined by the parameter ECCENTRICITY, which defaults to zero). Their initial location corresponds to apoastron.

Warning

By default, a central star of mass MSTAR (defined in src/fondam.h) is added at the mesh center. This default behavior can be superseded, for instance, to have the mesh centered on the center of gravity of the system, or to have a binary star at the center.

The value of MSTAR is used to initialize the orbits of the planets, in all cases.

Removing the default central star#

In some cases, you want to remove the default central star which is assumed to exist at the mesh center. In order to do so, so simply need to add the following line to your .opt file and rebuild:

FARGO_OPT += -DNODEFAULTSTAR

Planets which are stars#

Starting from version 1.2, FARGO3D detects in the configuration file of the planetary system which objects are massive enough to be considered as stellar, and deals with them in a specific way. Our arbitrary threshold, defined in src/fondam.h, is 0.05*MSTAR, which would correspond to a typical brown dwarf if MSTAR is one solar mass. This threshold can be redefined depending on the user’s needs. The number of stellar objects detected in the configuration file can be 0, 1 or 2. Beyond this number, an error message is issued.

The array below summarizes a variety of possible setups. The first line corresponds to the default behavior. The last two columns indicate whether the user should define NODEFAULTSTAR (in the opt file) and what value should be given to MSTAR (in src/fondam.h). Note that whenever NODEFAULTSTAR is defined, all stellar objects must be explicitly given in the configuration file.

Intended setup

NODEFAULTSTAR

MSTAR

One star at mesh center

mass of central star

One star. The mesh is centered on the center of mass of the star and planets.

Defined

mass of central star

Circumbinary disk

Defined

mass of binary

Circumstellar disk in binary system

mass of star at disk center

Note

What is the reason for defining manually MSTAR in the second and third rows of this table ? The reason is that this variable is used throughout the code to perform a variety of tasks, such as:

  • initialization of the planetary orbits

  • evaluation of the Roche radius (for potential smoothing issues)

  • evaluation of the planet’s orbital elements

so it is important that its value matches that of the mass of the object(s) at the disk center. Also, note that the orbit[i].dat output files contain correct values only for the first kind of setup.

When one uses the third setup in the array above, the code enforces that the mesh be centered on the center of mass of the binary (and not of the system), and the components of the binary start do not feel the gas or the planets. The indirect term, if included, arises from the acceleration of both components of the binary, weighted by their mass. See the setup binary for an example.

Initialization of the position and velocities of the stellar objects#

The initial position and velocity of the stellar objects differ from that of the planets, in the case in which there is no default central star.

If one stellar object is given, its position and velocity is set so that the center of mass of the planets and star is at the mesh center, and so that this point does not move with respect to the mesh.

If two stellar objects are found, the parameters of the second column of the config file, usually devoted to the distance to the star, takes a different meaning. For the first star, the value found is used as the binary period, and for the second star, the value found is used as the binary eccentricity. The binary system is initialized on the \(x\) axis, at apoastron, and the first star has \(x>0\). The mesh is centered on the center of mass of the binary. See the file planets/Kepler38.cfg for an example.

Indirect term#

Whenever the center of the frame is accelerated, an additional potential term, called the indirect term, is added to the potential, which corresponds to the fictitious force arising from the frame acceleration. This indirect term can be purposely discarded, by setting the parameter IndirectTerm to NO (the default is YES).

Initial eccentricity and inclination#

All objects are initialized with the same eccentricity, defined by the ECCENTRICITY parameter (its default value is zero), and they are initially set at apoastron. If a different setup is needed, some edition of the file src/psys.c is required. Alternatively, one may edit manually the planet files (planet[i].dat, see also Planet files).

In a similar manner, all objects are initialized with the same inclination, defined by the INCLINATION parameter (its default value is zero, and it is given in radians). The \(x\)- axis (\(x>0\)) on which the planets are set at the start of simulation corresponds to the ascending node: the planets are initialized with \(z=0\) and \(v_z>0\).

Convenience parameters: PLANETMASS, SEMIMAJORAXIS and ORBITALRADIUS#

In order not to have to edit the planetary configuration file (which can be useful to spawn easily a parameter space exploration), we have defined three convenience parameters: PLANETMASS, SEMIMAJORAXIS and ORBITALRADIUS. Their default value is 0. When the first two (PLANETMASS and SEMIMAJORAXIS) differ from zero, they are used to supersede the values given in the planetary configuration file for the first planet only. When the last one (ORBITALRADIUS) is set to a value different from zero, this value is used to rescale the semi-major axis of all planets prior to the run start. This parameter is, therefore, a dimensionless factor that can be used to expand or shrink the whole planetary system. This is mainly useful in parameter space explorations.

Correcting the shift of resonances#

The gas in FARGO3D is not self-gravitating: it orbits in the potential of the star exclusively. On the other hand, a planet which feels the gas gravity orbits in a slightly different potential, which is the sum of the stellar and disk potential. The Lindblad and corotation resonances of an embedded planet are shifted (with respect to a more consistent situation in which the gas would be self-gravitating) as a consequence of this mismatch. This has sizable effects on the measured migration velocity of a planet moving freely in the disk, as shown by Baruteau and Masset (2008). A workaround to this problem consists in removing the azimuthally averaged density to the density of each zone prior to the force evaluation: if the disk were unperturbed, the planet would exclusively feel the star’s gravity. In order to activate this workaround in the code, one has to add the following line to the`.opt` file and rebuild:

FARGO_OPT += -DBM08