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

Response to questions II

Hi, thanks again for your prompt responses!  I'll try and provide answers/plots here.

How did you define the phase angle you are plotting? Note that getting the phase angle of the instantaneous orbit in a frozen potential requires to obtain the map from (x,v) to (J,phi) (action-angle variables), though J is not required. Thus phi is not just some polar angle in (x,v) space.

I did not do such a transformation... the angle in the plots is defined as the arctan(v/x), so in this case it is just some polar angle in phase space.  Philosophically I understand that there is an arbitrary scaling in V, so that this angle is not well defined, but on the other hand, if you look at the plot I made of angle vs. x, the density at each x is the projection through all the layers, so the density does somehow depend on this non-meaningful angle.  Another way to say this:  if all the velocities were scaled, then the windup proceeds at a different rate as a function of time... so a given snapshot would have 1) different angles yes but also 2) a different total number of winds.  Therefore I am confused.

what do you think about this: 2010arXiv1010.3008V ?

Whilst this paper is rather confuse, I think in a way they also claim that
mixing (of energies for example) is not relevant or does hardly happen,
similar to our findings. However, I don't think they provide any explanation
for the final profile.

I read this paper, thanks for pointing it out.  I had not thought of Scott's objection, but I agree with him.  Specifically, I think the energy transfer between particles on various orbits depends a lot on their exact phase relationship as they bounce through the middle (or I guess periapse), and this is not evolving dynamically in their approach, which may make a difference in the final density distribution.

I am not sure I agree with your assessment that they are finding that the energies don't mix.  Maybe I am interpreting this wrong, but if you compare the 1st and 2nd panels of figure six, the different colored curves seem to me to be roughly proxies for initial and final energies (for a snapshot at the final time).   The first panel shows that while there is no mixing in the outskirts, it's definitely particles that started life in the middle that contribute most to the inner radii.  We are finding this too.   The fact that in the second panel the ordering of the colored curves (even in the outskirts) is not the same implies that the turnaround radius of particles at the end of the day does not seem to be a monotonic function of the initial radius.   Caveat: since I am thinking the initial potential energy as the dominating energy, it's possible that whatever they used for initial velocities messes up this logic.

I think I understand what you mean by mean-field effect. The way I would say it is that the scaling behavior only shows up as long as the system looks cold, i.e. as long as there are well defined streams in phase space.

Well, I think the surprise is that the scaling behavior DOES show up while the system is still cold and there are well defined streams in phase space...   and apparently stays that way after the well defined stream disappear.

I think the ISW effect is not a useful analogy for non-cosmologists (!)

Agreed, I won't use it in the paper  :)

The situation in three dimensions is that rho/sigma^eps is a power-law over a much larger range than either rho or sigma, if and only if eps=3.  Various people claim that this is because eps=3 makes rho/sigma^eps look like a phase space density, but this is a weak argument. It's not at all obvious what value of eps works best in one dimension but if its not eps=1 then we have a counter-example to the arguments about phase space.  

Here is a plot of $\rho/\sigma^{eps}$.  The curves have been offset by a factor of 100, so their curvature at high r is easier to compare.  Eps=1 is indeed not the exponent that maximizes radial range of the power law, but I am not sure that is significant, as I will explain below:

You'll note that eps=3,4,8 all seem to give equally good power law behavior, and there is a reason for this.  Jerry pointed out to me that the behavior of particles in the outskirts of the system look to be behaving isothermally, and if that is the case, $\sigma$ will go toward a constant.  Here is a plot of $\rho$ and $\sigma$ individually, and I have also plotted $r^{-1/2}$ truncated with something that looks like $e^{-\Phi}$, which indeed seems to fit well (in a chi by eye sense).   

The velocity dispersion $\sigma$ seems to go to a constant (if you ignore the last data point), and is furthermore fairly noisy.  If you raise a constant (with noise) to a high power, it amplifies the noise. Dividing by this noisy thing makes it harder to detect any rolling off of Q, but if sigma is actually FLAT all powers of sigma are doing the same thing to the curvature of Q... nothing save hiding it in the noise. 

Of course $\sigma$ is rolling slowly toward flat, and may be rolling in concert with $\rho$ to keep Q scaling the same way in the rollover region... we might be able to argue from the plot above that $eps=3$ is a legitimately a power law farther out than $eps=1$.  But I am not sure.

You say that you have plots of the contribution of each wind to the density and the energy. A related question is this: is the potential gradient near the center dominated by the "small" winds near the center of phase space, or the "large" winds?

This question is confusing, can you clarify?  The force from particles outside r is canceling, so only the potential due to particles inside r influence the dynamics of a particle at r, or am I misunderstanding something?  Sorry...

Response to questions

Hi, Here are some responses to questions you both sent.  Scott's in red and Walter's in blue.

I'm not sure I understand what you mean when you say that the scaling is not a mean-field effect...

Mean-field is probably not a good way to describe what I was trying to say (since that implies some kind of averaging).   I was thinking along these lines in using the term mean field:  Once the system winds so tightly that there are order 1ish particles per wind, then the coherence is lost, and the impact of the various winds are somehow "averaged" or "mixed" together.

What I was trying to say is that the scaling behavior shows up well before the coherence disappears.  In the energy plots, we saw there was a short period of violent relaxation followed by a long equilibrium followed by mixing and the gradual deterioration of the well ordered nature of the particle energies.  I am saying that the $r^{-1/2}$ scaling shows up in the earliest phase of the evolution of this system during the violent relaxation, so this is the portion that we need to understand better analytically.  Possibly understanding functionally how the energy is redistributed along the wind will be key to a better understanding.

Am I correct that you are finding that the particle with the k^{th} largest energy at t=0 still has the k^{th}largest energy at the end of the simulation (or close to it)?

Yes, this appears to be the case.  Interestingly, the same is true for the interior particles, its the particles in the middle energies that tend to get mixed up if the simulation is run too long.   Incientally, The energy gaps start to open pretty early, you can see the tracks diverging even as low as 50ATUs.   Also, mixing among the ordering of the energies of the middle particles correlates with larger fluctuation on the energies of the interior particles.  It may be that the phase relationship between the middle and interior particles has cycled all the way around until it has reached pi/4 again.  

How are you defining energy? (1/2)v^2 + the sum of the potential
energies for that particle?

Yes.

It would be interesting to plot initial rank in energy vs. final rank in
energy. This would show the effect quite dramatically. 

Sure, here is a plot for the same six timesteps.  This is for the particles that started life with initial position > 0.

You can't really see it above, but departures from monotonicity begin early, especially near the center.  Here are the same six snapshots, plotting final rank-initial rank.

Actually, I thought it was interesting to look at a few other times between 40 and 200.   I don't know anything about resonances, but these departures don't exactly look randomly placed. 

I think that's it for Scott's questions, now on to Walter's questions.


what is the "ISW effect"?

I was referring to the integrated Saches-Wolf effect, that we see in the CMB anisotropies.  Basically, when photons fall into gravitational potential wells they blueshift, and when they climb out again they redshift.  In our universe, potential wells on the very largest scales are decaying at late times, and more so if there is more dark energy.  Photons from the CMB traverse these decaying wells on their journey to us, and get a boost in energy because they don't redshift as much as they blueshifted.  We see this as slightly elevated power in the CMB anisotropy spectrum on very large scales. 

On p13, sigma and rho are power laws only out to about r~1 Is Q a power law over a significantly larger range in r? If so, for which exponent eps in Q=rho/sigma^eps?
 

I will look into this question further...  just to clarify, in the plot I mistakenly plotted $\rho/\sigma^3$.  Obviously when both $\rho$ and $\sigma$ are power laws, then Q will be for any eps, but you could be right, there may be a special eps for which the departures in power law in $\rho$ and $\sigma$ cancel each other out and Q remains powerlaw farther out in r.  I'll try to find out.   Actuallt eps=3 is looking pretty good!

on p15, another effect is the linear density along the "winds".  The winding stretches the initial curve and the linear density along the stretched curve is reduced, the more so the more it is wound up. Thus the inner "winds" contribute each less to the density. I wonder whether one can not quantify this better (you say "hard to model" on p21), via mass conservation along the phase curve.

You touch on a point that I have been struggling with for weeks.  I really wish I had a variable for the linear density along the winds.  Every time I try and go down that route, Jerry is telling me that there is no metric in phase space, so the length of the arc has no meaning, therefore I can't divide the number of particles by the arclength as I would dearly love to do.   This is actually how I discovered the energy thing, I was looking for a monotonic quantity along the arc, and hoping that the particles would sample this quantity in a predictable way.  I will think about mass conservation... I actually think that mass scrunches up along the arc, sampling region of high velocity and small position more densely and caustics more sparsely.

I have plots of the contribution of each wind to the density (and the energy) but they don't help me much, after a while each wind is contributing 1 or 0 particles to a particular radial bin, so I can only reliably trace the contribution of the few outermost winds.  Here is an example plot:

Here each "layer" is a wind where the layer level is switched as the phase angle goes through zero.   Red is outermost.  You might worry that yellow and pink are overlapping, but I think it's just an artifact of finite sampling, you'll note that there are just a few particles contributing to each bin in that region.   I wish I'd kept initial positions for the runs that were 80k particles, but I hadn't implemented it yet.  I guess I could run another, it's only 4 days.  I'll think about that.   I have these layer plots for energy too, but they are confusing.

As phase-space has no well defined metric, your definition of the phase angle is somewhat arbitrary. One may use the phase angle of the instantaneous orbit (the orbit of the particle if the potential were frozen at its current form).

Indeed, a fact that frustrates me daily.   My interpretation was that the phase angle I plot means exactly that, the phase angle of the instantaneous orbit.  Have I goofed that up?  I got that plot by accident, because I was using the phase angle to define the winds.   There is an even weirder plot that I made, that I totally don't understand but it could be significant.  Take a look at this:

This shows the final phase angle as a function of the initial positions (which for a while are a proxy for energy except in the innermost particles, where they start to mix up in energies).   You see the winds building up, but way before the point where the winds would be poorly sampled by the finite number of particles, you see the inner ones falling off the bandwagon, and hanging around at small position and high velocity.  (Red and blue are particles left and right of center.  T=23 should have been recentered or else was recentered poorly, sorry).   Anyhow, I have no idea what this plot means.

The "mixing of energies" seen at late times (p30) may be due to loss of numerical resolution: instead of the coherently time-varying potential of the winding spiral, we have inidividual particles scattering about.

Possibly, certainly at late times that seems likely.  The plots scott asked for, however, indicate that energy mixing actually commences well before coherant spiraling disappears. 

There is some funny effect (nicely visible on p33 at t=160) when the winding up is perturbed and a large gap appears. I also found this myself and have not yet found the reason for this behaviour (it occurs for both of my numerical algorithms and does not disappear when using "better" initial conditions). Perhaps there is some numerical instability in my code? Anyway, this effect seems to coincide with the onset of the energy mixing.

Yes, we wondered about the gaps.  I think it would be nice to answer whether this is an artifact or real.

Okay, that's it for now.  Long post.  Sorry.

The talk slides and a few extra things

Hi Scott and Walter,

    Sorry to decouple quite a bit in the last month or two.  I've done some work on our project since I last posted, there are a number of new plots.  I think the easiest way to catch up might be if you look through the talk I gave at our informal seminar here.   Here is a link:

http://public.lanl.gov/aschulz/talks/LA-astro23Feb.pdf

I had intended to send it to you for comments before giving the talk, but inevitably life intervened and I was working on it until the late the night before.  I think the main new thing I've found is that the r^-1/2 scaling is not really a mean field thing that arises after the system equilibrates, in a sense it is there from the beginning, which I suspect was not known before. 

I've recognized a few mistakes in the talk since giving it, Q in 1D should likely be rho/sigma,  just on dimensional grounds, and I don't think the equations I wrote down for N(w) make any kind of sense.  It was late at night.  

Also, I've added a bunch of slides after the conclusions, dealing with Energy evolution of the system.  These were not in the talk, but I thought you would like to see them.  The main point is that the particles behave like beads on a string, no matter how long the string gets, the beads do not pass one another, in that they keep their initial ordering along the wind.  Energy is a proxy for this ordering, because it remains a monotonic function of the initial position for an extremely long time.  This is significant because unlike the x and v, it evolves very slowly in time.  In a sense it governs how the "beads" are sampling different portions of the "string."   Energy might therefore be a good way to quantify how particles on each branch of the winding are contributing to the density at any given position.  I am still trying to work out how to do it, however.  One problem is that I don't have an expression for the potential energy, so my initial plan of solving to eliminate v is not working so well. 

Uniform density runs, sorted on initial position

I went back to the original 4 cases where I varied the initial conditions.  The initial density profile is uniform, and the particles lie exactly on a grid.  Each run has 10k particles, and I have repeated it 4 times to get error bars on the final density profile.  Here is a plot to remind you of what the different initial conditions were:

Runs 1 and 3 behaved very similarly in their evolution, forming a single object that winds up and collapses into the final virialized object.  Runs 2 (with the inflection point in the initial velocities) and 4 (which was initially expanding rather than contracting) were similar to one another in that they each formed 2 distinct objects that later merged to form the final product.  What you will see in this post is that the fate of the innermost and outermost particles depends greatly on the formation history, although the final density profile of the remnant does not. 

Here are the final phase space configurations (transformed so you can see the interior parts too) of the 4 runs.  Reddish brown particles started farthest out, and blue and green particles began life farthest in:

Here we see that for initial conditions 2 and 4, the ones in which there was a merger event, the innermost particles are driven to the far exterior of the object, and the innermost particles of the remnant started life in the outermost ring.   In contrast, for the monolithic cases, the innermost particles of the remnant come from a mixture of intermediate and inner rings, and just as we saw in the power law cases, the distribution of the inner particles is not isotropic in phase space.

Here are plots the initial vs. final positions and velocities of the particles:

These show that there is a very specific band of the initial particles that wind up in the middle of the remnant, in the merger cases.  Also, the "tail" that we were seeing int he velocity plots of the power law cases does not seem to appear here.

Finally, just to complete the story, here are plots of the final density profile of the 4 different runs.

 

Here it is evident that the monolithic cases scale as $r^{-1/2}$ over more decades than the cases with merger events, but even the cases where the inner particles come mostly from the outer layers scale quite convincingly (at least by eye) like $r^{-1/2}$.

IP sorting for all the power law runs

Got a little time during nap today to make more plots.  Here are plots comparing the fate of particles in different initial positions for four different runs:  $\rho \propto x^{-1/2}$ (you saw these plots yesterday), $\rho \propto x^{-0.8}$, $\rho \propto x^{-0.3}$  and $\rho \propto x^{+0.5}$.  I made the log and polar versions of the plots, but not the zoom ins.  Here they are, labels are at the top.

 

 

 

General trends:

1)  The steeper the initial power law slope, the farther out the origin of the innermost particles at the final configuration.
2)  The steeper the initial power law slope, the more asymmetric the final configuration is in phase space.
3)  If the density of the initial object increases with radius rather than decreasing, it is the outermost particles that wind up in the middle, and the innermost wind up around the ouside in phase space.
4)  For runs with a negative power law slope, although they maintain their density scaling in the intermediate range, it seems like almost none of the particles collapse in farther than the position of the innermost particle at the beginning of the run (only in the very shallow case do particles make it farther in).

Sorting by inital positions

I have added a feature to Walter's code, namely each particle carries about a memory of where it started it's life.   I have done several runs, but all of the following plots refer to initial conditions that were a power law density profile with $\rho \propto x^{-1/2}$.   I have defined 5 logarithmic bins in the particle initial positions.  Below is the phase space diagram, with different colors representing particles that started in different initial positions, the reddish brown ones began life the bin farthest out, and the dark blue began life in the innermost bin.  Each plot is zoomed in from it's predecessor (you can click on the plots to view them larger):

Just for illustration, here are the same types of plots, but for the initial conditions (note the aspect ratio is no longer square):

To get a better idea of how the different colors are distributed around phase space, I did a transformation on these data to fit all the information onto a single plot.  Below I have defined a distance in phase space which is $r=\sqrt{x^2 + v^2}$.  Then instead of plotting $(r,\theta)$, I plot $(r^{0.1}, \theta)$ so that we can see the inner parts. 

Here are the initial conditions in this view:

These plots demonstrate that the outermost bins retain their identity, but the inner bins get mixed in the ordering of their positions.  All of the velocities increase, but the inner bins tend to have only high velocities and small positions.  The phase space orbits for the outer particles and inner particles seem qualitatively different.  This result is robust to the phase of the winding:  I checked that it holds true for the last several time dumps, and also for different runs of the same initial conditions (different runs have different numbers of particles).

Another useful view is to look at scatter plots of $|x_i|$ versus $|x_f|$ and $|v_i|$ versus $|v_f|$.  Here they are, the black lines show $x_f=x_i$ and the analogous for the velocity.

You can see that the positions on average are collapsing inward, as expected, but the innermost particles do not collapse inward as much as some of their outer counterparts, in fact the very innermost ones seem to be thrown out of the central region.  The velocities are all increasing, with the innermost particles seeing the largest increase in velocity.  This basic structure is set up very early in the run, here is the first time dump after the start of the run, at 20 ATUs.

I was puzzled about that tail in the outer particles, I thought maybe looking at these plots in a slightly different way might shed some light.  Here are $x_i$ versus $x_f$ and $v_i$ versus $v_f$ plots, where again I just zoom in rather than doing the polar trick.  The axes have the same ranges as in the beginning of this post (except v_initial, which I scale down in the same proportions):

So you can see the velocity tail in a different form, in the first of these, I am still thinking aobut why it's there.  The innermost bins are doing something qualitatively different than the outer three, in my opinion.

There are still many more plots to make, I will start checking some of the other power laws.  Also, if I want to do warm runs, I need to repeat the modification to Walter's code, adding one more field (and since it's the last available bit, that is it, can't add any more data to carry around without doing something hard).    I will try and make some more plots this weekend, have good one.