Can of worms

In Walter's comments there was this observation: Figure 4, which plots a fit to the fiducial model, does not agree with the table, which also has a fit to fiducial model.  At first I thought it must have been something simple, like possibly I was not using the same domain of the data in the fit, or possibly fitting one parameter instead of two, something like that.  The answer turned out to be quite a bit more interesting.

It turns out, when we switched to the symmetric version of the code, and I re-ran all of the simulations from scratch, I forgot to halve the number of particles in the resolution study. So the simulations we have, instead of being [1k,2k,5k,10k,20k] actually correspond to (in our convention) [2k,4k,10k,20k,40k].  So it seemed like an easy fix.  Not so, the slope in the plot still did not agree with the 5k entry in the table.

I looked into why this was.  It turns out that I had made a resolution scan directory, and I had re-run the fiducial model (the 10k case) in that directory, separate from the fiducial model in the initial-condition scan.  Further inspection revealed that the two sets of initial conditions were not actually identical... one IC0 case had particle numbers (5003,5006,5009,5012,5015,5018) whereas the the other case had particle numbers (5004,5007,5010,5013,5016,5019).  Given the smallish dispersion between six sims of similar particle numbers, it seemed weird that the derived slopes should be so different for the realizations of the fiducial initial conditions.   Looking into this led me to make the following plot:

Odd that (given how many dozens of plots have been made) that we haven't seen a plot like this, but I think it is significant.  The x axis is time, and the y axis is the fit, done using the six simulations.  Red and blue are the two different sets of six.  Although the bootstrap tells us that for a given time step, the fit to $\rho(x)$ is very stable (i.e. always the same regardless of which subset of the $\rho(x)$ dataset we take, and stable to moderate changes in the range of data), the fit for different time dumps jumps around considerably.  Jerry thought he saw some kind of beating or periodicity in this plot, so I restarted at time 640, and dumped 16 times as often.  Here are the resulting plots:

Its unclear that these plots teach us anything more than the one above, the points are certainly not interpolating the data above, although I have checked that the at points that coincide, the fits agree precicely, and they do.

A long time ago, I suggensted using time averages, and made plots averaging several time steps.  Scott nixed it, here is the thread to remind you:

     The fit versus time plot suggested to me that the fits might stabilize
significantly if I were to average together several timesteps, as I do for the
clumpy warm cases.  Attached are some screen shots with results from this
procedure (average of 6 runs and 20 time steps).  One is the new table with fits
for the different resolutions, another is the plot where we vary the inner cut,
the third is the plot showing  the density scaled by the best fit.   These seem
cleaner to me, so let me know if you'd like me to use this method for all of the
fits...

I'm actually *not* enthusiastic about using the time-averaged plots that you showed me. The reason is that in these plots the error bars (the yellow bands) appear to be much narrower than the random bumps in the density plots (e.g., the bump at 0.15 in IC 2). I confess that I find this surprising---I think the yellow bands come from averaging several different simulations; why should they all show the same bumps and wiggles? Presumably this was true in the non-averaged plots as well, but didn't show up because the error bars were larger.

At this point I am not keene on redoing every plot again, so I am thinking of adding this plot, to give people a flavor of the kind of variability these systems exhibit....  Let me know what you think.  

Action plot etc.

I made the plot Walter asked for, but am hesitant to put it into the draft for two reasons 1) It's unclear what the "action" really means for times below around 40 or so, since when we compute it we assume the potential is stationary whereas in fact it is evolving, and 2) it doesn't tell us anything distinct from the energy vs. time, as far as I can tell.  Here is the plot:

I've made this plot before, but settled on energy instead of action for the draft, mainly because of reason 1).  It's your call though, we can throw this into the paper if you all agree it should go in.   Incidentally, I also took this opportunity to make a similar plot for one of the steep cases, with $\gamma=0.8$, for comparison. Here it is:

This demonstrates pretty clearly, that in comparison to the fiducial case, these actions are frozen almost exactly where they started.

I thought you might appreciate seeing the direct comparison for the energy plots, I know I've done this before but I couldn't remember if I'd shared with everyone.  These plots were the basis for my statements in the draft that the energy and actions for the steep cases evolve only minimally from the initial conditions.

Although the energies do redistribute a little more than the actions do, it's pretty clear that the density profile is not going to be able to change significantly, because the orbits do not have the opportunity to gain or shed energy from their original energies.  Another way to put this is that for the steep cases, the characteristic distance over which energy can be transported is much smaller than in the fiducial case.

Update

The plan, over a month ago, was to re-make all of the figures using one of the codes that stays symmetric, so that we would not have to do any re-centering shenanigans.   I had also wished to make some plots with the evolution of the gravitational potential, and the evolution of the actions, using the new action capability.  I have made many many plots over the last few weeks but the difficulty is the following.  I have 4 versions of the code, three of which are broken, and one of which does not stay centered.  Let us label them.

A) Walter's code, original 1Dgrav code that I was using since last fall - does not preserve symmetry
B) Walter's code, updated, Grav1D that he provided a couple of months ago
C) Walter's code, explicit symmetry enforced, included as an option in Grav1D distribution
D) Jerry's code, explicit symmetry enforced, includes both action and angle outputs

In addition there is Walter's action post-processing step. 

All the plots in the draft were made with code A.  A majority of these results, made with code A have been double checked with code D, however, code D cannot run to late times.  There is an intermittent bug that causes the system to explode.  This bug is difficult to track, because it happens after a long time, and if the particle state is dumped just before the explosion and the code re-started, the system does not explode in the same way at the same time.  At present, it does not seem to be an option to use code D to regenerate all the plots.

Code C is problematic, either I have misunderstood the way the initial conditions are meant to be set up, or it simply is not doing what it is intended to do.  Code C takes initial conditions on the positive x axis only, and assumes explicit symmetry for particles on the negative half.  I will email files for all initial conditions, but here is a plot of  the initial condition I put in (which is 1/2 the sign wave that I had put in before):

After about 20 time steps the system looks as follows:

Which is clearly not monolithically collapsing about x=0 and v=0.   I don't believe this is the intended behavior.

Code B, which I had decided to use instead of code A to generate all the plots (because in the runs I checked, it stayed centered to around 5 decimal places), is also not working but it took me much longer to realize it.  Unfortunately at the moment, I am having trouble finding the output that made me conclude that it was broken (I've been combinatorically comparing 12 runs from 4 codes, with 700 time dumps each, it's totally miserable) but at a minimum, code B is not staying symmetric even approximately for some of the runs, for example, here is an output from code B for the density profile, without the centering iterator.

   At this point I am getting frustrated because this process is not leading to any new science output and I have used up a truely unwise amout of time trying ot sort it out.  I'd like to move ahead on the paper with code A, which is the original outputs for the paper.  I believe code A because I have confirmed that at least in the early stages, code A and code D agree with one another particle by particle, and those codes are independent.   It's unfortunate because I will not be able to do some of the tests I had hoped, because the iterator is simply not good enough at finding the center to do them.

Just for completeness, here's what's up with code D, sometimes early, sometimes late, the particles fly off to infinity:

A few questions:

1) Walter, if there is something easy and obvious that you can suggest to fix code C, then I'd be willing to try one more time with that one.   The same goes for Jerry with code D.

2)  For the runs that are not centered, do the action integrals that you compute make any sense... are they still valid?

Both codes behaving strangely

After trying to get somewhere with the perturbation theory approach to the stability analysis, I've taken a break and decided to see if I can use the numerics to address the same question.  I have set up initial conditions as follows:

1) Pick a power law scaling for the density that I wish to study
2) Solve for the distribution of energies that is expected for an equilibrium solution with this scaling
3) Generate some random energies distributed according to this power law
4) Use these energies to generate positions, a) requiring that the orbits are equally destributed in time over a quarter orbit and b) using potential energies consistent with the expected power law scaling
5) generate velocities using v=sqrt(2E-2phi(x))
6) reflect the (x,v) pairs over the x and v axes to generate symmetric initial conditions
7) run the code to see if the scaling remains the same, or if it wanders to the attractor

Here is a plot of the initial phase space configuration for a $\rho \propto x^{-0.8}$.

This initial configuration, which has been constructed to be in equilibrium, actually explodes, the particles quickly go off to \pm inifinity essentially.   Below is an intermediate phase space configuration, that to me seems clearly wrong:

And a little while later:

My next step, since we have an independent code, was to see if it is behaving in the same way.  Here is the initial condition I fed to Jerry's code:

And here is what it looks like a little while later:

So they are both suffering from a similar problem, I think, although for Jerry, it is the velocities that are exploding (which must EVENTUALLY cause the system to explode, but I didn't run it that far). Some of the other power laws are not behaving quite so badly, but in terms of examining whether an equilibrium is stable or unstable (and running off to the attractor solution), I feel we need to understand what is going on here a little better. 

I will email the initial condition file that causes the bad behavior, this time it is definitely reproducible. 

Some Russians have done this calculation

In chatting with Andrei Gruzinov, he mentioned a paper that he'd read from 1995 in which these questions we are investigating are explored.  The one he showed me was all in Russian, and it was hard for me to tell what was going on.  However, today I dug up an English translation.  You can find it here:

Gurevich and Zybin

The whole paper is very nice to read, these people have some serious chops (although Andrei claims the cosmology part of it is a bunch of nonsense).  The relevant sections I think are 4 and 5, and the scaling of the density they find is in eqn. 79.   They have started with different initial conditions, but using an adiabatic invariant, they have derived the scaling of the density with x, which is not $\rho(x) \propto x^{-1/2}$ but rather $\rho(x) \propto x^{-4/7}$.  I believe their analysis does not hold during the initial period of violent relaxation, but is valid thereafter.

I threw their scaling (which is slightly steeper) onto a plot, to compare to the numerical results:

By eye theirs looks a little steep, but it's not significant compared to our error bars.  If we were going to try and say something about this quantitatively, all of a sudden this error bar discussion we have been having becomes really important.  Jerry pointed out that in the scheme I have been using, the dispersion in adjacent bins is extremely correlated, because if a particle leaves a bin, it automatically shows up in the bin directly next to it.  This means that treating the error bars as uncorrelated in a chi-square fit will result is an anomalously good fit.   A possible error model (which I can check numerically):  it is unlikely that a particle will ever "jump" a bin, that is, the dispersion always causes it to join the bin next to it.  In the same way we measure the variance in each bin between all the simulations, we can measure the covariance between all adjacent bins.  We can use a chi square estimator where the covariance matrix looks like the diagonal plus one off diagonal term either side of the diagonal.  This will cause the chi square values to be about 3x what they are now (which is more consistent with the lousy K-S test values we were getting once upon a time).  It may be that one of these two models will fit significantly better than the other.

Obviously finding this paper was exciting and also disappointing... I am still trying to figure out if what we have adds anything to this story, for example that the scaling is in place before the adiabatic approximation becomes valid (actually, I am not sure why they claim it becomes valid right after the first caustic, anybody understand this?)...  The different behavior of very steep initial density profiles is also something new (but it would be nice to explain the origin of this).  Anyhow, any thoughts about what I should do now would be greatly appreciated.

 

Trouble with energies resolved

Some time ago, when I had a discussion with Jerry about how to find the density maximum.  He was confused about why the density maximum was not located at the median position, and also what was causing symmetry to be lost.  This led him to write an independent code to solve the problem.  His code does not explicitly enforce symmetry (he does particle pushes on both sides of zero).  However as far as we have been able to tell, the density maximum remains at x=0 and the particle positions are symmetric (though we have not run to as great a time as I have run Walter's code). 

I have compared most of our results between the two codes, and the good news is, most of the key points are the same in both.  One of my favorites, though (about the energy) showed a discrepancy.  When I plotted the energy versus time plot from the output of Jerry's code (for the same initial conditions), I got the following:

This looks qualitatively different than what I have been getting with Walter's code:

This worried me quite a bit, because you'll note that the tendency of the particles to stay sorted in energy is not true in both versions.  I have spent about a week attempting to track this down, and today Jerry and I used runs with 8 particles each to try and sort out what was going wrong.  We found that although the positions and velocities agreed between the two codes, the energies did not.  This led us to total up

$\sum_{\rm particles} \frac{1}{2}v^2 +p$

And compare it between time steps.  Although the log file shows that energy is conserved (and I think it definitely is), this sum is not the same between time steps.  We discovered that the p being printed should be divided by a factor of two, in other words, if the force is defined as $(M_R-M_L)$ then the potential energy of the ith particle needs to be

$p_i =\sum_{j} \frac{1}{2} \left| x_i - x_j \right|$

This means all the energy plots I have made so far are wrong, and in particular, the monotonicity result is probably going to go away.  I'll remake them tomorrow, and we can see. 

On the bright side, I think that was the only outstanding problem that I knew about regarding the data and simulations.   But it's making me sad, even though the fact that you can SEE the phase mixing in the E vs. T plot is really pretty cool.  Sigh. 

A puzzling thing

Say you start with the Vlasov-Poisson Equation

$\frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + a \frac{\partial f}{\partial v} = 0$

where

$a=-\frac{\partial \phi}{\partial x}$

and

$\frac{\partial^2 \phi}{\partial x^2}= 4 \pi G \rho$

Now say you take a moment of this eqn... multiply by $v$ and integrate (which ultimately gets you to one of Jeans' equations):

$\int dv \left[ v \frac{\partial f}{\partial t} + v^2 \frac{\partial f}{\partial x} + va \frac{\partial f}{\partial v}\right] =0$

$\frac{\partial}{\partial t} \left[\left< v \right> \rho \right] +\frac{\partial}{\partial x} \left[\left< v^2 \right> \rho \right]+ \int dv \frac{\partial}{\partial v} \left[ a f v \right] - \int dv \,\, af = 0$


$\frac{\partial}{\partial t} \left[\left< v \right> \rho \right] +\frac{\partial}{\partial x} \left[\left< v^2 \right> \rho \right]+ \left a f v \right|_{-\infty}^{\,\, \infty} - a\rho = 0$

$\left a f v \right|_{-\infty}^{\,\, \infty} \rightarrow 0$ as long as f vanishes at $\infty$ (kindof strongly).

It's not 100% clear that this is okay, but if we assume we have reached a stationary solution then we can set $\frac{\partial f}{\partial t}=0$ and if we define $\sigma^2=\left< v^2 \right> $ then we are left with the following:

$\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) - a \rho =0$

This is just a statement that for steady state, the gravitational force inward has to balance the pressure force outward.  But here, things get a little weird.

Suppose we assume that our attractor solution is a steady state solution of this type.  Then we know $\rho=Ax^a$.  We can try plugging this in (this weird thing came up because I was trying to solve for the scaling of $\sigma$ that we should observe in the simulations if the system is obeying the V-P limit).

$\frac{1}{\rho}\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) = a = -\frac{\partial \phi}{\partial x}$

$\frac{\partial}{\partial x} \left[\frac{1}{\rho}\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) \right] = -\frac{\partial^2 \phi}{\partial x^2}= -4 \pi G \rho$

$\frac{1}{\rho}\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) = -4 \pi G \frac{1}{a+1} A x^{a+1}$

$\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) = -4 \pi G \frac{1}{a+1} A^2 x^{2a+1}$

This is a pretty strange place to be at.

Observation 1) It turns out, there actually is something special and different about the power law scaling $\rho \propto x^{-1/2}$.  The LHS is the gradient of the pressure, and the RHS is the gravitational force.  If $a=-1/2$ then the gravitational force is constant as a function of position from the center of the system.

Observation 2) There is the distressing fact that the sign of the coefficient is negative.  I have looked and looked, but I haven't found an error yet.  When I showed this to Jerry he pointed out that probably the constant of integration in the next step cannot be ignored, it's the pressure at the center.  If you integrate this guy one more time you get:

$\left( \sigma^2 \rho \right) = -4 \pi G \frac{1}{2(a+1)^2} A^2 x^{2a+2} + C$

Where here it's not clear if C is even finite, but if it is, I have no idea how to get it... it's $\rho(0)\sigma^2(0)$, whatever that is. 


Observation 3)  For $a=-1/2$, $\phi\propto x^{3/2}$ and is therefore growing with $x$.  This is opposite from the 3D case, where usually $\phi$ decreases with radius, e.g. $\phi \propto r^{-1}$.  Consider the force due to pressure:

$-\frac{\partial}{\partial x} \left( \sigma^2 \rho \right) = \rho \frac{\partial \phi}{\partial x}$ 

The RHS has a different sign in 1D and in 3D.  It looks like the pressure in 1D points in the same direction as the force, which is probably not news to anybody else but it does mean I am completely confused as to how we managed to arrive at a steady state in the first place.        


Observation 4)  Turns out, sigma is not a power law in 1D as previously thought, it's a sum of power laws

$\sigma =C_1 x^{-a} - C_2 \frac{1}{(a+1)^2}x^{a+2} $
where $C_1$ and $C_2$ are positive constants and $C_1$ might formally be infinite.

Incidentally, modulo this weird pressure stuff, I have an idea that the gravitational force being constant  is the reason that the $a=-1/2$ scaling is the dividing line between power law initial conditions that go to the attractor solution, and those that don't (and indeed probably the reason that $a=-1/2$ is the attractor).  When $a$ is shallower than $-1/2$ the force increases with radius, whereas if it's steeper the force decreases with radius.  Somehow if the force is increasing with radius it has enough umph to shove the profile toward the attractor solution, but if it's decreasing with radius it can't push it in any farther (especially since the object is already more "collapsed" than an $x^{-1/2}$ object).  Anyhow, I am dinking around trying to find a way to say this thought with math rather than this annoying handwavy argument.  Jerry is too... he will probably beat me (if it can be done at all).

Alexia