Tuesday, April 12, 2011

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. 

Friday, April 8, 2011

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