# Dissipative EDMD
###### tags: `Dissipative self-assembly`

---

---

### Owners (the only one with the permission to edit the main test)

EF, AP, FS, GF

---

---

### Lack of FDT

First we note that with the euler scheme $E$ goes to $T$ only when $\Delta t_N$ goes to $0$ (for a binary mixture or a monodisperse system).

The euler method doesn't give the FDT (https://www.intlpress.com/site/pub/files/_fulltext/journals/cms/2017/0015/0004/CMS-2017-0015-0004-a013.pdf):


The runge kutta method does give it a bit more:

But I opted for the verlet integrator because it gives FTD whatever $\Delta t_n$ is (we can calulate $<v^2>$ with the formula below to see this):

$$v_i\rightarrow cv_i +\sqrt{\frac{T_b(1 - c^2)}{m_i}}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)~~~~~~~c = e^{-\gamma\Delta t_N}$$

This problem is then solved, even for interacting particles. Indeed, with the verlet algorithm:


This solves the issue with the $q_4-E$ correlation at equilibrium.

___

### E and $q_4$ correlation

We observe a correlation between $E$ and $q_4$ out of equilibrium:


The rearrangement of the particles causes the collision frequency to decrease, which is also observed at equilibrium of course (this is a geometrical propertie of the configuration):


If the total energy increases, it means that the dissipation is minimized since the energy injected into the system is constant over time. This reminds me the minimum dissipation principle of Onsager for systems close to equilibrium. Interestingly **since the maximisation of free energy at equilibrium "induces a low collision configuration" and that out of equilibrium the system tries to minimize its dissipation (and thus collisions), both equilibrium and out of equilibrium can give the same crystalline structure**.
> [name=gfoffi] perhaps an interesting paper: https://hal.archives-ouvertes.fr/hal-03184723/document

This reminds me the minimum dissipation principle of Onsager for systems close to equilibrium.
If this minimal dissipation principle really holds, we might engineer some crystalline phase by playing with $\alpha_{bb}$, $\alpha_{ss}$, etc

Unfortunately, it doesnt seem that this principle really holds. For example, we can start from a crystal formed at equilibrium and add dissipation in order to melt it. If the minimal dissipation principle was really holding in this case, we would observe an energy anti correlated with $q_4$ which is not the case:

We also observe this correlation for QC8 for example:


The result seems robust, in the following, I started from a crystallized configuration and then increased the size of the box in order to melt it. The size of the box is increased from the left to the right and the top to the bottom:

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

The system might be driven by the minimization of entropy production though... Which is harder to test I think.

Now, the question is: how does the system behave when collisions are both the source of dissipation and the source of energy. If the minimum dissipation principle was holding, we might expect to not have this $E$-$q_4$ correlation, but will it be different in our case?

----

### Self-assembly for the langevin model:


As we can see below, we observe almost no change when $\Delta t_N$ goes from $\Delta t_N<\tau_{coll}$ to $\Delta t_N>\tau_{coll}$. Nonetheless, I expect a stronger effective attraction as $\Delta t_N$ becomes larger and that the phase diagram will shift toward higher densities (In any case that would be an artefact of the simulation, not a physical effect).

### Self-assembly for the $\Delta$ model:


For the $\Delta$ model we observe a shift toward lower values of density. Since energy is injected through collisions, in this case, we might have an "effective repulsion" as we can see on the radial distribution function.

I had hope that in the $\Delta$ model the $\Delta$ could play a role in the self-assembly. But it looks like that it is unimportant... However, I don't think it simply renormalizes the time of the simulation. I didn't manage to prove it using the formula and it is -I think- easy to find a counter example. Moreover, simulations seem to indicate that $\Delta$ can play on the max of $g(r)$.



_______

#### Correlation between $E$ and $q_4$ for the delta model


We still observe correlations between $E$ and $q_4$ for the $\Delta$ model. This is somewhat surprising since collisions also add energy into the system and are not only dissipative.

Something I don't quite understand is that the angle between $\vec r_i-\vec r_j$ and $\vec v_i -\vec v_j$ increases with $q_4$ while I would expect the opposite (Indeed $<\frac{(\vec r_i- \vec r_j) \cdot (\vec v_i - \vec v_j)}{(|\vec r_i- \vec r_j||\vec v_i - \vec v_j|)}>$ goes closer to $0$ while the system crystallizes"). The s-s collisions behave strangely, I'll do a rerun with $N$ higher to see what's up. Moreover the mean dissipated energy doesn't increases significantly (probably because I'm simulating only 400 particles)


## Attractive EDMD
We model the attraction between particles by square wells

### Attraction + equilibrium

Two parameters:
the length of the well $\sigma_U$ and its depth $U$:

A funny thing is that we recover the weird bump observed at $r/r_s\simeq 2.1$ both for the langevin dissipative simulation and for the attractive equilibrium EDMD:

We observe a shift toward lower densities when we add attraction which is the opposite of what we expected (but still in agrreement with our intuition when $r_U$ is large). I'll do other runs with $r_U$ smaller.


Same thing for $U\sim T$ and $r_U = 1.06r$.

It seems that the attractive potential stabilises crystals. I'll check later if infinite pressure phases that weren't present at the equilibrium finite simulations could be find with attraction.
_____

As of now, we onbserve a shift toward high densities for $U\sim T$ and $r_U = 1.03r$:


We'll try to see if we observe the opposite effect with repulsive potential (should be fine?)

### Self-assembly with a matrix coefficient of restitution.

The coeff of restitution is now a matrix with: $\alpha_{ss}<\alpha_{sb}<\alpha_{bb}$. This is kind of realistic since $\alpha$ is an increasing function of the effective mass according to Pöschel and Schwager:

With $\alpha_{bb}=1$, $\alpha_{ss}=0.9$ and $\alpha_{sb}=0.95$ I didn't find any new phase.

There is maybe a propensity for straight line of particle (but I guess we can find it also at equilibrium...)

Also, I didn't try at equilibrium, but we observe a loooot, of T1+HexL

And a strange QC12 I guess:

__________________

At equilibrium we don't find back these phase segregation, but we find these dimers:


We don't find a phase separation for a constant coefficient of restitution.

### Attraction + dissipation
If we add dissipation between particles so that they get stuck in the well because of the loss of energy during the collision, we observe a "self-assembly" with square crystals (probably because -as proposed by Andrea- upfront collisions are more dissipative (and hence more effectively attractive) than the collisions on the sides):

| With dissipation | Without dissipation (and a very weak thermostat so that particles are stuck) |
| -------- | -------- |
| |  |

We observe this square lattice structure only when $r_U>1.41r$ for energy minimization reason I think (6 vs 8 neighbors)

Playing with the coefficient of restitution between big big, small small and small big we can obtain such things:


{%youtube PLE1hBVAMr0 %}

But most importantly, it looks extremely cool:

{%youtube MKSmikOa01o %}

{%youtube ocgiDtErOg4 %}

### Pressure

in the article **Global equation of state of two-dimensional hard sphere system** it is argued that the pressure can be defined out of equilibrium as:

$$\frac{pL^2}{NE}=1 + (1+ \alpha)\phi g_+(\phi)$$

From what I tried with my EDMD code, at $\alpha = 0.8$, the relation doesn't seem to hold. Indeed, by measuring the pressure as a change of momentum against a wall and comparing it with the pressure obtained from the formula above we find a disagreement for large $\phi$. Since I'm lazy, I dont really take the value of $g$ at contact but slightly farer (usually at $1.03r_b$ depending on my histogram of points). This underestimates $g_+$ but in any case, taking this into account would just worsen the disagreement between the 2 curves (since $p_{virial}$ would be greater).

When $dt_N$ is low enough the pressure is accurately measured.


In fact, as we can see in [Statistical mechanics of fluidized granular media: Short-range velocity correlations](https://www.researchgate.net/publication/12025443_Statistical_mechanics_of_fluidized_granular_media_Short-range_velocity_correlations) an exact formula is (as will be derived below by Frank):

$$p = \frac{NT}{V} - \frac{1}{2V}\left\langle\sum_{i< j}\vec r_{ij}\cdot\vec F_{ij}\right\rangle$$

Now, the force between particles $i$ and $j$ is $\vec F_{ij} = −\dfrac{1 + \alpha}{2}\dfrac{|\vec v_{ij} \cdot \vec r_{ij} | \vec r_{ij}}{|\vec r_{ij}|^2} δ(t − t^{coll}_{ij})$

Interpreting the average as a time average (the real physical meaning), we can the dirac delta, this gives:

$$p = \frac{NT}{V} + m\frac{1+\alpha}{4V}\left(\frac{1}{t}\sum_{coll}\vec r_{ij}\cdot\vec v_{ij}\right)$$.

The average value can also be integrated using the two particles distribution and the hypothesis that $f^{(2)}(1, 2) = g(\phi)f^{(1)}(1)f^{(1)}(2)$ (velocity are uncorrelated) to give the first equation we saw in the article without citation. The hypothesis breakdowns at high densities. In fact, for granular systems, the meaningful $g$ to take might be only the precollisional one.

### Wall force

Let's consider only the $x$ component of all velocities, and the quantity:
$$\left\langle \sum_i m v_x x \right\rangle,$$
where $i$ sums over all particles, and the subscript $i$ is implicit for both $v_x$ and $x$.
This should clearly be a constant in the steady state:
$$0=\frac{\partial}{\partial t}\left\langle \sum_i m v_x x \right\rangle = \left\langle \sum_i m \dot{v}_x x + \sum_i m {v}_x^2\right\rangle.$$
In the event-driven simulations, we can break this down exactly into the moments and intervals in time when the total value of this sum changes. In particular, the first term can be written as a sum over the points in time where velocities change, and vanishes during the time intervals between them. The second term has negligible contribution during any discrete event, but adds a contribution for each interval between two events. So:
$$0= \lim_{t\to \infty}\frac{m}{t}\left[\sum_{col} \delta x^{col} \delta{v}_x^{col} +
\sum_{wallcol} \delta x^{wall} \delta{v}_x^{wall} + \sum_{kick}x^{kick}\delta{v}_x^{kick} + \sum_{all} \sum_i v_x^2 \delta t_{last}
\right].$$
Here, we sum over all events: collisions, wall collisions, and Langevin kicks. The last sum sums over all events, with $\delta t_{last}$ indicating the time interval since the last event. One could keep a running total of $\sum v_x^2$ that gets updated during each event to avoid a full recalculation every time. Alternatively, one could obtain the average contribution of the last term by simply sampling the total kinetic energy during the simulations, but note that the kinetic energy typically will go up during kick events and down in between, so the sampling should not be in sync with those events.

Note also that the second term is simply minus the average force on a wall times the distance between the walls $L_x$. Here, $L_x$ is defined as the distance between the centers of two particles located at the left and right walls, respectively. So:
$$F_{wall} L_x = \lim_{t\to \infty}\frac{m}{t}\left[
\sum_{col} \delta x^{col} \delta{v}_x^{col} +
\sum_{kick}x^{kick} \frac{\gamma}{m}v_x + \sum_{all} \sum_i v_x^2\delta t_{last}
\right].$$
Here, we have used the fact that the random noise kicks will average out to zero (they are position-independent). Note that all instances of $v_x$ are the values *before* they are changed by the events.

The question is now how the wall influences the different terms. It will have some effect on the first and last terms, but these should become relatively small in a sufficiently large system. The change in total kinetic energy or the collision contribution should scale with the wall surface area (i.e. its length along $y$ in 2D), but the terms themselves scale with the full system size. In contrast, any wall effect in the second term would scale with the system size as well: the number of kicks near the wall scales with the wall surface area and the $x^{kick}$ factor scales with $L_x$. So, if we explicitly measure the average of each of these three terms in a simulation, do we find agreement with the above equation?

**ANSWER:** Yes, we find good agreement between the left and right term. At equilibirum and out of equilibrium (in this case we add a factor $(1 + \alpha)/2$ in front of the collision contribution) for a perfectly reflective wall, the kick term is negligible and on average contributes to nothing. $x$ and $v$ are indeed uncorrelated when we use an equilibrium wall (FALSE BECAUSE I'M NOT TAKING INTO ACCOUNT THE BINNING), nonetheless $x$ and $E$ are correlated when collisions are dissipative:


The orange and blue points are different simulations because I'm stupid:

(the spreading of the points is due to the lower number of collisions at low density for a fixed time)


(again the points spread mostly because of the decreases of energy with $\alpha$ and thus the decreases of the number of collisions/time).

_____

When the wall is dissipative, $x$ and $v$ are correlated:



The formula seems to work even without the kick term. Perhaps because $x^{kick}v$ is the time derivative of $x^2$ and thus 0 in the steady state

_______________

### Pressure and S1 self-assembly

Motivated by [this paper](https://www.nature.com/articles/srep28726.pdf). We'll try to explore the pressure change during the self-assembly.

_______

pressure at equilibrium as a function of $\phi$ for $N=15000$ particles:

Still at equilibrium, we might observe a maxwell loop:

### Attractive Equilibrium self assembly

A lot of dimers:


Strange formations:

QC8 are easy to form:


### [16/11/2022]

When we start from the phase separated dimers formed with $\alpha_{bb}=1$, $\alpha_{ss}=0.9$ and $\alpha_{sb}=0.95$ it melts at equilibrium. This is especially signficant since usually, structure formed in a dissipative simuilation are stable at equilibrium (S1, QC12 or QC8). Most of the time, out of equilibrium, the crystals form in a region of phase space where they are dynamically arrested at equilibrium.


The same thing goes for out of equilibrium with a constant $\alpha$:


We still observe dimers for both of them, but it's not particularly clear.

______________

I've written a script to automatize the direct coexistence method:
{%youtube QawB-bGuUU4 %}

I can almost fix my density (at least when $L\rightarrow\infty$ since i've boundary issues).

### [23/11/2022]
With **LAMMPS** we also observe a decrease of $E$ with $q_4$ when we melt a crystal.


### [25/11/2022]

#### Trying to get the effective coefficient of restitution
In order to see what model is the correct one ($\Delta$ model or exponential one).

In order to do this, we make particles collide:


Blue one are still, red are moving with different velocities. We then calculate the coefficient of restitution with $v_x$:


The first issue is that $\alpha$ will depend on $XY$ collision angle, this approach doesn't take that into account. A second and more severe issue is that the $\alpha$ will also depend on the $ZX$ angle (the height at which the particle will collide).

For example, when particles are set in a configuration that lead to a synchronous effect (particles are colored by $z$ below). If we wait long enough, all particles will synchronize:
{%youtube ejNQ3FQLabs %}

and the coefficient of restitution will be trivial (since there is no transfert of the $Z$ momentum to the $XY$ plane):


Now, if we don't let the system "synchronize" (but still let it jump a bit so that it gains a $z$ velocity coherent with the amplitude used):


______

For the system that is not synchronous in the SS:


I suspect that, with this configuration, they are still too much synchronized (I used a slab of height 1.1*diameter, that play an important role also) since we can't observe a coefficient of restitution greater than 1 even for low velocities:


### [28/11/22]

We observe a loop but it is hard to get a proper value for $p_{coex}$:


### [29/11/22] (slamming my head against lammps)

I tried to reproduce Matteo's result with his script. But have a hard time doing so:

#### Synchronisation

Here are his results for synchronisation. I assume the shaker frequency to be $53$Hz.

| $r=2.5mm$ | $r= 1.25m$ |
| -------- | -------- |
|  | |

From what I found (by running a lot of simulations **with his script** and trying to see where I obtained the same results), for $r=1.25mm$, the region in blue below is not asynchro but syncho (as for $r=2.5mm$), the other part seems fine.


Also the top-left triangle of synchro takes way more time to synchronize than the other part of the phase diagram.

I'm starting with particles at 0 XY velocity and a random $z$ position.
___________

Now, For the same initial configuration as well as the same amp, Hz and height. It seems that my and his script give the same physics.

________

I have issue reproducing the same curve for the dead/living transition, but I think this is caused by different $r$.

### [30/11/22]

Pressure in a binary mixture:

Let's start from the virial theorem (though we could start from an stochastic equation as Frank did, but it would be a bit more painful, also, let's forget about the reservoir contribution). In a NESS, we have a fixed temperature:

$$p = \dfrac{NT}{L^2}-\frac{1}{2L^2}\left\langle\sum_{i< j}\vec{r}_{ij}\cdot\vec{F}_{ij} \right\rangle$$

$F_{ij}=m_i\left .\frac{dv_i}{dt}\right|_{coll-ij}$, using the collisions rules and the fact that a collision is instantaneous:
$$F_{ij}=m_im_j\dfrac{1+\alpha}{m_i + m_j}\dfrac{|\vec v_{ij} \cdot \vec r_{ij} | \vec r_{ij}}{|\vec r_{ij}|^2} δ(t − t^{coll}_{ij})$$

So finally:

$$p = \dfrac{NT}{L^2}-\frac{1+\alpha}{2tL^2}\sum_{coll-ij}\dfrac{m_im_j}{m_i + m_j}\vec{r}_{ij}\cdot\vec{v}_{ij}$$

When compared with the wall pressure for a binary mixture, we find perfect agreement between the two pressures.

Let's analyze the scaling of $P$.

#### m varied for a monodisperse system
First consider the case where we have only one type of particle:

$$pV = NT - \frac{(1+\alpha)m}{4t}\sum_{coll}\vec r_{ij}\cdot\vec v_{ij}$$

For hard spheres, $pV\beta/N$ depends only on the structure factor, so the second term of:

$$pV\beta/N = 1 - \frac{(1+\alpha)m}{4NTt}\sum_{coll}\vec r_{ij}\cdot\vec v_{ij}$$

should also be invariant by a change of $m$. At $T$ fixed, if we change $m$, $v$ will change by an order of $\sqrt{T/m}$, so by a factor $1/\sqrt{m}$, thus the dot product scales as $1/\sqrt{m}$ ([exactly: ](https://arxiv.org/abs/1211.1645) $\langle b_{ij}\rangle=-2\sigma\sqrt{{T\pi/m}}$ **at equilibrium**). This is not enough to cancel the factor $m$ in front of the sum, we also have to take into account that if the system slows down due to the decrease of velocity, the time will also slow down by the same amount, so the number of collisions/timescale will scale as $1/\sqrt{m}$. Thus we recover that $p$ at equilibrium is independent on $m$.

#### m varied for a polydisperse system

At equilibrium, since we have equipartition between particles. A polydisperse system will also behave the same way and $p$ will be invariant byu a change of $m$ at $T$ fixed. But this is not so trivial out of equilibrium where each type of particles will have a different temperature. For two types of particles $s$ and $b$:

$$pL^2 = N_bT_b + N_sT_s-\frac{1+\alpha}{2t}\left( \frac{m_b}{2}\sum_{b-b}\vec{r}_{ij}\cdot\vec{v}_{ij}+\frac{m_s}{2}\sum_{s-s}\vec{r}_{ij}\cdot\vec{v}_{ij}+\dfrac{m_bm_s}{m_b + m_s}\sum_{b-s}\vec{r}_{ij}\cdot\vec{v}_{ij}\right)$$

Each sum will scale differently with particle mass and temperature...

### [01/12/22]

Here are my results with Matteo's script. This is somewhat consistent with his results, knowing that the initial configuration might play a role on the degree of synchronisation (if we start close to the ceiling I expect particles to syncrhonise more since they might be able to bounce on it). But still, there is a difference between his runs and mine (in the LAMMPS code) I think.

Also, for an amp>0.04 we see that some system don't synchronize but still, particles touch the ceiling.


Here is with my code:


They are essentialy the same

### [13/12/22]

#### Average over collisions for the boltzmann equation
With Andrea, we were thinking about getting the curve of dead/active particles using the dissipation function for the exponential model.

As a reminder, $G$ is the rate of dissipation per unit of time (and probability not unit of area as said in the article [Hydrodynamic modes in a confined granular fluid](https://www.dropbox.com/preview/GRANULAR/ABSORBING/TO_READ/Hydrodynamic%20modes%20in%20a%20confined%20granular%20fluid.pdf?context=standalone_preview&role=personal)). If the system is homogeneous (opposite convention than what is used in the article):

$$\dfrac{\partial T}{\partial t} = G$$

Thus, in the steady state, $G(\phi, T)$ must be 0, allowing one to get $T(\phi)$, which gives us for the $\Delta$ model:


The shape of the curve is good but the temperature is not (I think). In any case, the ideal gas mean free path has been used in the formula which is not appropriate for $\phi>0.2$.

In the $\Delta$ model the energy is injected and dissipated at collision (which happens at a frequency $\omega/2$) and dissipated by continous friction so:

$$G = \frac{\omega}{2}\langle E' - E\rangle_{collisions} - \left.\dfrac{\partial T}{\partial t}\right|_{free flight}$$

Starting from the collisions rule we have:

$$\langle E' - E\rangle_{coll} = m\Delta^2 + m\langle|\vec v_{12}\cdot \hat r_{12}|\rangle_{coll}\alpha \Delta - m\langle|\vec v_{12}\cdot \hat r_{12}|^2\rangle_{coll}\dfrac{1 - \alpha^2}{4}$$

with $\langle|\vec v_{12}\cdot \hat r_{12}|\rangle_{coll} = \sqrt{\pi T/m}$ and $\langle|\vec v_{12}\cdot \hat r_{12}|\rangle^2_{coll}=4T/m$.

and

$$\left.\dfrac{\partial T}{\partial t}\right|_{freeflight}=0.5m\left(2\vec v \cdot \frac{d\vec v}{dt}\right) =-2(0.5m\vec v^2)=-2T $$

Which gives us:

$$G = \frac{\omega(\phi, T)}{2}(m\Delta^2+\alpha\Delta\sqrt{\pi m T}-T(1-\alpha^2))-2T$$

____

Now, let's find how to calculate over average.

$$\vec v_{12}\cdot \hat r_{12}=|\vec v_{12}|\cos(\theta)$$

where $\theta$ is the angle between the velocity difference and the one joining the center of the two particles colliding, normalized to unity.

since $\vec v_{12} = \vec v_1 - \vec v_2$ both taken in $\mathcal N(0, T/m)$ we have $\vec v_{12}\sim \mathcal N(0, 2T/m)$:

Now, as I'm lazy and *slightly* confused about cross sections, I plotted the histogram of $\theta$ instead of calculating it :) :


We see that $P(\theta)=-\cos(\theta)/2$ if $\pi/2<|\theta|<\pi$ (proof by: **we clearly see it**). Otherwise, no collision happen.

BLABLA to change
Now, if $Y = \cos(\theta)$, the distribution of $Y$ is:

$$P_Y(y) = \sum_{y = -\cos(\theta)/2}P(\arccos(y))\left|\dfrac{d\theta}{dy}\right| = \sum_{y = -cos(\theta)/2}\frac{-\cos(\arccos(y))}{2}\left|\dfrac{2}{\sin(\arccos(y))}\right| = -\dfrac{y}{\sqrt{1 - y² }}$$

Which agree with what I find numerically. It was just for fun since I don't need it..

In this case:

$$\langle|\vec v_{12}||\cos(\theta)|\rangle = \int_{0}^{\infty}dv\dfrac{mv^2}{2T}e^{\frac{1}{2}\frac{mv^2}{2T}}\left(2\int_{\pi/2}^{\pi}-\dfrac{\cos^2(\theta)}{2}d\theta\right) = \sqrt{\frac{\pi}{2}\frac{2T}{m}}\times\frac{\pi}{4}$$

Which is not what we were looking for. I an article I've found [a mean value for which we include a factor $\langle|\vec v_{12}||\cos(\theta)|\rangle$ directly (equation (11))](http://lptms.u-psud.fr/membres/trizac/Articles/Pre_Matthieu02.pdf) , probably to ensure that the velocity sampled at this moment is juste before the collision, but still, with this it doesn't look like I can recover the correct value.. That would give us a $cos^3(\theta)$ and a $v^3$ which don't give a $\pi$).
I've found the issue, see below.

#### Coexistence:


I ran my script for the coexistence method with 40k particles at equilibrium:


Left is the pressure according to density, mid is the pressure ''along'' the $x$ axis (by only taking into account the $x$ momentum transfer). Right is $q_4$.

This is pretty horrible, I think the script hasn't run long enough (which is surprising since I let it run during a lot of time...). For example, we still observe a crystal at the lowest packing fraction which was not the case when we started from a liquid configuration, this is clear here:

{%youtube 0fW7WK4LGKc%}

Also clear that the crystallization is not finished for high packing fraction:

{%youtube WVx3y_Y3ap0%}

#### First order transition with dissipation

I tried to see the first order transition when the system is dissipative starting from an initial configuration fully liquid, indeed, we observe the maxwell loop. When we normalize by the energy of the system we obtain a pressure close to (but higher) than the coexistence pressure of Etienne ($\alpha = 0.96$ $T = 5$ $dt_N = 0.2$).


### [14/12/22]

#### Correct calculations of the mean values for the dead state

Indeed, we have to take into account the fact not all particles can participate in a collision. The correct formula is:

$$\langle A\rangle_{coll} =\dfrac{ \int_{0}^{\infty}\mathrm{d}v\int_0^{2\pi}\mathrm{d}\theta\Theta(-v\cos\theta)|v\cos\theta| A \frac{mv}{2T}e^{-\frac{1}{2}\frac{mv^2}{2T}}}{ \int_{0}^{\infty}\mathrm{d}v\int_0^{2\pi}\mathrm{d}\theta\Theta(-v\cos\theta)|v\cos\theta| \frac{mv}{2T}e^{-\frac{1}{2}\frac{mv^2}{2T}}}\tag {1}$$

The $\Theta$ term allows only configurations that could lead to a collision and the $v\cos\theta$ term gives (somehow, I'm not so clear with this, it's explained in the book: granular gases) the good statistic of $v$ (the magnitude of the difference of the velocities). We supposed that the process is isotropic (which is not the case when there is dissipation). We also assumed a gaussian distribution function which breakdown at the transition and before the transition.

So finally, we have:

$$\langle A\rangle_{coll} =\dfrac{\displaystyle{\int_{0}^{\infty}\mathrm{d}v \frac{mv^2}{2T}e^{-\frac{1}{2}\frac{mv^2}{2T}}}\int_{\pi/2}^{3\pi/2}\mathrm{d}\theta|\cos\theta| A}{2\sqrt{\dfrac{\pi T}{m}}}\tag{2}$$

With $A=v\cos(\theta)$ or $A=(v\cos(\theta))^2$ we find the results of the article.

#### Exponential model

Let's get back to:

$$G = \frac{\omega}{2}\langle E' - E\rangle_{collisions} - \left.\dfrac{\partial T}{\partial t}\right|_{free flight}\tag{3}$$

Starting from the collision rule we have:

$$\langle E' - E\rangle_{coll} = - m\left\langle|\vec v_{12}\cdot \hat r_{12}|^2\dfrac{1 - \alpha_0^2e^{-2v/v_0}}{4}\right\rangle_{coll} = - \dfrac{m}{4}\left(\langle|\vec v_{12}\cdot \hat r_{12}|^2\rangle_{coll}- \alpha_0^2\left\langle|\vec v_{12}\cdot \hat r_{12}|^2e^{-2|\vec v_{12}|/v_0}\right\rangle_{coll}\right)\tag{4}$$

The first mean value is already known:

$$ - \dfrac{m}{4}\langle|\vec v_{12}\cdot \hat r_{12}|^2\rangle_{coll}=-T\tag{5}$$

The second term has to be calculated (:( I'm tired of integralsssss)):

$$\left\langle|\vec v_{12}\cdot \hat r_{12}|^2e^{-2|\vec v_{12}|/v_0}\right\rangle_{coll} = \dfrac{\displaystyle{\int_{0}^{\infty}\mathrm{d}v \frac{mv^4}{2T}e^{-\frac{1}{2}\frac{mv^2}{2T} - 2v/v_0} }\overbrace{\int_{\pi/2}^{3\pi/2}\mathrm{d}\theta|\cos(\theta)|\cos^2\theta}^{4/3}}{2\sqrt{\dfrac{\pi T}{m}}}\tag{6}$$

The gaussian integral looks horrible with mathematica...

### [18/12/22]

$$\displaystyle{\int_{0}^{\infty} av^4e^{-\frac{1}{2}av^2 - bv}}\mathrm{d}v=\dfrac{\sqrt{\pi/2}(3a^2+6ab^2+b^4)e^{b^2/(2a)}\mathrm{erfc}(b/\sqrt{2a})}{a^{7/2}}-\dfrac{5ab + b^3}{a^3}\tag{7}$$

with $a=m/2T$ and $b=2/v_0$.

Here is the plot of the mean dissipated energy for one collision when the system is at temperature $T$ for a set of $a_0$. We correctly predict that for low temperatures, the system will gain energy by collisions and will lose energy when the temperature is high. The point at which $\Delta E = 0$ is the steady state temperature (for now, we haven't yet added dissipation while moving (the $\gamma v$ factor)).


### [20/12/22]

Without drag the steady state temperature is independent on $\phi$ since the only source of energy of the system is collision dependent. In this case (**without drag**), as in the $\Delta$ model:
$$G=F(\phi, T)\Delta E(T)$$

Where $F$ is the frequency of collision and $\Delta E$ is the energy change for one collision. The steady state is given by a balance of the energy source(s) so $G = 0$.

$F = 0$ is only realized at $T = 0$ (no collisions means zero temperature at finite density). So if we want a non trivial steady state temperature we have to look for the zeros of $\Delta E$ which can't depend on $\phi$. So for model which only add/remove energy into the system by collision, the steady state temperature is density independent.

We could calculate $T_{ss}$ and check it with the value obtain in simulation, but I don't know the usefulness of such a thing.

Also, this is great because it indicates that the energy of a binary system can be maintained at a constant level regardless of the density.

End of the useless rant...
___

Now, let's find the steady state temperature when drag is added.

$$G =\dfrac{\omega}{2}\left[\alpha_0²m/6\sqrt{\dfrac{m}{\pi T}}\left(\dfrac{\sqrt{\pi/2}(3/9(m/T)^2+3(m/T)b^2+b^4)e^{Tb^2/m}\mathrm{erfc}(b/\sqrt{m/T})}{\left(\dfrac{m}{2T}\right)^{7/2}}-\dfrac{5/2(m/T)b + b^3}{\left(\dfrac{m}{2T}\right)^3}\right)-T\right]-2\gamma T \tag{666}$$

The situation is not as simple as for the $\Delta$ model. For the $\Delta$ model the dissipation $G$, I think, had only at most one zero (without counting the trivial $T = 0$).

For the exponential model:

- At low $\phi$, $G$ is monotonous and never touches $0$ except at $T = 0$ (blue).
- At the critical $\phi = \phi_c$, $G$ touches one time the $0$ (at $T\neq 0$) (orange).
- At high $\phi$, $G = 0$ has two non trivial solutions one stable and one instable (green).

Left is $G$, right is $U$ with $G = -\partial_x U$

To know which of the two points is the stable one, recall that $G$ is the dissipation/energy injection rate. If $G>0$ the energy injection rate is $>0$ and the opposite if $G<0$. Now, take, the rightmost zero of the green curve (let's call this temperature $T_c$). If we disturb the temperature by a tiny negative amount (at $T_c-\delta T$), we will place ourselves in a region where $G(T_c-\delta T)>G(T_c)$ which means that the temperature of our system will effectively rises and we will reach again $T_c$. The same thing goes for a tiny positive perturbation (this is just linear perturbation theory of a fixed point...). For the left fixed point, we have the opposite effect. Which means that the left fixed point is unstable while the right one is stable.

I'll have to check the stability of the fixed point of the $\Delta$ model but I think this is fine...

Here is the phase plot

I didn't check, but I think the temperature and the critical packing fraction are very wrong. For he packing fraction, this is ok since we didn't used the real collision frequency but the ideal gas one. For the temperature this is a bit more troublesome.

_____

Here is a comparison of the theoretical phase diagram:

And the Numerical one:


____
Also, one can find the critical surface:


The same things can be done for the $\Delta$ model, I'll do them later.

### [02/01/23]

#### Coexistence

Simulations of coexistence are still running but are taking some times since they are long in order for them to reach the correct steady state (and with a lot of particles to reduce finite size effect).

This time, the initial configuration is this one:


A liquid configuration is created, as before, and a perfect S1 crystal is generated within a square area. The crystal is then deformed by elongating it in the y direction by streching it of $2dL$. In prior simulations, a gap of size $2dL$ was left at the top and bottom of the crystal, leading to the formation of liquid at the top and bottom of the box upon initiation of the simulation.

#### absorbing state

Here is the phase diagram for a varying $v$:


The larger v is, the closer we are to the experimental one:


______
For $\gamma$ this doesn't work as well, when we have too muchd dissipation, the theorical phase diagram is plain wrong (IO forgot the axis, $y$ is $\gamma$ and $x$ is $\phi$):


experimental: 

We can expect it to be a little better with real collision frequency.

#### Virial pressure for collision based energy injection model

We have:
$$p = \dfrac{NT}{L^2}-\frac{1}{2L^2}\left\langle\sum_{i< j}\vec{r}_{ij}\cdot\vec{F}_{ij} \right\rangle$$

$F_{ij}=m_i\left .\frac{dv_i}{dt}\right|_{coll-ij}$, using the collisions rules and the fact that a collision is instantaneous for the exponential model we have:
$$F_{ij}=m_im_j\dfrac{1+\alpha(v)}{m_i + m_j}\dfrac{|\vec v_{ij} \cdot \vec r_{ij} | \vec r_{ij}}{|\vec r_{ij}|^2} δ(t − t^{coll}_{ij})$$

and get back the same result as before.

Now for the $\Delta$ model, the collision rule is different and so will be the force. To the usual collision rule for the $i$ particle, we add a term:

$$-2\Delta\dfrac{m_j}{m_i + m_j}\dfrac{\vec r_{ij}}{|\vec r_{ij}|}$$

So we obtain:

$$p = \dfrac{NT}{L^2}-\frac{1+\alpha}{2tL^2}\sum_{coll-ij}\dfrac{m_im_j}{m_i + m_j}\vec{r}_{ij}\cdot\vec{v}_{ij}+\dfrac{\Delta}{tL^2}\sum_{coll-ij}\dfrac{m_im_j}{m_i + m_j}|r_{ij}|$$

the term added by the $\Delta$ is independent on the angle of collision, as expected.
### [03/01/23]

#### $\Delta$ model

Here is the phase diagram calculated from kinetic theory for the $\Delta$ model:


Versus the "experimental" one:


It fails at high density which is expected since the collision frequency used for the computation is the one computed starting with the ideal gas.

____

The behavior at the "transition" is also different:

The same parameters are used but the $x$ axis is different for both plot.
Theory:

Experiment:

With the mean free path theory, the "critical point" was $\alpha$ independent, this seems to be more in line with what we obtain by simulations.

### [06/01/23]

#### Collision frequency

Since I'm a massive idiot, here is the new freq vs phi, where there is almost no difference between eq and dissipative systems:


Note the change of the cristallization location:


_____

**WRONG**
Here is a plot of the excess kurtosis for simulations without $\gamma$:

$$a_2 = \dfrac{1}{2}\dfrac{\langle v^4\rangle}{\langle v^2\rangle^2} - 1$$

It seems, alright to use a gaussian approximation (not sure about the sign of the excess kurtosis though... It's strange) except maybe for the exponential model.


Of course, with $\gamma$ the distribution will go even farer from a gaussian.

_____

update to the wrong thing above:

Expo model:


Delta model:

______

Now that we are good with that, let's see if we can refine the results, taking into account the new $\omega$:

### Delta model

| Theoretical | Simulation |
| -------- | -------- |
| | |

good agreement, we even see the change of shape of the transition curve below $\Delta = 0.275$


### Expo model

| Theoretical | Simulation |
| -------- | -------- |
| | |

great agreement (the weird patterns are caused by my root finder which is terrible). Way better than what we had with the ideal $\omega_{ideal}$:


________________

| Theoretical | Simulation |
| -------- | -------- |
| ||

________

| Theoretical | Simulation |
| -------- | -------- |
|||


The high kurtosis prevents us from having a nice result here.

### [10/01/23]

#### Can we make real phase transition with the $\Delta$ model?

By real, I mean $T = 0$ for a finite $\phi$. I had in mind that changing the friction $dv/dt = -\gamma \mathrm{sign(v)} \left|v\right|^{2\nu - 1}\tag{10}$:

$$G = \frac{\omega(\phi, T)}{2}(m\Delta^2+\alpha\Delta\sqrt{\pi m T}-T(1-\alpha^2))-2^{\nu}\gamma T^{\color{red}\nu}\tag{11}$$

could lead to a real transition:

Indeed with $m = 1$ and $\omega = 2\Omega \sqrt{T}$:
$$\Omega(\phi)\sqrt{T}(\Delta^2+\alpha\Delta\sqrt{\pi T}-T(1-\alpha^2))-2^{\nu}\gamma T^{\color{red}\nu} = 0\tag{12}$$
implies that
$$\Delta^2+\alpha\Delta\sqrt{\pi T}-T(1-\alpha^2)-2^{\nu}\gamma'(\phi) T^{\color{red}\nu - 0.5} = 0 \tag{13}$$

For $\nu > 0.5$, $T = \infty$, for stability, the RHS must goes to $-\infty$ and at $T \rightarrow 0$, the RHS has a positive value for a finite $\phi$. There must be a 0 inbetween (different than 0. Actually there are sometimes 2 **stable** zeros:


). So no real phase transition can happen in the $\Delta$ model for this values of $\nu$ since $T_{ss}$ is always different than 0 (except at $\phi = 0$).

For $\nu = 0.5$ we have either a first order transition or no phase transition. It depends on the value of $\Delta^2 - 2^{\nu} \gamma'(\phi)$ (and the other parameters). Since a second like order transition would require an arbitrary small $\phi$ solving the equation, it would make $\gamma'$ explodes, we can't have it, here.

Here, for example, we have a first order phase transition for $\nu = 0.5$ with the $\Delta$ model:


Same thing for $\nu < 0.5$

$\nu = 0.5$ corresponds to $\dfrac{dv}{dt} = -\gamma$, this is dry friction.

We can try to see if this is what we observe in lammps. Namely, if when the transition is discontinuous, we have dry friction, and otherwise, the transition is continuous.

In order to obtain a clean graph, we need to start from a small velocity, in the case $v\sim 0.5$ (what units???? lammps says lennard jonnes but it doesn't make any sense, perhaps 0.05m/s" (I hate this thing so much :) ). We can place ourserlves at any parameter. Here is the plot of certain realizations for amp = 0.08 and different $h$. $h = 1.55$ corresponds to an asynchronous situation while the other ones are synchronous. The $x$ axis is the timestep and the y is the velocity. We start from a z randomized.


Here is the mean value of a lot of realizations:

Here is a theoretical dry friction velocity in logscale:


It seems that at high velocity, we always have some kind of dry friction and at low velocity, viscuous one.


This doesn't look promising
_____

For the exponential model, analytics are a bit harder... Numerical computations seems to say that we can only obtain first order phase transition.

#### Trizac's model

Trizac model considers a random coefficient of restitution.

In his article, he requires $\langle\alpha^2\rangle = 1$ so that $\langle E' - E\rangle = 0$ in order to have a steady state. We can give a little bias such that:
$$\langle E' - E\rangle > 0$$

Now, we're not sure that the system converges to a finite temperature. This is the main point of collision based energy injection model, we can add a velocity dependent cutoff that imposes $G(T\rightarrow \infty)= -\infty$ **whatever** the drag model used is that lead to a well define steady state for all $\phi$.

For trizac model, we obtain (with $\omega(T, \phi)/2 = \Omega(\phi) T^{0.5}$):
$$G = \Omega(\phi)T^{3/2}(\langle\alpha^2\rangle - 1) -2^{\nu}\gamma T^{\nu } \tag{14}$$

$T = 0$ is always a solution, the other 0 are given by:

$$\Omega(\phi)(\langle\alpha^2\rangle - 1) -2^{\nu}\gamma T^{\nu - 3/2} = 0 \tag{15}$$

If $\nu<3/2$ we see $G(T = \infty) = \infty$, the system explodes, also, we can prove that in this case, the zeros are unstable. It means that we are looking for a drag that reads: $dv/dt = -\gamma \mathrm{sign(v)} \left|v\right|^{\beta}$ with $\beta>2$. Not very realistic.

In any case, this model only gives second order like phase transition, that are not very interesting.

### [13/01/23]

#### The synchro $\Delta$ model with linear drag:

**IN ALL THE PLOTS BELOW $\tau_s\gamma$ is wrong. Every value has to be multiplied by $10^{-4}$**

Andrea proposed a way to incoporate the degree of synchroneity in the $\Delta$ model.

We suppose that after a time $\tau_s$ without collisions, a particle start to be synchronized. Two particles synchronized will have a collision with a $\Delta = \Delta_s$ while collisions with at least one async particle will have a $\Delta = \Delta_a>\Delta_s$. For example, if we put $\Delta_s = 0$, we enforce the fact that synchroinized particles colliding won't transfer Z momentum in the XY plane and that the collision will only be dissipative.



_____
#### Behaviour as $\tau_s$ goes to 0.

We go from a first order to a second order transition at finite $\tau_s$


#### Behaviour with non zero $\Delta_s$

Here we have a non zero $\Delta_s$,

Below, $\tau_s\gamma = 50$ $\Delta_s/\beta = 2$ and $\Delta_a = 3$

For these cases, the green curve interpolates the two simple $\Delta$ model curves
Below, $\Delta_s/\beta = 1$ and $\Delta_a/\beta = 3$


With a $\Delta_s/\beta = 0.1$, $\tau_s\gamma = 200$


Zoomed

Here the transition is different when $\Delta$ goes from 0 to some very small finite value


While here it is not ($\tau_s\gamma = 50$) for the same values of the deltas:

### [16/01/20]
We see that, as expected, when the inverse of frequency of collision goes to $\tau_s$ we recover the upper curve.

See also:


_______________________________________

This is expected but I was expecting to be able to simply add a heavyside function into the root finder in order to simulate the switch from one delta to another but it does not seem possible. For example, if we place ourselves at $\phi = 0.45$ above and color particles according to their synchronicity, we see that it's not just: "every particle is either synchro or asynchro".


_______

I also tried to obtain a(n unphysical) jump like this:
because I was afraid that the model could give something similar, but some tests seem to indicate that for reasonable values of parameters, this kind of thing don't appear.

________
Also the Delta model REALLY goes to 0 for a finite value of $\phi$ (dots) while the theoretical prediction doesn't (plain curve).

### [17/01/23]

#### Theory of the synchro $\Delta$ model

With the synchro $\Delta$ model, we can keep the usual expression for $G$ but use a $v$ dependent $\Delta$:

$$\Delta(v) = \Delta_s\Theta (v) + (\Delta_a - \Delta_s)\Theta(v - l/\tau_s)\tag{20}$$

With $l$ the mean free path. This type of expression foreshadows a significant amount of pain associated with the erfc function. :) Using a continuous function like a sigmoid or a Fermi-Dirac function could help alleviate the pain caused by this horror. Anyway...

Recall that: $$\langle A\rangle_{coll} =\dfrac{\displaystyle{\int_{0}^{\infty}\mathrm{d}v \frac{mv^2}{2T}e^{-\frac{1}{2}\frac{mv^2}{2T}}}\int_{\pi/2}^{3\pi/2}\mathrm{d}\theta|\cos\theta| A}{2\sqrt{\dfrac{\pi T}{m}}}\tag{21}$$

And ($m = 1$):$$ E' - E = \Delta^2 + |\vec v_{12}\cdot \hat r_{12}|\alpha \Delta - |\vec v_{12}\cdot \hat r_{12}|^2\dfrac{1 - \alpha^2}{4}\tag{22}$$

Let's go term by term:

##### First term:
We suppose that $\Theta(a) = \Theta(a)^2$, which is definitively not right at $a = 0$ but I decided that I didn't care and that it should be ok nonetheless. Taking the square of equation (20) and using $\Theta(v - l/\tau)\Theta(v) = \Theta(v - l/\tau)$:
$$\Delta^2 = \Delta_s^2\Theta(v) + (\Delta_a^2 - \Delta_s^2)\Theta(v-l/\tau)\tag{23}$$

So averaging gives:

$$\langle \Delta^2\rangle = \Delta_s^2 +(\Delta_a^2 - \Delta_s^2)\dfrac{\displaystyle{\int_{l/\tau}^{\infty}\mathrm{d}v \frac{v^2}{2T}e^{-\frac{1}{2}\frac{v^2}{2T}}}}{\sqrt{\pi T}} = \Delta_s^2+(\Delta_a^2 - \Delta_s^2)\left(\frac{l/\tau}{\sqrt{\pi T}} e^{-\dfrac{(l/\tau)^2}{4T}}+ \mathrm{erfc}\left[\dfrac{l/\tau}{2\sqrt{T}}\right]\right)\tag{24}$$

Sanity check:
- if $\tau = 0$, the system synchronizes immediately, and indeed: $\langle\Delta^2\rangle = \Delta_s^2$.
- if $\tau = \infty$ the system never synchronizes, and indeed: $\langle\Delta^2\rangle = \Delta_a²$.

##### Second term:

$$\langle|\vec v_{12}\cdot \hat r_{12}|\alpha \Delta \rangle = \alpha\Delta_s\sqrt{\pi T}\tag{25} + \alpha(\Delta_a - \Delta_s)\overbrace{\dfrac{1}{4}\sqrt{\dfrac{\pi}{T}}}^{\mathrm{cos } + norma}\displaystyle{\int_{l/\tau}^{\infty}\mathrm{d}v \frac{v^3}{2T}e^{-\frac{1}{2}\frac{v^2}{2T}}}$$
Which gives:
$$\langle|\vec v_{12}\cdot \hat r_{12}|\alpha \Delta \rangle = \alpha\Delta_s\sqrt{\pi T}\tag{26} + \alpha(\Delta_a - \Delta_s)\dfrac{1}{4}\sqrt{\dfrac{\pi}{T}}\left(4T + (l/\tau)^2\right)e^{-\dfrac{(l/\tau)^2}{4T}}$$

Sanity check:
- if $\tau = 0$, the system synchronizes immediately, and indeed: $\langle|\vec v_{12}\cdot \hat r_{12}|\alpha \Delta \rangle = \alpha\Delta_s\sqrt{\pi T}$.
- if $\tau = \infty$ the system never synchronizes,and indeed: $\langle|\vec v_{12}\cdot \hat r_{12}|\alpha \Delta \rangle = \alpha\Delta_a\sqrt{\pi T}\tag{27}$.

##### Third term:

No dependance on $\Delta$ (thank god...), it's the usual collision term:
$$\left\langle |\vec v_{12}\cdot \hat r_{12}|^2\dfrac{1 - \alpha^2}{4} \right\rangle = T(1-\alpha^2)\tag{28}$$

----------

#### Gluing the terms:

$G = \dfrac{\omega}{2}\left(\Delta_s^2+(\Delta_a^2 - \Delta_s^2)\left(\frac{l/\tau}{\sqrt{\pi T}} e^{-\dfrac{(l/\tau)^2}{4T}}+ \mathrm{erfc}\left[\dfrac{l/\tau}{2\sqrt{T}}\right]\right) + \alpha\Delta_s\sqrt{\pi T} + \alpha(\Delta_a - \Delta_s)\dfrac{1}{4}\sqrt{\dfrac{\pi}{T}}\left(4T + (l/\tau)^2\right)e^{-\dfrac{(l/\tau)^2}{4T}} - T(1-\alpha^2)\right) - 2\gamma T\tag{29}$

We manage to make the $\Delta$ model as horrible as the exponential model!

________
The results are "meh":
For $\Delta_s =0$ we always get a first order phase transition


Which seems to be false here:

______

Sometimes, we can't get a second order phase transitio:



This is especially true when the transition happens at a packing fraction for which the classical delta model with $\Delta = \Delta_s$ **is at a low temperature.**

______

If the transition happens at a high packing fraction (with a non zero $E$) it can give a second order phase transition:


For the purple curve, this is really a first oder PT. For the violet one, a bump is created, as with the first order phase transitions, but it's created above 0. So really, the PT is a 2nd order one.

| purple curve | green curve |
| -------- | -------- |
| ||

### [18/01/23]

#### Focusing on the fat tail

As noted above for the synchro $\Delta$ model, at the transition it seems that we have a very fat tail and that this might be the very source of disagreement between theory and simulations (maybe we have some kind of phase transition with the kurtosis as an order parameter). Let's see this in more details.

We start from this

For example, here, we use very high values for $\Delta_a$ and $\tau_s$ so that the transition can happens at a "low" packing fraction. In this case the agreement is better (probably because the transition happens at a point for which the simple $\Delta$ model with $\Delta = \Delta_s$ would give a low energy):


But the tail is very fat, excess kurtosis = 0.36 (for the plot above, at $\phi = 0.35$ the excess kurtosis was equal to 0.99):

____________
Pushing the values farther, we find again a disagreement:

And a huge tail (excess kurtosis = 1.35!!!!):

In linear axis:

_____
##### RECAP
From what I understand:
- the fatter the tail at the transition, the greater is the deviation from the theory (it is pretty resilient. Excess kurtosis < 0.5 seems to still gives good results).
- A gaussian model seems to just interpolate the two curves of a simple delta model.
- To have a good agreement between theory and simulation, it seems that if the value of $E$ at $\phi_c$ corresponding to the delta model when $\Delta = \Delta_s$ or $\Delta_a$ should be close together. By that I mean that the transition of the **synhcro** model has to happen at a $\phi_c$ such that the energy of the **simple** delta model satisfies:
$$E(\phi_c, \Delta = \Delta_s) \simeq E(\phi_c, \Delta = \Delta_a)$$
This seems necessary but not sufficient.

_______

Let's try to reproduce this fat tail with a simpler model.


In this case, at $\phi=0.35$; 26.78% of the collisions have a $\Delta_s$ and the other one have a $\Delta_a$. Let's see if we can recover the curve with a probabilistic $\Delta$ model:
They are the same, it means that the fact that the even if the value of the kick depends on the free flight time (which creates correlation between velocities) it not important for the appearance of the fat tail.

### [23/01/23]

Here are better, coexistence phase at equilibrium:


At each density, I have a set of point with a streched config in the $y$ direction. We see that $p_x$ is less dependent on this streching. NICE!
Here is the minimum of these points for these two curves:

$p_x$ is always greater than $p$ NOT NICE!

### [30/01/23]

The straight line is Etienne's preditction for the critical pressure

Here, i've plotted the total pressure for the y spacing which corresponds to the lowest $p_x$

Both plot are a bit above the prediction but I think that it's ok...

_______

We have some issues on the spacing dependence though... We don't see a clear minimum of $p_x$ at a given spacing (here we are at $\phi = 0.849$):


and $p$ is non monotonous...


_______________

#### About the Synchro model
I've plotted the synchro and asynchro velocity distrib and they are not gaussian...


I don't think there is much we can do theoretically with this model. We would need to take into account the difference of velocity distribution and differentiate the two species in the theory. Which is kinda hard. We don't even know their temperature separately... In any case there is too much different things, it's too messy.

_________

### [03/02/23]

#### Modeling with two populations

We have a population of A(synchronized) particles and one of S(ynchronized) particles.

$P_B = 1-P_A$. The evolution of $P_A$ is:

$$\dfrac{dP_A}{dt} = W_{S\rightarrow A}P_s - W_{A\rightarrow S}P_A$$

* $W_{S\rightarrow A}$ the rate at which particles change from state S to state A:
$$W_{S\rightarrow A} = \omega(T)P_A$$ The number of collisions by unit of time time the probability that the other particle involved in a collision is Asynchronized. Technically, since we have two population the collision rate depends on the state of the particle. But let's forget that!

* $W_{A\rightarrow S}$ the rate at which particles change from state A to state S. Let's assume that the number of collisions $n$ a particle undergo between 0 and $t$ follows a poisson process. Then $$P(n) = \dfrac{(\omega t)^ne^{-\omega t}}{n!}.$$
The probability of having no collision is:
$$P_e = e^{-\omega t}$$
Particles A that don't collide during $\tau_s$ move to S, so: $$ -W_{A\rightarrow S}= - \textrm{proba of no collision}\times \textrm{rate of transfer }A\rightarrow S = -e^{-\omega \tau_s}\dfrac{1}{\tau_s}$$
If particles don't collide and can freely move to S because the frequency of collision is very low compared to the frequency of synchro, we recover $-1/\tau_s$ at the 0th order. The first order approx gives: $-W_{A\rightarrow S}= -\dfrac{1-\omega\tau_s}{\tau_s} = -\dfrac{1}{\tau_s}+\omega$. We see that collisions effectively reduce the rate at which we go from A to S (because of the term $+\omega$)

Now that we know the equations that describe the change of population, we have to couple it to the equation that describe the energy of the total system because the frequency of collision depends on it. Again, we will just look at the energy of the whole system and not the energy of each population separately. It would be waaaaay better, but I don't really know how to do that. On top of that, if we did that; we would have 3 the collision frequencies between AA, between SS and (harder) between AS...

The evolution of T is given by:
$$\dfrac{dT}{dt} = (P_A^2 + 2P_AP_S)\langle E' - E\rangle_{\Delta_A} + P_S^2\langle E' - E\rangle_{\Delta_S} - 2\gamma T= (P_A^2 + 2P_AP_S)G(\Delta_A) + P_S^2G(\Delta_S)$$

$P_A^2 + 2P_AP_S$ corresponds to picking at random either two A or one A and one S. NOW, the main issue, is again the fact that we consider only one temperature: $T$ while $G_A$ should depend on $T_A$ I think... But maybe it's not that much of an issue, results will tell us.

_______________

Here are the two coupled equations we have to solve to get the evolution of the population:

$$\dfrac{dT}{dt} = (2P_A - P_A^2)G_A(T) + (1 - P_A)^2G_S(T) \\
\dfrac{dP_A}{dt} = \omega(T)P_A(1-P_A) -\dfrac{e^{-\omega \tau_s}}{\tau_s}P_A$$

______

blue is theory and purple is simulation

$\Delta_S = 0$, $\Delta_A = 0.03$, $\tau = 1$

$\Delta_S = 0$, $\Delta_A = 0.03$, $\tau = 6$ WRONG SEE BELOW

-----

We can now use a "theoretical frequency" of collision which breakdowns at the transition:

In this case, the frequency of collision is: $$\omega(\phi, T) = 4\pi rg^+\sqrt{\pi T/m}$$

with $g^+$ the RDF at contact. Which can be found from the pressure, in 2D:

$$g^+ = \dfrac{Z - 1}{2\phi}$$ with $Z = PV/(NT)$ and
$$P = \dfrac{1 + \phi^2/8 - \phi^4/10}{(1-\phi)^2}$$
[Source here](http://www.sklogwiki.org/SklogWiki/index.php/Solana_hard_disk_equation_of_state)

As a side note, for very inelastic gas, we now that this will breakdown because the pressure is corrected by a factor $1+\alpha$. So the same formula can't hold for eq and out of eq.

_____________

#### Comparison with the other theory


Concerning this plot, there is an issue, $\tau$ seems to control the time to reach the steady state at the transition point. Which is verrry large. If I let the equation of population run for a longer time, I don't have a continuous phase transition anymore but a jump.

Here is the plot with a long final time for the population equation:

We see that the new and old theory are almost the same. For these parameters.

For some parameters the old one is even better:


I think that overall, this population model doesn't give better result than the other one. I'll have to add a $T_A$ and $T_S$ to population equation. Indeed, it doesn't seem to fail because of the breaking of molecular chaos or even maybe because of the non gaussianity but because each specie, at the transition has a very different behaviour from the other one that cannot be captured by the total temperature of the system.

Here is a naive try at a theory with two temperature but it doesn't work... I'll try again later

### [06/02/23]

#### Random $\tau_s$ or constant?

Instead of taking a constant $\tau_s$ we can take it from an exponential distribution $\tau_s \sim \dfrac{1}{\tilde\tau_s}e^{-t/\tilde\tau_s}$:

#### Explaining in more details the two temperature construction

Let's write the change of energy in an given interval of time large enough so that we can take the mean values because of the high number of collisions. In this interval there are $N_{aa}$ collision between AA particles, $N_{ss}$ collisions between SS particles and $N_{as}$ collisions between AS particles.

* The change of energy for AA and SS collision is easy since it doesn't mix the two species, we can just use the same kind of term as before: $\delta E_a = N_{aa}\langle E'-E\rangle_{T_{a}, \Delta_a}$ and $\delta E_s = N_{ss}\langle E'-E\rangle_{T_{s}, \Delta_s}$

* Now, there is an issue with AS collisions! Because, we start with one A particle and one S particle but end up with two A particles. So, the contribution of the collision to the energy goes to $E_a$ but we have to be careful to retrive the energy of the (back before the collision) S particle that is now (after the collision) an A particle to the S energy. Something along this way might work, we have two particles; one labeled by 1 which was in the state $a$ before the collision and one labeled by 2 which was in the state $s$ before the collision: $$\delta E_a = N_{as} \dfrac{m}{2}\langle(v'_{1, a})^2 + (v'_{2, a})^2 - (v_{1, a})^2 = N_{as} \langle E' - E\rangle_{T_{as}, \Delta_a} +N_{as}\dfrac{m}{2}\langle (v_{2, s})^2\rangle_{T_s}\\ \delta E_s = -N_{as}\dfrac{m}{2}\langle (v_{2, s})^2\rangle_{T_s}$$
Now, what do I mean by $T_{as}$? Recall that:$$\langle E' - E\rangle_{coll} = m\Delta^2 + m\langle|\vec v_{12}\cdot \hat r_{12}|\rangle_{coll}\alpha \Delta - m\langle|\vec v_{12}\cdot \hat r_{12}|^2\rangle_{coll}\dfrac{1 - \alpha^2}{4}$$ with $v_{12}$ the relative velocity between the two particles involved in a collision. Here, both particles have a gaussian distribution, but one is synchro (so folowing a gaussian) with temperature $T_s$ and an asynchro with temperature $T_a$, since we are just summing random variable, the distribution of $v_{1, a; 2, s}$ is a gaussian with a temperature $T_a + T_s$ so $T_{as} = \dfrac{T_a + T_s}{2}$ (factor 0.5 is there for consistency with my notation).

* Now, let's convert this $\delta E$ into nice differential. It is simply a matter of changing $N$ into a frequency. Since we are now dealing with two species, the frequency of collision is not simply $\omega(\phi, T) =\Omega(\phi)\sqrt{T}$, we also have to take into account the probability of each species. We will suppose that the positions of our particles are not correlated to their specie (which is definitively false :), and might play a very important role but I don't know how to go further...). So the frequency of collision between AA and SS species is easy: $$\omega_{SS} = P_S^2\omega(T_S)\\ \omega_{AA} = P_A^2\omega(T_A)$$
For $AS$ collisions, we will suppose that the frequency depends on the mean relative velocity between particles, so: $$\omega_{AS} = 2P_AP_S\omega(T_{AS})$$

* Packing all of these together, we get: $$\dfrac{d T_a}{dt} = \dfrac{\omega_{AA}}{2}\langle E' - E\rangle_{T_{A}, \Delta_{A}} + \dfrac{\omega_{AS}}{2}\left(\langle E' - E\rangle_{T_{AS}, \Delta_{A}}+T_S\right)-2\gamma T_A\\ \dfrac{d T_s}{dt} = \dfrac{\omega_{SS}}{2}\langle E' - E\rangle_{T_{S}, \Delta_{S}} - \dfrac{\omega_{AS}}{2}T_S-2\gamma T_S$$

We are not done yet! We have to write the evolution equation for $P_A$.

* $S\rightarrow A$: Each time a particle A collide with a particle S, the S particle become to A (so one half of them). So we just have a term $\left.\dfrac{d P_A}{dt}\right|_{S\rightarrow A}=\dfrac{\omega_{AS}}{2}$ This corresponds to the same term as the one found for the previous equation with only one temperature.

* $A\rightarrow S$: At a frequency $\tau_s$ each A particle will become S **unless** it is involved in a collision. Particles $A$ will be involved in a collision at a frequency: $\omega_{AA} + \omega_{AS}$, so the probability of having no collision between a particle A and an other particle in $\tau_s$ is: $e^{-(\omega_{AA} + \omega_{AS})\tau_s}$. So: $\left.\dfrac{d P_A}{dt}\right|_{A\rightarrow S}=-\dfrac{e^{-(\omega_{AA} + \omega_{AS})\tau_s}}{\tau_s}P_A$

* Which gives us: $$ \dfrac{d P_A}{dt} =\dfrac{\omega_{AS}}{2} -\dfrac{e^{-(\omega_{AA} + \omega_{AS})\tau_s}}{\tau_s}P_A$$

Doesn't work :)


I've issues with the frequency of collision. $\omega$ is twice the number of collisions that happens in a system divided by 2 (ie: $\omega = \lim_{t\rightarrow \infty }N_{coll}/t$ with $N_{coll}$ increasing by 2 at each collisions). Since it really is the number of collision ONE particle feels during a second (so doublke counting). Sometimes I'm using $\omega$ when maybe I should use $\omega/2$... I don't think so, here is my rule of thumb:
When the computation is about a collision as an even in the system, I use $\omega/2$ (like for the energy change). When the collision is about a change of property of a GIVEN PARTICLE, I use $\omega$ because this particle will undergo collisions at a frequency $\omega$.

------

### [03/04/23]




LAMMPS SIMUILATIONS

comparison between 1k and 5k

__________________________

### [13/04/23] Swimming in a sea of unknown... and drowning

I ran multiple tests with the data at hand. But TLDR: I need better data, simulations are running and the cluster is being stormed as it should.

Some papers talk about finiote size analysis using only the metastable/transient value of the order parameter such as this one: https://arxiv.org/pdf/cond-mat/0309165.pdf#page=11&zoom=100,334,902 They also use binder cumulant but it did'nt work extremly well for me!

Other articles:
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.48.1710

(same curves as us:
)

For now, every state point has been run just once and no average is not (if not stated otherwise).

#### Time to reach the absorbing state:

legend is number of particles.

linear scale

log scale

In linear scale, we can think that a bit below 0.1505 all curves collpse to the same time to reach the steady state and above, the curves are highly dependent on N. But, in log scale it doesn't look like it's the case (time to reach seems proportional to $\sim e^{\phi}L^\mu$). We need more simulations.

In any case, I think it is interesting to check the time to reach the steady state according to $\phi$ since this time is related to the correlation lenght of the system. Indeed, the system can reach an abosrbin stat when $\xi\simeq L$. Since $\xi$ follows a power law with N and packing fraction, we might also observe a power law behaviour for the time to reach the dead state if the transition is of second order!

_____
Also, coherent:


I think these figures already show that we are in a second order phase transition. With first order phase transition, the correlation lenght doesn't scale as a power law with $\phi$ and $N$. Here we really see a strong dependence on the time needed to reach the absorbing state with both parameters...

#### Fluctuation and E

E might be staying almost constant while std(E) seems to decreases (as $N^{-0.5}$, blue line). In fact, for a usual phase transition we expect to find $E(\phi_c) - E(phi)\sim N^\mu$. If $\mu<-0.5$ we should be good. In any case, here, it seems that the separation between E and std(E) is growing
Here E/std(E) seems also to be growing:

Especially when we don't look at extreme points close to very active state: 

Now, there is the issue that the suceptibility (size of fluctuation), $\chi = N(<N_a^2> - <N_a>^2)\sim |\phi - \phi_c|^{-\tilde\gamma}$ with $N_a$ the number of active particles (which is almost the same thjing as $E$) is supposed to diverge closer to the transition point. This doesn't seems to be the case... Perhaps due to $\gamma$. But for this I need more data. For now, it looks like that, mostly, the fluctuation are not extremely gaussian like. If it was: $Var(N)\sim N^{-1}$ and $\chi\sim N^0$ implying that we couldn't observe a divergence of $\chi$ :

N = 20k

N = 18k

We see a small bump at low value, this is caused by "failed" abosrbing state. A configuration that almost lead us to an absorbing state where the dyanmics was really slow and then, an other particle was touched and the dynamics could start again. THIS is propably the key to observe a diverging fluctuation size.

Maybe, the critical point is way lower than what I expected.

#### Evolution of number of active particles

At least this figure is ok :^ ⁾. We see two trends and eventually the critical packing fraction!

### [18/04/23]: SIMPLE DELTA MODEL

A convincing argument is the one of increasing $\chi = N(\langle\rho_A^2\rangle - \langle\rho_A\rangle^2)$ ($\rho_A$ the number of active particles. For some reason, I have better statistics with this quantity than with $E$) as a power law close to the critical packing fraction. The exponent derived from a curve fit is 1 close to the mean field one for most critical system (strangely) and the critical packing fraction derived from the fit is close to the one derived from the mean free path argument.



We also observe an increase of the time needed to reach the steady state:
It is extremely abrupt

And not a power law. I think that's because, in fact, all these points are above the critcal packing fraction which should be around $\phi=0.148$, a point at which we would probably see a collapse of all curves to almost the same value.

We also see a similar behaviour than the directed percolation


### [10/05/23]

| CONTINUOUS | DISCONTINUOUS (I hope so) |
| --------------------------------------------- | --------------------------------------------- |
|  |  |
|||
|||
|||

Discontinuous is jumping at a way higher ratio of E/std(E) to the absorbing state but simulations are lasting a different amount of time (I'm not comparing equal collision time).
For the double delta model. We see the whole system being active
But, now, I'm worried that both transitions are continuous....

### [06/07/23]

Memo: only the rule A->S and without the rule S->A. In this case, d N_A/dt = W_(A->S)N_A and so, by looking at the exponential decay of N_A(t) we can see which W is the right one:
phi = 0.2
phi = 0.5

checking poissonianity: 

### Real EDMD model mess :' )
The model I was simulating was not the one we were thinking.

Before, each collision would make the particle async (so S-S coll would make both particles transition to A).

While with the "correct" model, collision would async the particle only if it was a A-S collision (for S-S collisions, particles would still be sync after the collision).

Now, if we simulate this new model, we obtain:


Where the first jump is a stable wavy dynamics:
{%youtube D4OKufsfGqg %}

Which was observed in the wrong model, but it was not stable and was a prelude to the dying of the system:

{%youtube 6pE4F7OpRco %}

Moreover, in the wrong model, this wave dynamics observed at the critical packing fraction was a state that was attained after a lot of time while here, for the values tested, it seems pretty easy to get to it (the system falls immediately in this config in almost all cases tried).

Interestingly, by lowering the coefficient of restitution we cannot get rid of the waves but we can smooth out the energy. Above, we were at $\alpha = 0.95$, here is $\alpha = 0.75$ and $0.6$:

| Column 1 | with theory | without theory zoomed|
| -------- | -------- | -------- |
| res = 0.75 |  |  |
| res = 0.6 |  (theory predict continuous) |  |

high $\alpha$ messes up the theory because of high heterogeneity and non gaussianity. Note also that We are at a small $\tau_s$. I think that at high $\tau$ the theory might be better.

And the dynamics really looks like what we observe with lammps:

LAMMPS above $\phi_c$:
{%youtube mDmBXJ3uYYw %}

LAMMPS below $\phi_c$:
{%youtube wEX2yUNsaGc %}

EDMD res 0.6 above $\phi_c$:
{%youtube EOA_7n7Mkn8 %}

EDMD res 0.6 below $\phi_c$ (we dont see it dying because it takes a lot of time):
{%youtube A27c428eWqs %}

BONUS: EDMD res 0.75 above
{%youtube znPcFkOqaG0%}
___________________________________

With a random $\tau_s$ event at high $\alpha$ we can smooth out the energy (we still see a change of trend I think):

$\alpha = 0.95$

But the waves are still very much there (in any case I think we understood that these waves were something physical but they are kind of hidden by the small effective coefficient of restitution in lammps).

______

The issue with small coefficient of restitution is that it messes up the theory of the simple delta model:

$\alpha = 0.75$:


(N = 5k, explaiinng why the MFP argument might be a little wrong)

I'm also a bit afraid that for a small coeff, the transition is second order, but to be sure, I have to do longer simulations on a wider range of N and $\phi$:


We could mix low $\alpha$ and random time and it seems that the transition is similar in scaling but this might be caused by a weird change of the energy at the beginning with a random time that could hide from us the true steady state :
{%youtube TFPkD4ORAUA %}

### old model:

Maybe the wrong model is not that bad since the phenomenology of lammps is recovered with it (just need to justify S-S -> A-A)

For the correct model:
**A-S** -> **A-A**, **A-A** -> **A-A** and **S-S** -> **S-S**.

$$\dfrac{d P_A}{dt} = -we^{-\omega \tau}P_A + \omega (1-P_A)P_A$$

A strange thing is that the theory of the correct model with the data of the wrong model agree pretty well:

#### data of wrong model + theory of correct model

For the wrong model:

**A-S** -> **A-A**, **A-A** -> **A-A** and **S-S** -> **A-A**.

$$\dfrac{d P_A}{dt} = -we^{-\omega \tau}P_A + \dfrac{\omega}{2} 2P_AP_S + 2\times\dfrac{\omega}{2}P_SP_S= -we^{-\omega \tau}P_A + \omega (1-P_A)$$

#### Theory of the wrong model with the data of the wrong model:


The theory for the wrong model predicts a continuous transition (whatever $\alpha$ is) while with the wrong model we observe a discontinuous transition.

The steady sytate of the wrong model is given by:

$P_A e^{-\omega\tau} = (1 - P_A)\Rightarrow P_A = \dfrac{1}{1 + e^{-\omega \tau}}$ or $P_A = whatever$ when $\omega = 0$
Which can be 0 only when $\omega= 0$ otherwise it's minimum is 0.5 (For the other model: $P_A = 1 - e^{-\omega\tau}$ which can continuously go to 0 with $T \rightarrow 0$)

| | wrong model simu | wrong model theo prtediction |
| --- | --------------------------------------------- | --------------------------------------------- |
| E |  | |
| P_A |  |  |

We can also try to see if the measured steady state $P_A$ in the simulation is close to the theoretical $P_A$ but using the frequency of collision measured!


Not really

### [20/07/2023] 3 AM doubt over kinetic theory

In the kinetic theory, we use: $\Delta(v_1, v_2) = \Delta \Theta(v_{12} - l/\tau)$ while we should be using: $\Delta(v_1, v_2) = \Delta \Theta(\max(v_{1}, v_{2}) - l/\tau)$! Fortunately there's not much difference between the two when we average. Here $$\langle\Delta\rangle_{coll} = \dfrac{1}{\sqrt{\pi T}}\int dv_1dv_2|v_1 - v_2|\Delta(v_1, v_2)f(v_1)f(v_2)$$

With $f$ gaussian.


Unfortunately, I don't know why

### [28/08/23]: A universality class different from the manna model?

Cates model is believed to belong in the CDP/Manna universality class. Here are the results for our model:

With perhaps a slighlty better $\phi_c$
It's hard to tell which critical exponent is better. I would argue that $\beta = 0.7$ is better.

_____
For the evolution of the temperature with time at the critical packing fraction, the agreement between the manna model and ours seems pretty robust:

NONETHELESS, with the fluctuation, we are far from the manna model:


And this can't be explained only by a wrong $\phi_c$. Perhgaps using $T$ and not the numbver of active particles causes this issue! Moreover, it is known that by adding activity to dead particles we can change the universality class https://journals-aps-org.ezproxy.universite-paris-saclay.fr/pre/pdf/10.1103/PhysRevE.105.L032602 however, our model doesn't fall in this new universality class. But the idead that activity can change it is interesting. Especially, because in our kinetic theory, we also have term non analytic in $T$ as they do.

### [04/09/23] Where to go next?

- **Universality class of the continuous transition?** The simulations would not be that long I think... But do we care? An interesting point is that we seem to be closer to the Manna universalitry class than the universality class investigated in *Absorbing phase transitions in systems with mediated interactions*. The difference between their model and the cortes one is that, in their case, dead particles have a small "activity", they can move a bit. Which in turn adds a term in $\sqrt{order~parameter}$ in the landau type free functional which lets them escape the manna model/Directed percolation universality class (which is the one of the cortes model). I dopn't think this is what we have in our case, but since we don't really have "dead/active particles" in our model (because we can't have overlap) I don't know. I like the fact that our model with kinetic theory also contains non analytic term. Which might hint to a new universalirty class (the kinetic theory breaks down at the critical packing fraction, but I think that should be ok).

- **Kinetic theory with spatial and/or time dependence**. What Olivier proposed, with Andrea Puglisi. Seems cool! When I realized that my code was not doing the right thing, I tried a lot of configuration and a lot of parameters. There's really a zoo of spatial structure and weird evolution. For example very strong waves (stable) or an initial decrease of the energy with the system almost completely dying to in the end, REHEAT LIKE HELL (like 3x the minimal energy observed). So there's a lot to unveil, but maybe too much :) It seems very complicated.

- **Granular osmosis**. Frank was talking about checking the existence of a chemical potential out of equilibrium I think, using the osmosis setup. Andrea talked about a presenbtation at a conference abouyt granular osmosis and a theoretical prediction with onsager relation I think? But I don't really want to mess with LAMMPS that much. Neither to look at a perfect mapping between lammps and EDMD to have the exact same phenomenology as Gabriel was looking at. Maybe pure EDMD would be enoughg, at least at the beginning. I like also the idea of extednging a bit the equilibrium theory of osmosis to out of equilibrium situation with kinetic theory or usual langevin/Fokker plank equation.

- **Crystallization**. Last year, about this project, we talked about an out of equilibrium phase diagram of crystallization but it seems a bit unimportant and/or boring.

- **True cristallization/Breakdown of Mermin–Wagner theorem** but investigated with our model at high density (from Andrea: *Two-Dimensional Crystals far from Equilibrium*).

- **Parallel EDMD** At some point we very quickly discussed about letting me implement a parrallel EDMD (taken from a paper) which could be fun... Maybe not fun, but interesting :)

### [06/09/23]: Mapping Hertz mindlin to hard spheres with ML/array of values

The idea is to try to model the model given by LAMMPS in EDMD. We need to find $f(...)$ such that an instantaneous collision given by the usual hard sphere collision rule:

$$v_{i, j}' = v_{i, j} \pm f(...)$$

reproduces as much as possible the dynamics onbserved with lammps. We also need to incorporate the angular momentum $\omega$. As a first try, we forget $\omega$ and place ourselves in 2D.

In this context, $f$ depends on $\vec v_i - \vec v_j$ and on $\vec r_i - \vec r_j$ (as well as on the physical parameters of the material used for the beads, blablabla). We can perform the LAMMPS simulation for $\vec r_i - \vec r_j // \vec u_x$ and by rotation infer the result of $f$ when in EDMD the collisions won't be parrallel to the x axis (since $|\vec r_i - \vec r_j| = cst$).
_____
In lammps, we simply initialize the system using two particles Black and Blue at contact:
The blue one is at stop, and the black one is moving toward the blue one. For a large range of velocity of the black particle, we see how the velocity of the blue particle changes. Its velocity change is $f$ for a given $\vec v_i - \vec v_j$ (which in this case is simply the velocity of the black particle) since its initial velocity is 0.

| INITIAL | FINAL |
| -------- | -------- |
| | |

Here is a plot of the experience. In black we have curves of $v_x$(left pannel) and $v_y$ (right pannel) for the black particle and the corresponding velocities for the blue particle.

(it gives very cool figures:

)

Here, $|\vec v|$ was fixed and only the angle of collision changed. We have to do this thing with a ton of different $|\vec v|$. This done (not done at the time of writing :) :) ). We can numerically map $\vec v_i - \vec v_j$ to $\vec f(\vec v_i - \vec v_j)$ using a neural network or a table of values. Let's go for the neural network because it's trendy.
Here: $\vec v_i - \vec v_j=\vec v_{black}(t = 0)$ and $f=\vec v_{blue}(t = \infty)$. As a sanity check we can verify that $v_{black}(t = \infty) = v_{black}(t = 0) - f$

the error is defined as $\vec {error} = \vec v_{black}(t = \infty) - \vec v_{black}(t = 0) - \vec f$

_______

After training, we get a prediction of the collision:

Somejhow, it even manages to generalize this is for a velocity that is x100 larger than the maximum velocity on which it was trained:
The curve looks correct...

We simply have to use this prediction in our EDMD algorithm (and don't forget to rotate the vector so that $dy = 0$):


### [07/09/23]

An issue is that the algorithm reach a non trivial steady state while it should be purely dissipative.
Here is the ratio of the post collision with the pre collision velocity WITH LAMMPS we see that the function looks like it diverges with $v\rightarrow 0$. Physically, we know that it should go to 1 at 0.

The same happens for slippy collisions:

After learning on these data, here is the prediction of the neural network trained for velocities as small as 0.01:


________
Off course training it on smaller values allows us to push smaller energy for the steady state.

### [08/09/23] quasi 2D

For quasi 2D, it's not that harder (so 3D). By rotational symmetry, we can keep the initial model which was trained on $\vec r_i - \vec r_j // \vec e_x$ with $\vec v_i - \vec v_j \in \vec e_x \vec e_y$, the plane spanned by $\vec e_x$ and $\vec e_y$. We now have to rotate every vector so that the input we give to the neural network corresponds to a data on which it was trained. Then we rotate back the result to the initial configuration of the collision:

First we need to rotate $\vec r_i - \vec r_j$ so that it becomes parrallel to $\vec e_x$:

We define: $$\theta_1 = \mathrm{atan2}(y_j - y_i, x_i - x_j)~~\mathrm{ and }~~\theta_2 = \mathrm{acos}((z_j - z_i)/|\vec r_j - \vec r_i|)$$
And then apply some rotations:

$$R_y(\pi/2 - \theta_2)R_z(-\theta_1)$$
To our vectors distance and relative velocity. the $z$ rotation kills the $y$ components of the vector distance and the $y$ rotation gets rid of the $z$ component which leaves us with a vector distance // to $\vec e_x$.

The transform velocity is now: $\tilde{dv}$

Then, we need to rotate one last time around the $x$ axis so that the relative velocity is $\in \vec e_x \vec e_y$. This is achieved by rotating $\tilde{dv}$ by an angle: $$\theta_3 = \mathrm{atan2}(\tilde{dv}_z, \tilde{dv}_y)$$.

After that, we give $\tilde{dv}_x$ and $\tilde{dv}_y$ to our neural network. Then perform the opposite rotation on the to have the resulting velocity change. And tada (with a small vibration):


_________

The algorithm is not very stable because it doesn't like almost perpendicular/slippy collisions. Moreover, there's still no angular momentum.

But in any cae, by adding small friction at roof collision (for example: $v_{x, y} -> 0.9 v_{x, y}$ at collision with roof or bottom plate). We can (off course) get an absorbing state:

Funnily enough, we see breathing like pattern, like what I observed with EDMD (but not with lammps):
And this pattern, I think, give an energy curve similar to the one found in EDMD with two different regimes, one after the transition and a second one way after it...

____
Wat's next is adding angular momentum, which is straightforward (I think, but kind of long...)

### [11&12/09/23] ELUCIDATE THIS HYPERUNIFORMITY QUESTION ONCE AND FOR ALL

I learn from Andrea that there is already an analytical mapping of a granular soft sphere model to a hard sphere in **Model for collisions in granular gases**. So the ML is paused.

______

More on manna model

______

#### HYPERUNIFORMITY

Every paper I could read basically said that Manna model universality class implied hyperuniformity close to the transition:
(in **Hyperuniformity of critical absorbing states**)

In these systems, the critical region is believed to be hyperuniform with $S(k\rightarrow 0)\sim k^{0.45}$. This extends slightly above and below $\phi_c$ for a small range of $k$. A bit above $\phi_c$ we can observe $S(k)\sim k^{1}$ for a small $k$. It is nevertheless believed that $S(k\rightarrow0)=cst$ above $\phi_c$ (by some, the author of the paper finding the $S \sim k$ says that it's the $S\sim k^{0.45}$ that is the transient..):
(in **Noise, diffusion, and hyperuniformity**)

See for example:
(in **Hyperuniform Density Fluctuations and Diverging Dynamic Correlations
in Periodically Driven Colloidal Suspensions**)
(note the connection between fig a. and the mean free path argument)

or
(in **Hyperuniformity of critical absorbing states**)

Here is what i got with my system:

It's similar if we don't mind the small k which are due to (I think) averages. In any case, it's expected to be hyperuniform until a cutoff $k$ that diverges as $\phi\rightarrow\phi_c$.

Here are results for $\phi$ closer to $\phi_c$
The legend is $\phi/\phi_c$. It looks quite hyperuniform..

Strangely, in our case, it seems that we have a transition from a regime 0.45 to 1 even if $\phi < \phi_c$

Let's see different size:


I can't tell if the bump at the end is due to poor averaging or if it's genuine.
At least with $\phi$ closer to $\phi_c$ we can get a longer $k^{0.45}$ behavior. So I think, this is geniune!!!!


______

STRANGELY ENOUGH: it seems that our model (with or without drag) exhibits hyperuniformity in the active phase:


This is the case also in some other model but only when the center off mass is conserved by the dynamics (which is not the case of our dynamics due to the drag, but maybe it is insignificant at high densities because the main source of energy change is collision...):
(in **Noise, diffusion, and hyperuniformity**)

Here is what i got:



Here is the eq to compare:

If the active state is indeed a hyperuniform state, we can expect that the solid is a true crystal with long range order. Even without drag.

_________
I'M pretty much conveinced that we are indeed in a true crystal at high friction.

This is the curve for MSD starting from a crystal (almost perfect, we'll come back on that later) with strong drag!

We see that most MSD plateau at the same value except a few. These ones, are not due to phonons but to defects diffusing. My initial crystal is not perfect and sometimes we can a defect diffusing (this is now fixed, it was expadning at $t = 0^+$) until the steady state is reached (the plateau):

For comparison, here is the MSD vs time of an equilibrium system:

(also with defects problem)

_____
DONE:

## unfortunate sunday evening correlation:

__________

[19/09/2023]
Ran Ni model:

We also observe an exponent closer to 0.2...
Fluctuations, really happen on a huge time scale. That might be caused by a bad sampling:


_____ Puglisi article:

For example: https://arxiv.org/abs/1103.0166

His setup:

He wants to model his system with a langevin equation:

$\dot v_i = -\gamma_{plate} v_i + F_{coll} + \eta_{plate}(t)$

Where, for now, the noise is only due to the shaker $\langle\eta_p(t)\eta_p(t')\rangle = 2T_{plate}\gamma \delta(t - t')$.

($_{plate} = _p$)

We'd like to model the collisions as a noise because it's hard to deal with $F_{coll}$.
But, there is a profound difference between the two type of noises. The plate gives noise as a normal thermostat. But, collisions, if they are to be modelled as noise, must conserve momentum. As in model B, we take the divergence of the noise instead of simply a noise. For example in fluctuating hydrodynamics:

with the noise chosen to satisfy the fluctuation dissipation theorem:
Which in turn leads to:

We see that the noise for the velocity field depends on $k$.

(more sources:https://arxiv.org/pdf/cond-mat/9706079.pdf https://arxiv.org/pdf/cond-mat/9810251.pdf)
We could continue with the hydrodynamic equation (it was the way Ran Ni choose, but they don't have an equivalent of the platge noise, they only have the drag). Instead let's go back to our simple langevin equation and add this noise (now, we talk about the transverve velocity because it's better describe by the above langevin equation
):

$\dot v_\perp = -(\gamma_{plate} + \nu k^2)v_\perp + \eta_{coll}(t) + \eta_{plate}(t)$
with $\nu$ the viscosity. Inspired by the hydrodynamic equation and the Fluctuation Dissipation Relation:

$\langle\eta_c(k, t)\eta_c(-k, t')\rangle = 2T_{c}\nu k^2 \delta(t - t')$
with $T_c= T_g$ the granular temperature.

Looking at the $<v^2>$ we obtain: $$\langle v_\perp(k)^2\rangle = \dfrac{\gamma T_p + \nu k^2 T_c}{\gamma + \nu k^2}$$
In the case of non dissipative collisions: $T_c = T_p$, we recover the equipartition theorem.
This shows that there is a transition from from $T_p$ at low $k$/large scale to $T_g$ for high $k$/small scale/collision scale.

A similar expression can be found for the density field. If the transverse velocity field follow a free overdamped langevin equation. The density field follows a underdamped harmonic langevin equation
From this equation we find the same kind of expression for $\langle\delta \rho (k)\delta \rho^* (k)\rangle$ as for $\langle v^2\rangle$.
We see that with a rough vibrating plate, we can't have hyperuniformity.

### [21/09/23]

For EDMD. I'm now sure that we indeed have long range order with delta + drag.
Top is non eq and bottom is eq in crystal.

__________

For lammps, I'm trying to find a set of parameter that would give us hyperuniformity in the dense state. It proves to be harder than what I expected. Since we need strong drag from the roof we need high amplitude or high frequency and small height.
It's pretty easy, with these range of parameters, to find something that develops strong inhomogeneities:

zoomed:

Strangely, the steady state reached is pretty energetic (the energy doesn't go to 0).

### [22/09/23] && [25/09/23]
Simply to be sure. I check with pure $\Delta$:
AS expected, no violation of Mermin Wagner !

_____

Interestingly, the MSD behavior for the MSD is more or less the same, whatever the density is!


__________

I did a scan of parameter with lammps to see if I could observe hyperuniformity in the dense region. Here is a plot of $S(q)$ for different frequency, amplitude and h at $\phi = 0.75$ for 100k particles.


No parameters lead to hyperuniformity.

It's not a perfect indicator but here is a plot with every simulations:


A point corresponds to a simulation and its color corresponds to value: $S(min(q))$. Black is high value (Strong clusterization) and white is small. We see that a high frequency facilitate clustering. We also see that we're better off with a small height.

BUT, the real issue is that particles are synchronizing at these densities! I thought it wouldn't be a problem at densities this high and thus decided to look at high frequencies. We can look at the energy. It's decreasing and leading us to clusters.

I will now restrain myself to zones known to be asynchronized.

here is an example in the async zone:


not very promising.

### osmosis proof of concept with EDMD

{%youtube tJjDKlO0Zyk %}

### Osmosis computation for the difference of density.

In the steady state with no flux, the solvent density is:
$$\dfrac{\partial (D(x)P)}{\partial x}=\dfrac{f}{\gamma}P\tag 1$$

With $D$ the coeff of diffusion **of the solvent** that depends on $x$ because there is a gradient of temperature and $f$ the force felt by the particles because of the wall and collisions between the solute and the wall. $\gamma$ is an effective coeff of friction.

$$f = -V n_{solute}\dfrac{\partial U_{barrier}}{\partial_x} = V T_{solute}\dfrac{\partial n_{solute}}{\partial_x}\overbrace{=}^{????}V T_{solute} \Delta n_{solute}\delta(x -L)$$

where we supposed that there is no solute on the left. The barrier is at $L$ (box size $2L$x$L$). We assumed an equilibrium distribution for $n\sim e^{-U/T_{solute}}$. We add a $V$ so that the number of solute cointributing to the force is: $V n_{solute}$ (hope it doesn't depend on the concentration of solvent)

Equation 1 has the solution:
$$P(x) = \dfrac{P_0}{D(x)}e^{\displaystyle{\int_0^{x}\dfrac{Vf}{\gamma D(x)}dx}}$$

With $P_0$ a normalization factor that I'm too lazy to compute.
Which gives:

$$
P(x) = \left\{
\begin{array}{ll}
\dfrac{P_0}{D(x)} \textrm{ if } x<L\\
\dfrac{P_0}{D(x)}e^{\dfrac{VT_{solute}\Delta n_{solute}}{\gamma D(L)}} \textrm{ else}
\end{array}
\right.
$$

At equilkibrium, we have equal temperature everywhere:

$$\int _L^{2L - a}P(x)dx-\int_0^LP(x)dx = \Delta n_{solvent}$$

With a taking into account excluded volume due to solute. I have to think. This excluded volume should give rise to diffusion but i'm not sure that is does