# LAMMPS absorbing state

h = 1.47:


###### tags: `Granular Simulations` `granular`
### [27/02/23]: Amp vs Height phase diagram

Simulations are run at 53Hz.

Height and amplitudes are in particle diameter.

We will focus on the region where particles always touch the ceiling.

The first simulation was done at $N = 5000$ and $\phi = 0.2$. This is a bit too low:

#### $\phi = 0.2$


Here are the dynamical degree of synchro and %of particles touching the ceiling

We see three things:
* At high amplitude and low height, the particles touch too frequently the ceiling, thus the absorbing state.
* A sharp transition at high height and low amplitude as expected since we are in a synchronous zone.

* A potentially continuous transition in the zone of the asynchronized particles.


The discontinous transitions makes sense because we find the absorbing state at high amplitude since the number of collisions with the ceiling increases and so does the tangential drag.

But, the continuous one happens at high amplitude which is a bit strange... It might be an issue with the time to relax but it doesn't seem so (this is the evolution of the energy in a region where synchronization happen):


(the energy scale is different but it is because I'm an idiot, we can trust the tendency of the curve if not my poor choices of python scripting. Also the steps are not real steps..........).

Taking the log of the imshow gives us this:
Which means that what we observe above is not the continuous transition IN the async region but more the transition from the synchronized region to the asynchronized one.

Now, let's wait 2 days fpor the whole phgase diagram :)

### [27/02/23]: Amp vs Height vs $\phi$ phase diagram/ Giant fluctuation?

#### Not a real transition?

Indeed, the "transition we see" is really a transition from the synchro to the asynchro zone:


In the async regime, the system may be always active or the continuous transition might just be a smooth one since we have activity (very looow activity) even at low packing fraction.

_______

#### Hyper non uniformity/Giant fluctuation?

When we place ourselves here:

here is what we get, an homogeneous active state. We can see sometimes agregates forming (top/center left or bottom right) but they get immediately destroyed in an explosive manner even though we are in the sync zone and so discontinuous transition. Thus, we kind of observe cascade of energy even in this zone of the phase diagram where the transition is discontinuous but this is not as clear as in the async zone as we will see later.

Same thing if we place here:
Same thing but less explosive


Now, if we place ourselves here:

Nothing moves (off course) and we find this system with (I think...) huge fluctuations/hyper non uniformity: 


Probably because we are a tiny bit above the mean free path required to make the system active and so each particles just "almost touch" their neighbours but not enough so that the state could become active. I don't think we would see this thing experimentally since the absrobing state of thge experience is not a real one...

I tried to plot $\sigma^2/r^2$ according to $r$. It looks like that the system is going to 0 like a hyperuniform system (the opposite of what I expect...), but maybe we don't have enough particles to prob the $r \rightarrow \infty$:

Here is an active state in the async zone, the dynamics is way more separated into dead and active particle than what we observe in the sync zone:


Here is the velocity distribution (remember that we are probably close to the transition packing fraction but not exactly):

In absolute velocity we don't see two peaks as a sign of two populations:

Even for a system very very sparsely active, we only obtain this kind of curve. We might see two regimes but it's not clear.

___________________
Here is the full phase diagram:

|  | | |
| -------- | -------- | -------- |
|  |  |  |
|  |  |  |
|  |  |  |
|  |  |  |
|  |  |  |

We can't see correctly the continuous traznsition because the discontinuous one is too BRIGHT!

Here is the phase diagram without the upper (synchro) zone (but still with the lower part). No better:

|  |  |  |
| -------- | -------- | -------- |
|  |  |  |
| | | |
||||
||||
||||

More simulations are running at a low packing fraction I hope that we can see the transition here

### [02/03/23] Nothing works today!
* Tried to come early, the RER B putted me through hell.

* Tried to compute the structure factor, failed miserably for $q\rightarrow 0$ which is important to study correlation lengh.

* Tried to find two populations at the critical packing fraction, didn't work.

* Tried to find a transition at low density in the asynch zone, found this: $\phi = 0.1$

* Reran the experiment to see if the particles were still not touching the ceiling while being in the sync region and in an active state, found that actually all particles touch the ceiling, can't figure out why the first one didn't gave this result (particles have a normalized radius of 1):

*
| New experiment (touching the ceiling) | Old experiment (not touching the ceiling) |
| -------- | -------- |
| | |

* At least, now, we know that we don't have to care about this possible issue anymore. Particles that touch the ceiling in the dead state also touch it in the active state.

* Moreover, particles won't synchronized if there is no roof.

### [03/02/23]: Video, packing fraction plot and critical packing fraction probing

POrobably not extremely interesting, but at high density, at a given amplitude, the energy of the systems might converge to the same point.

| One point of convergence for this range of h | Two points of convergence for this range of h|
| -------- | -------- |
|  |  |

For the second picture, we are always in the async zone thus, the change of behaviour is not caused by such effects.

We can see on the highest packing fraction imshow that there is some kind of bifurcation going on:


The zone with dash corresponds to a region in which the final energy is the same.

ANYWAY!

__________

I tried to see if we could obtain a real continuous phase transition (by that I mean: $E(A_c) = 0$ with $A_c$ inside the async zone and $E(A_c^-)>0$) in the asynch zone but couldn't:

________

#### Discontinuous A = 0.06225, h = 1.95: BIG JUMP at the transition


| Before $\phi_c$ | After $\phi_c$ |
| -------- | -------- |
|  |  |

At $\phi_c⁺$

Close to this point, we also found (~ h = 1.98, A = 0.06):

| Below $\phi_c$ | above $\phi_c$ |
| -------- | -------- |
|  |  |
Where we can see a lot more these agregates (even in the active state)

_________________________________________

#### Discontinuous A = 0.07, h = 1.76: MID JUMP at the transition


| Before $\phi_c$ | After $\phi_c$ |
| -------- | -------- |
|  |  |


______________________________________
__________________________________________
__________________________________

#### Continuous A = 0.0715, h = 1.41

Continuous:

in log scale (not going to 0 because of the simulation time I hope, with this continuous transition we often see a tiny bit of activity below what I would call $\phi_c$):


| Before (what I would call) $\phi_c$ (here $\phi = 0.17$). This is a gif| After $\phi_c$ |
| -------- | -------- |
|  | |

(sorry the videos suck, I messed up windows of observation)
Might be hyperuniform (or not)


_______

| Continuous | Discontinuous |
| ----------------- | ---------------- |
| In the async zone | In the sync zone |
| Definitevly a massive increase of time scales at $\phi_c$ | Probably not a divergence of time scales at $\phi_c$. |
| Kind of uniform (maybe hyperuniform) structure at $\phi_c^-$ (particles are organizing in the max distance to neighbour?)|A lot of aggregates close to $\phi_c$ (below and above). Particles "almost touch" their neighbour below $\phi_c$ but stop before hitting one|
| Is currently not found while varying the amp| Is found while varying the amp|
|$P(\textrm{abs}(\vec v)\rightarrow 0)\neq0$ (really a population of dead particle)|$P(\textrm{abs}(\vec v)\rightarrow 0)=0$ (like a usual MB) |

### [06/03/23]: Kuramoto vs variance
So far we were discriminating the synchronization of the system using:

$$\displaystyle{S = \dfrac{\textrm{var} (z_{cm})}{1/N\sum_{i=1}^{N}\textrm{var} (z)}}$$

Where the variance is taken over different snapshot at the steady state.

A somewhat better synchronization parameter would be an equivalent of the orientation bond parameter: we assign to each particle a value $\theta_z$ linked to its position between 0 and $2\pi$ then define:

$$S' = \left| \dfrac{1}{N}\sum_{j = 1}^N e^{i\theta_z}\right|$$

if all particles have almost the same value of $\theta_s$, $S'$ goes to 1. Otherwise, it goes to 0. $\theta_z = z$ would be a good choice but we'd prefere some parameter that could differentiate between particles going down or up at the same $z$. A simple choice is:
$$\theta_z = \underbrace{\pi\dfrac{z - \textrm{min}(z)_t}{\textrm{max}_t(z - \textrm{min}(z)_t)}}_{\textrm{between 0 and }\pi}\times\textrm{sign}(v_z)$$

which goes from $-\pi$ to $\pi$.

$$S' = \left| \dfrac{1}{N}\sum_{j = 1}^N e^{i\theta_z}\right| ~~\textrm{ with }~~\theta_z =\pi\dfrac{z + \textrm{amp}}{\textrm{h + amp}}\textrm{sign}(v_z)$$

which goes from $-\pi$ to $\pi$.

Here is the difference:

| Continuous | Discontinuous |
| -------- | -------- |
|  |  |

For the continuous transition, no difference. For the Disconcitnuous one, I don't know which one is the correct one. If the sytncghronbization really goes to 0 or go up again.

The kuramoto order parameter is nice, because we don't need multiple snapshots since we don't use a statistical parameter like the variance.

### [08/03/23] Critical packing as a function of amp and height


bump, etc, simulations running for more

### [09/03/23] Explaination of the agregates

For the continuous phase transition, the transition is driven by $\gamma$, the friction with the roof. But for the discontinuous transition, the transition is driven by $\Delta_s = 0$. Which means that the absorbing state is caused by dissipative collisions (and $\gamma$ but mostly dissipative collisions between synchronized collisions). The aggregates we observe are probably the signature of an inomogenehous free cooling that can also be observed with the synchro delta model as long as $\alpha$ is low!

### [13/03/23]

As found by Andrea:
h = 1.76 A = 0.088
This transition happens in the transition zone between synchronized and async particles
We observe something very similar with what we find with **the theory** with $\Delta_s \neq 0$. That is, the first half of the curve (before $\phi = 0.33$) follows a continuous phase transition and the second half continues as a sharp transition (discontinuous one that couldn't be). Note however that this kind of thing hasn't been observed in the SIMULATIONS but only as a false prediction of ouyr theory (for now).
Here is an example with the delta model:

_________

| h = 1.55 |
| -------- |
|  |
|  |

### [17/03/23]

Z energy of the active or dead state (they are the same).

SO, in fact, the energy changes with h. We can even see the transition from async to sync.
I also checked since I had doubt, the async state is indeed a stationnary state/ particles won't synchronized after a long time.

Also at the interface between the async et sync zone we observe a verrry large synchronization time:

So, our plan of using one $\Delta$ for every $h$ falls short. But still, we can approximate the collisions in lammps by a simple collision rule:

$$\vec v' = \vec v + \dfrac{1+\alpha(|\vec {\delta v}|)}{2}(\vec{\delta v0} \cdot \vec\sigma)\vec \sigma$$

Let's place ourselves in 2D with an x and z axis, integrate the z axis and forget about the z axis:

$$v'_x = v_x + \dfrac{1+\alpha(|\vec {\delta v}|)}{2}\left[\delta v_x\sigma_x + \int \delta v_z\sigma_zP(\sigma_z)P(\delta v_z) d\sigma_zd(\delta v_z)\right]\sigma_x$$

For the angle distribution, it's I think easy. Since $|\vec \sigma| = 1$ we have $\sigma_z = \sin(\theta)$ with $\theta$ the angle between $x$ axis and $\vec \sigma$. Moreover, we have:
$$P(\theta)d\theta=P(\sigma_z)d\sigma_z\Rightarrow P(\sigma_z) = P(\theta)/\cos(\theta)=P(\arcsin(\sigma_z))/(\cos(\arcsin(\sigma_z))) = P(\arcsin(\sigma_z))/\sqrt{1 - \sigma_z^2}$$

Let's supoose that all angle of collision are equally probable, which is definitively not the case for $h\simeq 2 rad$. We find:

$$P(\sigma_z) = \dfrac{1}{\pi\sqrt{1 - \sigma_z^2}}$$

Let's make this harder! Now, consider than the particles still collide with a random uniform $z$ but the angle of collision is constrained by the height of the box like this:


Then the max angle is $\pi/2 - \theta$ and the min one is $-\pi/2 + \theta$. With $\theta = \arccos(h - 1)$ (dont forget that our h in the simulation is expressed in units of diamter) for h between 1 and 2. Sanity check:
* If h = 1, particles can only collide frontally and indeed the formula gives us $\theta = \pi/2$. If $h = 2$ particles collide in every possible direction, and we get $\theta= 0$.

This lets us correct the formula above:
$$P(\sigma_z) = \dfrac{1}{(\pi - 2\arccos(h - 1))\sqrt{1 - \sigma_z^2}}$$

With $\sigma_z\in [-1 + \sin(\arccos(h-1)), 1 - \sin(\arccos(h-1))]$, we could develop the sin(arcos) ...
This looks extremely suspect...

_________


orange is h = 1.5 and blue is h = 1.8, A = 0.08


We might want to have a distribution of velocity that depends on $z$...

### [27/03/23]

With lammps, as in the simple delta model. The time scale doesn't diverge at the transition point with $N\rightarrow\infty$. Here we are at the transition point of a cointinuous transition:

The dynamics is very slow compared to when we are far away from the transition but it is independent on N! For the first order phase transition, we observe the same thing.

A nice thing here, is that we observe two effects. First, the massive decrease which is caused by dissipative collisions, then the slow decrease which is caused by the friction with the plate.

Also, observing a power law behavior (at the transition) would be kind of stran,ge since the only mechanism that could lead to this behavior is the free cooling ($E\sim 1/t^2$) which is in part hidden by the viscous drag ($E \sim e^{-t}).

### [30/03/23]
In blue and orange is plotted the Z energy of dead particles in the absorbing state

| Synchro | Asynchro |
| -------- | -------- |
| | 

For the synchro we see particles bumping into the ceil and the floor and consistenly gaining energy.
For the asynchro, we see that sometimes particles can't reach the roof (my first analysis was wrong, I was not taking into account the fact that the plate was moving. Still, they are "more" touching the ceiling than the particles which had a position in the phase diagram which I cutted because I said they were not touching the ceiling) Since these ones couldn't even reach z = h - radius. While in our current async zone, the particles reach z = h - radius everytime but maybe not z = h + A - radius which is the max z a particle could reach). Which induces a delay into the position of the particles compared to others that could touch the ceiling and thus cause asynchronization.

### [20/04/23] and [24/04/23]: why are populations coarsening? Single particle evolution and LAMMPS comparison with EDMD

#### Doubly synchronized state

What we thought was an async state is in fact a state in which we observe two sync populations of particles (both doing the same motion but phase shifted between population, see below).

Here is a snapshot of a system which went to the absorbing state close to the transition point. Color is $z$ position.

We can see that there is a correlation in opulation belonging. Clusters of populations red are clearly visible for example. This is due to the fact that collisions between particles of group red and blue add energy into the system while collisions between particles belonging to the same population will be purely dissipative (since they are syncronized). Then when some populations start forming, particles in the same populations start feeling an effective attraction. Here, the formation of these clusters was interupted by the system going to the absorbing state due to **finite size effect**. But, at what we would call $\phi_c$ without any synchronization effect, when time scale diverges. The system would have as much time as it wants to reach a state in which we have a huge blob of population 1 and a huge blob of population 2. The energy decreasing over time because collisions, in average are lesser and lesser effective to add energy into the system. At this point, every particles (except at the boundary between the two populations) will be in sync with their neighbour and the transition will be first order.
This assumes that the number of active particles goes to 0 as $\phi\rightarrow\phi_c$ in the situation for which we forget about the sync effect. This also assumes that the power law behaviour of the system won't destroy this argupment based on local arguments which might not hold with such long range interaction (which might be smeared out by the synchronization phenomena since the transition would juste be a failed continuous one).

When the synchronization time (which is large) is greater than the time needed to relax to the steady state. The system is going for a second order transition. But as soon as we get closer with $\phi$ to $\phi_c$ so that the time to reach the steady state get greater than the time to synchronize which is more or less equivalent to the point at which we see the number of active particles decreasing (since the particles need to be dead to synchronize in one or the other population); the system will probably coarsen into two populations and probably hide the second order transition behind this phenomena.

##### note to Andrea's remark
I think it's normal to not observe the system clusterize in density at $\phi = \phi_c^-(L\rightarrow \infty)$ as in the sync zone where the transition is async and the dissipative collisions, close to the critical point, make the system behave as a cooling granular gas.
In the doubly sync zone, the system is purely purely driven by dissipative collision, very close to the critical point. And on top of that, only at $\phi \rightarrow \phi_c^+$ since below $\phi_c$ the system is driven by collisions with the roof.
In fact, above the critical packing fraction we can already see similar clustering effects (the effect would be greater if I tried to go closer to $\phi_c^+$) while still dying but with a L large so that the population can coarsen:

| Color by population | Same snapshot where we can see density clusters |
| -------- | -------- |
|  |  |

#### Single particle trajectory

Frank, as always, was not wrong, in the doubly synchro region, two populations of particles oscillate with period equal to twice the one of the plate (the peak at 1/2 in frequency domain). And the population one is phase shifted with the second one:
Orange is population 1 and Blue is population 2.

| Time domain | frequency domain|
| -------- | -------- |
|  |  |

I can get back the same type of behaviour with the simple model (penetration of plate by a strping force and viscous drag):
**h = 30./10.
amp = 3./11.
w = 1.
viscous = 1.
k = 1001.
g = -0.1**
But this behaviour is really not standard, most of the time we observe synchronization (with one or two bump with the bottom plate for one hit with the top plate) or a purely chaotic behaviour:
**h = 20./10.
amp = 11./11.
.
.
.**

#### Comparison between LAMMPS and EDMD
LAMMPS simulations were done with 10k Particles

| LAMMPS | EDMD |
| ---------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------------------------ |
|  Not enough run, not enough time to equilibrate, not enough particles, phi interval too large |  |
| |  very resilient to treshold change! We are definitively not in the same kind of interval for EDMD and LAMPPS|
|  Fluctuations go down with phi ????? As in [this](https://arxiv.org/pdf/2105.15027.pdf#page=5&zoom=100,422,489) |  |
|  |  |

Longer simulations are running, but they are taking a lot of time :( .

# THERE IS PROBABLY AN ASYNC STATE, IM AN IDIOT

I'm plotting the histograms of positions (THANKS ANDREA) at some time for a pair of parameter.


Deep in the async region, I was not in the black paert but a bit above where we can have two populations or 4 (really reminiscent of chaotic map).
Maybe we haven't reach equilibrium... AGAIN... But the simulations lasted a long time.

### LAMMPS IN ASYNC ZONE VS DELTA MODEL

| LAMMPS IN ASYNC | EDMD |
| -------- | -------- |
|  Vanishing of variance at $\phi_c$. For directed percolation, we find that the variance just jumps (critical exponent equal to 0) in the mean field approx |  Divergence of variance at $\phi_c$ |
|||

I'm doing other simulations. I think that we could also observe divergences of variance. But with longer time since we see a lot of variation for the varance of the lammps simlulation but the window of simulation might be too small to probe it.

I'm running larger lammps simulations to check things. And also on a larger phi range.

### [04/05/23]
25k particles both
| LAMMPS | EDMD |
| -------- | -------- |
|bad statistics even with 60Go of dump :( ||
|||
|| longer time of observation |
|||
|||

___________

LAMMPS:
{%youtube jgH5psfT0P8 %}

EDMD:
{%youtube aFgIEymH6Bk %}

EDMD accelerated:

I haven't reach eq for LAMMPS close to the critical packing fraction...

### [05/05/23]
EDMD average over 300 realizations

### [24/05/23] Data Collapse (Andrea's incursion)

We try to study the behaviour of $\phi_c$ as a function of A and h. Before a quick summary of heatmaps:

SynchParam
I will probably cut the two lowest height

SynchRate (i.e. inverse of synch time fitted from time evlution of SynchParam, set to zero for simulation in which the final SynchParam is lower than a threshold)
Here threshold is 0.98

CriticalPhi

Non monotonous behavior of $\phi_c$ as a function of amplitude for different h

Sort of Monotonous behavior of $\phi_c$ as a function of amplitude for different h
Vetically ordered by h

We can also study $\phi_c$ in the asynch region:
all points colored by Synch parameter
averages and errorbars

TO DO: check with dissipation rate. Possible scenario: within a signle h, in the synch zone, all is governed by SynchRate. Can we explain different behaviour for different h with the fact that for larger h particles dissipate less energy? We could attempt a data collapse with the mfp argument

### [11/07/2023] Data Collapse II (Andreas's incursion)

In the following I skip the first two height where strange things happen (1.2, 1.26).

I tried to do the data collapse with the dissipation rate fitted from the exponential part of the energy decay but is giving strange results. It rescale well a two groups of curve but not all togheter (not really visible by next figure). Maybe there is something physical maybe it is just an artifact of the automatic fit script.

Anyway I think there is a better eay to estimate the dissipation rate from h and the average kinetic energy on z $E_z$. Indeed instead of fitting the energy decay in a particular regime, we can just say that the dissipation rate is related to the amount of collision between the floor and the ceiling sp it will be something like $\gamma \propto \sqrt{E_z}/h$. Then $\Delta$ can also be taken proportional to $\sqrt{E_z}$ since the energy converted into $xy$ is a fraction of the vertical one. A this point we can estimate the critical packing fraction from the naive mean free path argument $\phi_c^{mfp}=\pi(\sigma_L/2)^2\gamma^2 / \Delta^2 \sim \pi(\sigma_L/2)^2/h^2$ so a good idea is to rescale the measured $\phi_c$ with $\phi_c^{mfp}\sim \pi(\sigma_L/2)^2/h^2$. This is what we obtain
$phi_c$ VS amp
$phi_c$ VS synch freq
Rescaling with estimated $\phi_c^{mfp}$
It is not really a data collapse but because I think that between low and large $h$ there are two different functional forms.

Plotting only curves at high $h$ makes all more understandable


For the paper maybe I'll use only the fisrt and the last of these three...

### [24/05/23] LAMMPS/EDMD comparison with finite size effects

Energy and fluctuation size as a function of $\phi$
| | EDMD | LAMMPS |
| -------- | -------- | -------- |
| DISCONTINUOUS |  | |
| CONTINUOUS | |  |

Size of the gap as a function of number of particles

| | EDMD | LAMMPS |
| -------- | -------- | -------- |
| DISCONTINUOUS |  ||
| CONTINUOUS |  | |

_______

energy vs $\phi$ in log log with best phi_c I found.

| | EDMD | LAMMPS |
| -------- | -------- | -------- |
| DISCONTINUOUS |  |doesn't make sense...|
| CONTINUOUS |  | |

_______

Fluctuation vs $\phi$ in log log with best phi_c I found. Bad stat for the continuous one because I'm an idiot :')

| | EDMD |
| ------------- | --------------------------------------------- |
| DISCONTINUOUS |  |
| CONTINUOUS | |

ARRIVAL TO THE ABSORBING STATE AT THE CRITICAL PACKING FRACTION

EDMD
{%youtube 6pE4F7OpRco %}

LAMMPS
{%youtube wEX2yUNsaGc %}

Note the inhomogeneous density field at the end for LAMMPS because the effective coefficient of restituton must be pretty small and we can effectively obtain an inhomogeneous cooling state:


### [29/06/23]: Logistic map and chaos

With lammps:


I AM AN IDIOT: HERE f = 1 NOT 53
With model where particles can penetrate the plates:

With instantaneous collisions with the plate:


increasing dissipation by going down in the array:

| Hard plate | soft plate |
| -------- | -------- |
|  |  |
|  |  |
|  |  |

END OF f BEING WRITTEN AS 53 EVEN THOUGH IT S 1

h = 1.6, w = 53 (so like in lammps) hard plates:

WITHOUT GRAVITY:

At $g = 0$, assuming that the plates are not moving (but still giving velocity in a sin way) we find an equation for $v(t)$ that looks like the standard map equation. Explaining why the system is chaotic.

### [30/06/23] Math

#### Assuming plate is not moving and g = 0

g = 0, plates positions are 0 and h, the plates is assumed to give a velocity $Aw\cos(wt)$ to the particles. $t_i$ is the time of the $i$th collision and $v_{i\rightarrow i+1}$ is the velocity between the $ith$ and the $i+1th$ collision.

$$v_{0\rightarrow 1} = v_0~,~t_0 = 0$$
$$v_{1\rightarrow 2} = v_{0\rightarrow 1} - (1+\alpha)(v_{0\rightarrow 1} - A\omega\cos(\omega t_1))~,~t_1 = t_0 +|v_{0\rightarrow 1}|/z_0$$
$$v_{2\rightarrow 3} = v_{1\rightarrow 2} - (1+\alpha)(v_{1\rightarrow 2} - A\omega\cos(\omega t_2))~,~t_2 = t_1 + |v_{1\rightarrow 2}|/h$$
$$v_{3\rightarrow 4} = v_{2\rightarrow 3} - (1+\alpha)(v_{2\rightarrow 3} - A\omega\cos(\omega t_3))~,~t_3 = t_2 + |v_{2\rightarrow 3}|/h$$
$$...$$
$$v_{n\rightarrow n+1} = \alpha v_{n-1\rightarrow n} + A\omega (1+\alpha)\cos(\omega t_n)~,~t_n = t_{n - 1} + |v_{n-1\rightarrow n}|/h$$

This looks very much like the standard map 

The attractor I'm plotting above are position and velocity at $t = 2\pi m\omega$ with $m$ an integer. So there's still a bit more work to do.

We can also show that a particle boucing on a single plate sinusoidally driven will diaply a chaotic motion: https://vrs.amsi.org.au/wp-content/uploads/sites/92/2020/01/ceddia_julian_vrs-report.pdf

https://sci-hub.ru/10.1103/physreve.48.3988

https://iopscience.iop.org/article/10.1088/1742-6596/246/1/012003/pdf

It is left to show that if we had a plate above, the chaotic motion is still there... I guess that we can still show the (in)stability of the periodic motion in this case as it's done in the sci-hub.ru article.

### [17/07/23]: LAMMPS synchronization mechanism with respect to collision

As a proof of concept, collisions between synchronized particle can desynchronized them:



Nonetheless, this is a highly non physical example since for this simulation the collisions happen at a very high velocity (~ 20 times what would be the average velocity on an absorbing simulation). And the particles are launched without any angular momentum which might explain why the first collision causes almost no desynchronization while the second one does. What is clear, at least, is that, collisions between ''almost'' synchronized particle can easily desynchronize them (because of angular momentum or very small $dz$ that get amplified strongly by the collision). The tricky part is to estimate how much this can play a role in a realistic simulation and how much can this effect be taken into account into a simplfied model (it's clear that a S-S collision between almost synchronized and synchronizd would be vastly different in lammps because the small $dz$ difference might cause a massive change in the postcollision dynamics)

________

Acccording to what we saw earlier, two effects might play a role:
Angular momentum and small asynchronicity in the synchronized state, let's see what we observe with particles going slower (averge v of order of what we have in absorbing simulation) so that particles can bounce with a non zero x velocity before colliding and then gain angular momentum. We will all launch them at the same velocity to see if there's any important effect due to a very small asynchronicity.



Here, with small velocity, particles still desynchronize a bit, but synchronize back on a time scale way shorter than the synchronation time starting from a random configuration. We also see that each row has its own synchronization which means that the small asynchronicity difference at the beggining is not playing any role.

**This is definitively linked to the angular momentum**
When the two rows are separated by a large amount (so that we can get angular momentum from collisions with the plate before colliding between particles):

| | plot | gif |
| --- | -------- | -------- |
| Collision with the roof before collision between particles |  |  |
| No Collision with the roof before collision between particles |  |  |

When the two rows are separated by a small amount (so that no angular momentum is gained before the particle particle collision):

Almost no change

____

We can also look at the effect of collisions when particles are in the verge of synchronizing:


We see that collisions do not really asynchronized the system. It's just asynchronizing them on a short time scale. After this, the system resumes its synchronization almost as if nothing appened. There is not much change compared to a state in which we don't give any XY velocity to the particles and let them synchronize as they wish:


I think it tells us that the synchronization time $\tau_s$ is not only linked to the time needed to synchronize from a random state but also the time the system, being almost synchronized, need to synchronize back after a collision. This depends on the degree of synchronization of the two particles colliding (of course) but in a slighly surprising way!

If we make them collide sooner in their synchronization process we obtain this:

The fact that they collided didn't change a thing to their global evolution (WTF???). Maybe if they collided back after with an other particle. Or maybe it's because I'm in an idealized case but here the particles should have some angular momentum so :) :) :)

_________________

At least, if we make the following protocol:

- We let both row synchronize.
- At a time t' we randomize the z_velocity of the second row to desynchronize it and send the particles of row 1 onto them (still synchronize)
- During the whole time we measure the sync parameter of the 1 row. During 0 and t' it should increase. During t' and the collision time it should increase. And at the collision time, we expect it to go to a small value because the completely async particles should completely asynchronize the particles colliding with them.

Luckily, this is what we observe:

What we understood from this experiment and the last is that: synchronization can't really be modelled by a A and S state. But more so by a continuous state. Because when 2 particles collide it seems that they both end up in the least synchronize continuous state of the two particles colliding. For example, if 1 has sync parameter 0.8 and the other 0.7. they might both end up, with a sync parameter 0.7. BUT, we also saw that collision between two sync particles could, in theory aynchronize the system, but it looks like this can happen if there's a lot of collisions in a short amount of time since even though they might desynchronize shortly after the collision, they will synchronize shortly after that (if no other collision) with a time scale shorten than $\tau_s$ the time scale for synchronization starting from a completely radom state!!!!