Hey all,

Here is a run down on the things I have been working on. Pretty much I have been unit testing the initialization code. I have been matching everything generated to theoretical quantities. I will put up some plots next time around of my fits. But so far I have been able to fit the radial and angular distributions. This means that the particles are being placed in the correct place.

What I am currently working on is fitting the potential and the velocities. These are especially important because if they are off, then the evolution of the particles will be wrong.

We want our dwarf galaxy to be in virial equilibrium. Mathematically speaking, we want the total kinetic energy and potential energy of all the particles to follow: 2K/P = 1 where K is the total kinetic and P is the total potential energy. Thus, if the particles have bad velocities they would not be in virial equilibrium, i.e. the ratio will not equal 1.

I am close in fitting the potential, but it is not exactly matching the theoretical value. So there is work to be done on that front.

Also, I am working on a way to test the velocities as well, so that is also a work in progress.


I think it is funny to note that my test codes combined is now longer than my dwarf initialization code.

More to come soon.

Cheers,
Sidd

Hey All,

Sorry for the long delay in posting an update. So lets just jump right in. (Also, apologies if I start getting a bit too detailed.)

Disclaimer: the plots shown here are most likely not in a publishable quality/format (would have to better label things like axis and titles, etc). I am including them here only to give an idea of what I am working on.

I mentioned in the last blog post that the distribution matches the radial distribution. In 3d space, we need 3 coordinates. In our code, we assign radii to our particles by sampling the density distribution we want, in this case a plummer sphere. We also assign angles, the azimuth and polar angles from spherical coordinates. From these angles are are able to calculate the required x,y and z coordinates for that particle. We have checked these distributions and they match theory. These are picture below. The particle distributions are binned in the histogram and the theoretical distribution is shown with a line.



Above are the initial distributions for a two component dwarf separated into its light and dark components. The binned data is from the simulation and the curves are theoretical distributions.



Above are the initial spacial angular distributions for the combined two component sphere.


We then tried to do the same with the velocity distribution. As I mentioned in the last post, this is especially important. If the velocities are assigned too high, the particles go flying off into space. If they are assigned too low then they collapse into the center of the sphere. So we definitely need to make sure they match theory. One small problem, we do not know what the theory curve looks like.

In order to remedy this, we look at special cases. Our code right now allows us to produce a two component plummer sphere. However, in the special case, one of those spheres has 0 mass and no particles are assigned to it. Therefore, using the 2 component structure, we produce a 1 component sphere. How does this help?

In order to assign velocities, we sample something known as the distribution function. This function is very complicated. This function allows us to find the probability of finding a particle at a particular location and velocity by integrating it over space and velocity. Therefore, if one were to integrate it over all space and over all velocities you should get 1 (or 100% probability), since we know the particle exists and is somewhere with some velocity. The general version of this equation involves an complicated integral shown below:


Above is the distribution function. The epsilon in the integral limit and integrand denominator is the energy of the particle. The rho is the density distribution we are using, and psi is the potential. These are shown below:


Where the a's are scale radii for the light and dark components, the M's are the masses of those components, G is the gravitational constant, and r is the radii.


We cannot solve this analytically because of our two component model. However, for a single component plummer sphere, the integral is readily solved and can be found in various texts. Therefore, we test our code by having it produce a single plummer sphere and comparing it to a single component model. In fact, doing this opens a new realm of testing.

Here is the single plummer distribution function:




The first set of tests, before anything else, was to make sure we are calculating the integral correctly. In order to solve the integral, we expand the integrand in terms of r. For the limits, we find the radius at which the potential is equal to that energy (using a root finding algorithm. Originally written by Jake Weiss). Therefore, the integral is actual a function of r and v of the particle and is integrated over r' as shown below:



The first test of the integral is to hold v constant and alter r and calculate the integral. Then compare it to the single component equation when doing the same thing. This is shown below:


The red curve being the one calculated from our two component model
The next test was to hold r constant and alter v. This is shown below. Again, red is our two component model.


Originally, the two curves, calculated and theoretical did not quite fit. This was because of the integration function we were using as well as the form of the integrand. The integrand has a singularity in the denominator. As the curve moves away from the singularity, it quickly tends to zero. We wanted this section of the integral to be calculated more accurately. Simply increasing the resolution of the integral solved the issue but made the initialization procedure extremely lengthy. Therefore, we decided to put in an adaptive resolution so that the integral would be calculated very accurately where the integrand "mattered" and not as accurately where it was practically zero. This led to a decent initialization time (about a minute to 2 on my laptop) and the above plots.

Below is a plot of the integrand. The part that is actually integrated is to the right of the spike:


Finally, back to the velocity distributions. First off, just to quickly get them out of the way, we were able to match the velocity angular distributions to theory:





Testing the velocity distribution was a bit more complicated. Since we know the single plummer distribution function, we can create a sample set of velocities by sampling that equation. We can then compare that distribution to the one created from our two component model (remember, it is still being used to generate a single component sphere, so they should match within random seed differences). We were able to see that they match:


The histogram is our calculated distribution and the green line is the test distribution.


As for the potential I mentioned in my last post, we were able to find multiple things. What we were looking to do was match a theoretical value of the total potential energy of a two component plummer sphere against what we actually calculated. I do not have plots for these, but we were able to find that the theoretical values did closely match what we get from our generated particles. To be more specific, we calculated the total potential energy using the following method:

The first term after the equal sign is the normal method of calculating the total potential energy of a group of particles. By changing this into an integral, we can use our theoretical distributions.


Finally, we were able to (and was actually the first thing done after the last post) to get the virial ratio to be equal to 1. We calculated this ratio two ways. In both ways we calculated the total kinetic energy by summing over the kinetic energy from each particle. In the first way, we calculated the total potential energy using the following method:

This takes the mass of each particle and multiplies by the theoretical potential we are using to calculate the potential energy.

The second method uses:

This uses the classical method of calculating the potential energy, by calculating the total potential energy each particle gets from every other particle (and divided by 2 to compensate for double counting).

Thankfully, both methods yielded the same answer, and the virial ratio, 2K/U, was found to equal 1 which is what we want for an object in virial equilibrium.


So, now we know everything is being initialized correctly. However, we want to make sure they stay that way. Therefore, we replot each of the above plots over time. (My code that does this stability test produces 112 plots, all of which are useful!). Ideally, all of the distributions would stay mostly the same over time. We find that the angular plots do. However, the radial plots show a tiny bit of movement:


Here are the radial plots initially, after 1gy and 4gy of run time:





and the velocity distribution:





As you can see, there was some movement at 1gy. The movement was a bit worse previously which was due to two things. The first was that some plummer sphere mass and radii parameters are intrinsically a bit unstable. The second was that the parameters being sent into the calculation of the time step were being sent in incorrectly. After these two things were remedied we are able to get the above stability plots. However, this is still too unstable for our liking.

What we are currently looking at is if the timestep equation is correct. We are comparing this against established code such as NEMO, which has the ability to create a single plummer sphere and evolve it. However, I do not believe it will be much longer before we are able to release a new version of NBODY.

Sorry for the long post. I do not have any potatoes for you. Also, I did not check the grammar, so forgive me for any errosr.

--
Cheers,
Siddhartha Shelton