# Dissipative EDMD (Raphael)
###### tags: `M2`

---

---

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

EF, AP, FS, GF
## MEMO to do list
* Validate $\Delta$ and exponential model with lammps simulations doing the statistics of the collision. Do we measure a $\Delta$-like/exp-like rest coeff. in simulations where we observe continuous/discontinuous transition?
* Derive the rate of dissipation/energy injection G for the exponential model
* Try to go beyond maxwellian assumption of the pdf. We can actually measure it in simulations (Rapha did it) and experiments. We have probably a kind of Gaussian more peaked in zero (something like eq. 5, 7, 13 of https://iopscience.iop.org/article/10.1209/0295-5075/102/14002/meta?casa_token=A2Jpy0Qt5c8AAAAA:UeK72-gSYFQCEqagS_9UUepveFuk2JC6EG7-BXxbaa-vpzlWkrxa6uLOpiiiMMqeaFRM1rAlu-Q )

### Note about Andrea's comment

As Andrea told me, there is in fact, whatever $\phi$ an other solution to $G_{tot} = 0$ than $T = 0$ (I missed it because I wasn't looking at $T$ very small, shame on me). So I calculated $T_{SS}(\phi)$ numerically from $G_{tot} = 0$:


The shape looks good (particularly the linear tendency at $\phi>phi_c$) but the critical point and the values are wrong. I used the ideal gaz free mean path, it would be different with the real free mean path (I have to test, but I've lost my table of mean free path so it's a work for another day!)

# Dynamical phase transition

### Theory?

Naively, we could try to obtain the critical packing fraction this way:

We suppose that at the transition almost all particles are stopped. In order to activate the system, one particle has to be able to travel the mean free path of the system. Since all particles are almost dead, the velocity after the collision is $\Delta$ which means that the particle travel, in a infinite amount of time a distance $\Delta/\gamma$. Thus, at the transition we expect $\Delta/\gamma \simeq l$ with $l$ the mean free path.

For particles at rest (except one), the mean free path in 2D is given by:

$$l = \frac{1}{4r\rho} = \frac{\pi}{4}\frac{r}{\phi}\Rightarrow\phi_c\simeq\frac{\pi}{4}\frac{\gamma}{r\Delta}$$ with $\rho$ the surface density of particles.
It gives the good order of magntitude since we have $\phi_c = 0.17$ for $\gamma = 0.01$ or $\phi_c = 0.52$ for $\gamma=0.03$. But it's not that great... I guess we would need to take into account fluctuations of density in order to correctly predict $\phi_c$... (for future me: more infos on report)


color plot = simulations

I also tried using an approach based on [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). Starting from the rate of dissipation/energy injection $G$ (they say that it's in 1 over area unit divided by time but I don't think so for the area part?) and the fact that in a Steady State, $G = 0$. They give:

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

with $\omega = 2\rho\sigma\chi\sqrt{\pi T/m}$ the collision frequency which I would correct in: $\omega = <v>/l_{dyn}$ with $l_{dyn}$ the mean free path when all particles move (which is not the case. A factor sqrt(2) doesn't change much the result anyway) $l_{dyn}=l/\sqrt{2}$ and $<v> = \int_0^\infty m/T v^2e^{-mv^2/(2T)}=\sqrt{\pi T/(2m)}$ giving:

$$\omega = 4r\rho\sqrt{\pi T/m} = 4\phi/r\sqrt{T/(\pi m)}$$

I've a factor 2 on top compared to their result and I don't know what $\chi$, the static pair correlation function at contact is (in any case it wouldn't change much the result).

Since in our simulations, we have a $\gamma$:

$$\frac{d T}{dt}= -\frac{1}{2}m\gamma v(2v) = -2\gamma T $$

In total we look for: $G_\gamma=G + 2\gamma T=0$ at the steady state for a given $\phi$.

I search the first $\phi$ at which a solution different than 0 exist and I find that:


The derivation assumes a maxwellian distribution which is definitely not what we expect at the transition. We can make the predicted $\phi_c$ closer by changing the treshold on $E$ on the simulation curve to a lower value. But the shape is not good and in any case, the treshold, in this case, doesn't make any sense (visually, the transition doesn't happens at this value of $\phi$).


### The critical packing fraction using a better mean free path

The mean free path used at the top of the document (**for the first derivation**) supposed that the gaz was rarefied. There are ways to derive a better one, but most articles do this in 3D and I need to study a bit (a lot...) more to do it in 2D.

Nonetheless, we can find it numerically:

and using it (with $l = \sqrt{2}l_{dyn}$), we have a very very good agreement between "theory" and simulations since the two surfaces almost perfectly overlap.

Of course there is a treshold on E, so this is a bit arbitrary, but still... Here is a difference between the 2 curves for one of the best treshold on E (that makes sense, since the treshold is really small):


Now, in order to refined the analysis, we would need to take into account correlations.

### Different "equilibration" time

### Auto correlation

The velocity autocorrelation is independent on $N$ and dependent on $\gamma$ at the transition. We can't use it to find a diverging time scale. Or so it seems.


### energy

The determination of $\tau$ is based on the characteristic time necessary to reach energy "equilibrium" for the $\Delta$ model (this is nice, because it didn't work with the exponential model). As expected, $\gamma$ exponentially kills the time divergence
The points higher than expected are caused by an initial energy already close to the stationary value. The determination of $\tau$ in these cases is hard. Some funky things are happening at high packing fractions.

### Scaling


We obtain the same results with a MSD treshold based approach
### Questions

In fact, it looks like the $\Delta$ model could reproduce what we observe with the 2nd order phase transition observed numerically (by eyes, the dynamics looks the same and the curve might be also smooth at the transition). A mapping between the two models shouldn't be extremely hard:

$\gamma$ can be obtain from the stopping time of a free particle in LAMMPS.
$\Delta$ can be obtain from the (mean) difference of velocities before and after the collision.
$\alpha$ I don't know since the collision is velocity dependent with LAMMPS... (Is is that much important?)

Still, even for this model, it doesn't looks like the equilibration time depends on $N$:


### Critical point of the exponential model

Fraction of active particles **at the transition**:

Forget the "during 2 snapshots t=1s" below, I forgot to remove it, it doesn't make sense. Energy and velocity are instantaneous quantities
_______

It's kind of hard to say if the transition is of first or second order:

With more time to equilibrate, it seems that it's a first order phase transition:

Since we can't escape the absorbing state and we are not at $N\rightarrow \infty$ I think the exact phase diagram is hard to obtain.

In any case, I don't think it describes the second order phase transition observed numerically (and eventually experimentally). Both the gif are taken at the transition, the first one is for the exponential model, the second of is for the constant kick model:

The model with a $\Delta$ seems closer:

### Other model
With the model presented in the article found by Andrea (where we add a kick $\propto\Delta$ constant at each collisions). At the first sight, it looks like we obtain a second order phase transition but in fact, the curve doesn't change with $N$ and the growth is not really that of a power law, with a sharp transition at a critical $\phi$:


I also checked, even with a gargantuan time of simulation, the curve doesn't change.

I tried with different $\gamma$ and $\Delta$ but we still obtain the same shape :(

#### Finite size effect in the exponential model?


Not much....

### Model:

$\alpha = \alpha_0e^{-v/v_o}$ the coefficient of restitution where $v$ is the difference of velocity between the two particles colliding and $\alpha_0>1$
$\gamma$ the drag.

The model contains 3 free parameters.

We observe a fine phase transition, whatever the treshold used is:


The point of using the MSD is that we observe a confinement when $\phi\rightarrow\phi_{max}$

|Treshold MSD = $r/5$ in an interval of 1sec| Treshold MSD = $r/10$ in an interval of 1sec |
|:---:|:---:|
|||

---

### Some rough phase diagram

To plot the phase diagram we don't use any treshold, we only plot $D$ or $E$, for constant $v_0 ( =100\gamma r_b)$ and $\gamma$:


Some values:


----

For constant $\alpha_0 (= 1.5)$ and $\gamma$:


Some values:

----

For constant $\alpha_0 (= 1.5)$ and $v_0 (=100\gamma r_b)$:


Some values:

_____

I tried to obtain the divergence of the relaxation time, it's highly dependent on the treshold:
Treshold MSD = $r/10$ in an interval of 1sec:

Treshold MSD = $r/5$ in an interval of 1sec:

Treshold MSD = $r/12.5$ in an interval of 1sec

Treshold MSD = $r/20$ in an interval of 1sec (we can't really trust the results, but the general shape seems alright):

It doesn't seem that increasing $N$ increases the height of the peak... It's a bit troublesome. Perhaps, it's treshold dependent (in the sense that we could observe a "true" divergence only when we have the good parameter to determine wether the particle is dead or alive).

_____

### What to take as a parameter to determine whether a particle is dead or alive

#### $D$ might be good...


Reaching an absorbing state:
Reaching a SS:

But it doesn't lead to a peak for the time to reach the SS (it gives something akin to the TRESHOLD MSD $r_b/5$).

#### What about $E$ ?

The energy is a great parameter to determine whether a particle is dead or alive. But it's harder to use it in order to get a timescale.
If we start from a very high energy, we have a nice exponential but the accuracy is not high enough to extract a time scale since the decrease is very abrupt.
If we start from a low energy, **below the transition and for large $\alpha_0$**, it does this funky thing:


{%youtube 2vuEChJD0lY %}

The shape seems pretty robust, we always obtain this curve (even the precise values) for a large number of particles. And no timescale can easily be extracted from that.

**Above the transition and for tiny $\alpha_0$ (which corresponds more to the experimental case)**, it looks better (before the transition, the convergence to the absobing state must be driven by large fluctuations or by this local blob that we can see of the video that lead(s) to the linear decrease of $E$. It's a bit weird still since, if the particles are freeflying: $dv/dt = -\gamma v$, the energy should decrease as an exponential...):


But as with the MSD, it's $N$ independent.

As for now, I didn't find a parameter that made $\tau$ increase to infinity in the thermodynamic limit.

### Velocity at the transition
For system dilute enough, at the transition (the for $\alpha_0 \simeq 1.5$), the velocity distribution is close to a gaussian.

But, when the transition happens at a low density; we obtain this kind of thing:

Funny remark: with the phase diagram at $v_o$ and $\gamma$ fixed, the value of $D$ (or $E$) was increasing with lower values of $\phi_c$. In this case, $D$ at $\phi_c$ has decreased, compared to what we add for the first phase diagram. We expect that $D(\phi_c)$ increases till a certain point and then decreases. IF THE TRANSITION IS OF FIRST ORDER EVEYWHERE. Perhaps the system didn't thermalize since it takes a lot of time at low density...

Thinking about it, the absorbing phase transition is really sensitive on $N$ since, as long as we are not in the thermodynamic limit, the system will reach the absorbing phase in a finite time (somewhat exponential in $N$ I think, so it might be alright). Perhaps it disrupts a bit my predictions.

### Hyperuniformity at the transition

WRONG, DONT MIND, BUT I WANT TO NOW WHY
I plotted $\frac{\langle N^2\rangle - \langle N\rangle^2}{r^2}$ (where the average is on initial points $r_0$ where we calculate $N(r)$, with $r$ the distance between the particles and $r_0$) to see whether the fluctuations were suppressed:


It doesn't look particularly hyperuniform:

I also obtain the same curve with random point, so I guess I'm missing something.

# Recap of the 10/05's meeting

- Need to check a bit more if the absorbing state thing shows hyperuniformity (and I've been highly enlighted by the fact that the order parameter is basically the diffusion coefficient by the way...)

- The fact that it looks like we can't self assemble S1 at high dissipation is interesting. Might push forward in this direction (do something that I quite didn't get with the MSD now that I think of it...). We didn't talk about it but I think that the lesser $\alpha$ is, the higher we need $\phi$ to be. Maybe at some point we can't have a system dense enough to self assemble.

- Maybe check the evolution of MSD with time

- The jump in energy was caused by the fact that I sometimes took the snapshot of the system after the kick (high $E$) and sometimes before (low $E$).

---

## Background
Here you find the folder "GranularIntro" where some reference on granular materials. In the README 4 some guidelines for the reading.
https://www.dropbox.com/scl/fo/pfwt910len8rvafj9a27d/h?dl=0&rlkey=w754ok6bwjdt16qv8xkngf0nm

## Plans

## To Do

- Start simulations for an S1 configuration with $\alpha=1$ (keep sytem size relatively small)
- Explore $\alpha <1$

### [25/04/2022]
#### Simulation parameters and notes:
* $N = N_b + N_s$ (N big + N small)
* $r_s, r_b, m_b, m_s$ radius and mass
* $L$ size of the square box
* $\phi$ packing fraction
* $\gamma$ the friction coefficient with the plate.
* $T_b$ the temperature of the thermostat/kick. The kick is idealized by $\eta$ the gaussian noise with variance $\langle\eta(t)\eta(0) \rangle = 2\gamma T_b\delta(t)$. It should make sense to do this as long as the frequency of kicks is >> than the frequency of collisions.
* $\Delta t_n$ the time between each kicks/update of the velocity (should also be << time between collisions if we really want a langevin noise. It's an issue for performance at high density. Perhaps it's not necessary to simulate the kick as a "perfect" brownian noise).
* $\alpha$ the restitution coefficient (between 0.7 and 1 for the simulation to proceed smoothly with a reasonnable noise, if $\alpha$ is to small the particles start to clump together. It drastically increases the number of collisions/real time and creates overlap (inelastic collapse?))
* **MODEL 1:** $$v_i\rightarrow v_i - \frac{\gamma v_i}{m_i}\Delta t_n + \frac{\sqrt{2\gamma T_b\Delta t_n}}{m_i}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$
* **MODEL 2:**$$v_i\rightarrow v_i - \gamma \Delta t_n v_i + \sqrt{\frac{2\gamma T_b\Delta t_n}{m_i}}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$
#### The algorithm
For now I can only generate a square lattice of small and big particles.
Then: usual EDMD with noise, at each $N\Delta t_n$ we add a kick to all the particles or a subset of them. In the subset case, with $\alpha = 1$, the granular temperature doesn't converge to $T_b$ but the FDT should still holds with a different effective thermostat temperature. Indeed, $T_b$ is not anymore the temperature of the thermostat felt by the system. I was expecting a scaling of this effective thermostat temperature with the square root of the fraction of particles we update but no... See below

The update of the velocities each $\Delta t_n$ is done as follows:

$$v_i\rightarrow v_i - \gamma \Delta t_n v_i + \sqrt{\frac{2\gamma T_b\Delta t_n}{m_i}}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$

This scaling with the mass **ensures the equipartion** at $\alpha = 1$ and is said to reproduce experimental measurements https://arxiv.org/pdf/cond-mat/0501488.pdf (p.18 binary mixture) and https://www.researchgate.net/publication/10962384_Driven_low_density_granular_mixtures where two different models are given. We'll see soon that it fails to predict the repartition of the kinetic energy given by LAMMPS.
With this scaling, $\gamma\sim t^{-1}$ should be the characteristic time (no mass appears as in the standard langevin equation), the SS might be reached at the same speed whatever the mass is (does this imply something cool?). "Explosion" might appears (of course I wasn't able to reproduce it....) if $m_s/m_b\sim0.001$ but I guess it's not too harmful for our orders of magniturde.

The scaling with $\Delta t_n$ ensures that the physics doesn't (at least I hope...) change when we change the time between kicks (as long as it's smaller than the time between collisions of course).

The term $\gamma$ also prevents a collective drift. Here is the trajectory of some particles when $\gamma=0$ and a non zero noise.

The experimental $\gamma$ is really tiny. It's not clear that we need it. We could either discard it and put walls or we could conserve total momentum by kicking particles by pairs (and giving them opposite momenta).
Note added after the simulations are done: the characteristic time with the LAMMPS simulations seems extremely large. So $\gamma$ should be way smaller than what I tested.

:::success
- Test it for 2D hard hardspheres case and constant coefficient of restitution
:::

In what follows: $T_b = 1$
#### FDT


Pretty noisy but looks good since $T_b = 1$ (here we update **all velocities** each 0.1 time). The characteristic time is good too.

----------
Now, let's kick one particle out of 100 each $\Delta t_n= 0.01$. We obtain a very weird coupling between $T_g$, $T_b = 1$ and $\gamma$ (scaling choice of $\Delta t_n$?). The SS still varies as $1/\gamma$ however... No apparent FDT (we might define an effective temperature). Huge time scale.

----------

Same thing but we pick N/10 random particles at each kick.
Converges to $1.5\neq T_b$ with high $\gamma$.

----------

Same thing but with $N$ random particles at each kicking time (some are kicked twice and some none). We get back the FDT.

#### Energy at $\alpha \neq 1$
In what follows; at each kick we update all particles.

Even when $\alpha\neq 1$ is constant, $T_g$ at the SS depends on $\gamma$ (expected since we are not at equilibrium).

#### Maxwell Boltzmann?


#### Radial distribution function: Equilibrium vs out of equilibrium

Low packing fraction


______
High packing fraction

____

Binary mixture at N_b = N_s (only showing $g(r)$ of the big particles)

#### Binary mixture and equipartition


This is what we obtain for the model: $$v_i\rightarrow v_i - \gamma \Delta t_n v_i + \sqrt{\frac{2\gamma T_b\Delta t_n}{m_i}}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$

_____


and this is what we get with $$v_i\rightarrow v_i - \frac{\gamma v_i}{m_i}\Delta t_n + \frac{\sqrt{2\gamma T_b\Delta t_n}}{m_i}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$

_____


Opposite of what we observe in LAMMPS. In LAMMPS the ratio goes to $\sim 1$ with the packing fraction for a low enough noise!

Allow to discriminate models? Depends on packing fraction too... So Packing fraction + Ratio + $\alpha$ might give $\gamma$ ???

### [27/04/2022]
#### α>1
Test with: $\alpha>1$, $\gamma \sim 0.1$ and $T_b = 0$. Either absorbing state or divergent energy->no SS:

For example, for this packing fraction: $\alpha = 1.57$:
{%youtube uBZ6z9FOz5A %}
The explosion of velocity is due to fluctuation of density.
The energy is drained to 0 when $\alpha$ is below 1.57 for this packing fraction.

#### First S1 tests:
With $\alpha < 0.9$, $0.001<\gamma<0.1$, $\phi = 0.79$, $N_b=N_s$, $r_s/r_b = 0.42$, $m_s/m_b = 0.06$ and additive disk we observe a phase separation **only with the first scaling where we divide by $m$ the noise and the viscous coeff** (model 1)


From what I saw, the model 2 (with no scaling with $m$ for the friction and a $1/\sqrt{m}$ factor for the noise) doesn't lead to a demixing. But it gives $E_b/E_s>1$.

Starting from: 

|MODEL 1 | MODEL 2 |
|:---:|:---:|
|||

________

### [29/04/2022]

#### S1 formation without noise at $\alpha = 1$
Additive disks took too much time to self assemble. This is the best simulation with $\phi= 0.82$ and $r_s/r_b = 0.41$ for a ~20 hours run:
|Start | End |
|:---:|:---:|
|||

So from now on, I will run simulations with non additive disks.

_______

With non additive disk, the tiny particles can easily sneak under the big ones and we see S1 forming after a couple of minutes (phew...):

------------------------------

#### S1 formation with noise AT EQUILIBRIUM ($\alpha = 1$): $\phi = 0.85$, $T_b = 1$ and $r_s/r_b= 0.457$

~1h30 of simulation:
##### MODEL 1

|$\gamma = 0.003$ | $\gamma = 0.006$ |$\gamma = 0.01$ |
|:---:|:---:|:---:|
||||
##### MODEL 2

|$\gamma = 0.003$ | $\gamma = 0.006$ |$\gamma = 0.01$ |
|:---:|:---:|:---:|
||||

MODEL 2 might be better than MODEL 1 (hard to say on 1 simulation). $\gamma$ low might be better too (which also implies a tiny noise).

#### S1 formation with noise OUT OF EQUILIBRIUM ($\alpha \neq 1$): $\phi = 0.85$, $T_b = 1$ and $r_s/r_b= 0.457$

It doesn't look like we experience demixing with MODEL 1 in this case (maybe dense enough).

For now, it seems that any $\alpha<1$ tends to destroy the order. For example, starting from a S1 crystal obtained by a simulation with $\alpha = 1$ and starting from it with the same parameters except $\alpha = 0.95$ we get:

|0 collision | $10^8$ collisions |
|:---:|:---:|
|||

The same thing happens for MODEL 1 and MODEL 2.

We also have a similar effect even without $\gamma$:
$$v_i\rightarrow v_i + f(m_i)\sqrt{2 T_b\Delta t_n}\zeta~~~~~~~\zeta\sim\mathcal{N}(0,1)$$
with $f(m_i)\sim m_i^{\beta}$ and $\beta = -1/2$ or $-1$ or $0$. Perhaps we need a wall?

### [02/05/2022]

#### $q_4$ with noise OUT OF EQUILIBRIUM ($\alpha \neq 1$): $\phi = 0.85$, $T_b = 1$ and $r_s/r_b= 0.457$
Last line is $\alpha = 1$
**MODEL 1:**


**MODEL 2:**

#### Different $\alpha$ for large particles and tiny ones

For example $\alpha_{bb} = 0.2$ and $\alpha_{bs} = \alpha_{ss} = 1$

________________

Starting from a crystal with MODEL 2, all $\alpha=1$ except for one:


Looks quasi stable when the only dissipative collisions are the ones between small particles. It's also the only case where it looks like we recover the equipartion (pretty unexpected) even though the system is "strongly" out of equilibrium since $E/(NT_b)\neq 1$ (but we are not that far from one compared to the others simulations below)! The closer we are from equilibrium, the more stable we seem to be.

The drop of energy at 60000s corresponds to a loss of the global 4-fold symetry!




Inversion of $E_b/E_s$, perhaps it means that to we need a different restitution coefficient for small and big particles (in order to reproduce what we see with LAMMPS).

#### Comparison of MODEL 1 and MODEL 2


Equipartition is broken for MODEL 2 when the 4-fold symetry is lost (at $t>400000$s)!!!

### [04/05/2022]

We observe S1 formation at high packing fraction/high size ratio (and high $T_b$ even though it doesn't look like it's that much important, simulations to test this are running).


It's not crystal clear with $q_4$ but it is with ovito ($\phi = 0.88$ and $r_s/r_b = 0.46$):
|0 collision|~10^9 collisions |
|:---:|:---:|
|||

However the structures don't look stable ($\phi = 0.86$ and $r_s/r_b = 0.48$) perhaps due to the "low" packing fraction. S1 might be able to form at higher packing fraction and be stable, but the crystallisation might take a lot of time:

|t = 300000s|t = 350000s |
|:---:|:---:|
|||

_________________

The temperature doesn't affect the convergence to S1. Nonetheless, simulations with a small temperature are faster (less collisions). So effectively, it's best to keep a small temperature:


________________________

### [05/05/2022]

With the right packing fraction we obtain almost a perfect S1 (the time scale is also almost surely wrong ($t_{max}\sim 450000s$) but it's not that much important, curse my lazyness and my crappy python script).

$\phi = 0.8716$ and $r_b/r_s = 0.4767$:

_______

Until now, all simulations were done at $dt_N=0.01$ (because I thought we needed a time between each kick much smaller than the time between each collision):


We can increase it and still get good results.
_____

We might get something akin to a quasicrystal, for example (it's really really not clear):

$\phi = 0.8$, $r_s/r_b= 0.55$, $N_b/N = 0.3$

Not the exact frame we saw on the previous figure, but the same run:

________________

$\phi = 0.87$, $r_s/r_b= 0.5$, $N_b/N = 0.3$


### [06/05/2022]

For 10k particles at s1:

For 5k particles:

________

The value of $\alpha$ is also important (S1 parameters and $T= 10$):


For $\alpha$ too low we don't observe any S1 formation (perhaps because in this case we need a higher density? Or perhaps because the energy of the system is not high enough to get out of local extrema (of whatever quantity governs the phase transition)) and for an $\alpha$ too large the simulation takes forever because it leads to an increase of the number of collisions. Without friction it's not too troublesome because the temperature doesn't play any role in the dynamic and more collisions == "faster time" == crystalisation happens in a shorter *real* time. But the situation looks different out of equilibrium (and perhaps at $\alpha = 1$ with $\gamma \neq 0$ but it's a bit strange to me...).

It's a bit annoying because it means that we might have to fine tune the parameters for EACH phase we want in order to obtain the formation of the wanted crystal. Perhaps not.... Nonetheless, it means that if we want to explore an even more damped system (which would mean farer from equilibrium? I'm not so sure...) we will have to find the new appropriate parameters.

### [09/05/2022]
Some configurations and their k-space:
https://www.dropbox.com/scl/fo/y4cozri23i1trmyfl8a81/h?dl=0&rlkey=e8bqtuqxmdcd2g2mmlsi9fcui

A lot of them might not have equilibrated now that I look at some of them...

The S1 doesn't seem that hard to form. Perhaps, the dissipation just "shifts" to a higher density (or pressure? Is it quantifiable?) the non dissipative phase diagram : 

We also observe, I think, QC12 (but it's definitively not clear... At least to me):

A weird one too:

There is also this, I'm not sure what this is (the real space looks like the usual Hex):
and perhaps this ?:


Also a clear QC8 I think:

______

I tried to observe the S1 crystallisation with a tiny $\alpha$ but I can't make it auto-assemble, whatever the $T$ or the $dt_N$. The $dt_N$ tested might be too high in the case of a tiny $\alpha$.

### [11/05/22]

I'm trying some things with the MSD and walls. In the meantime here are all my QC12 structure factors: https://www.dropbox.com/scl/fo/xzkurngafeewk4uebr63k/h?dl=0&rlkey=xtdyl7o93osnh02oj7tsf7waz. I also wanted to see if the initial condition and rare events were really important and I got this for S1:


Which is a bit troublesome...

_____
In my case, the temperature is lower at the walls (BPC on the $y$ axis and reflective wall with a coeff of rest = $\alpha = 0.95$ on the $x$ axis).

And in the case of binary mixtures
(S1 parameters:)

To simulate the shaking of the walls, perhaps, I could transform them in a thermostat (adding noise to the velocity of the particles each time they kick a wall).

### [12/05/22]

Configurations forming S1 in the dissipative case are completely frozen in the non dissipative case! That's really funny to watch, nothing moves, even in 100k seconds..

____

Walls don't prevent the formation of S1 (the islands of S1 form at the center, same thing that what we observed with LAMMPS from what I understood).
$\phi = 0.868$ and $r_s/r_b = 0.457$


As we saw earlier, since the evolution of the configuration is highly dependent on the set of random numbers drawed during its simulation, it's hard to reach any conclusion or to compare the formation of S1 with or without walls, but it seems that it doesn't change a thing. At least for large systems.

____

It doesn't look like we can observe QC8 with $\alpha = 0.85$ (which surprises me a bit...) even though, with the naked eye, the configurations look similar to the usual QC8 (maybe not enough strips):


Here are all the configurations: https://www.dropbox.com/scl/fo/s6dy2fqtnuv02mt7ts8nj/h?dl=0&rlkey=2y81sd2l82yrjydt2nb5drao3

_____
At first I plotted with a tiny $dt$:

$$D = \dfrac{\langle||\vec r(t+dt) - \vec r(t) - \langle \vec r(t + dt) - \vec r(t)\rangle ||^2\rangle}{dt}$$

But it looks like I'm still in the ballistic regime of D:

and

It probably means that the mean free time is greater when the system forms a crystal.
____

### [13/05/22]
#### Diffusion and how does the crystal form?

I plotted the unwrapped coordinates of the particles after some time for a coeff of restitution equal to 1 (first figure) and equal to 0.95 (second ond). In both case we have a thermostat. I tried to have an average energy equal for both simulations but it's not really the case....


Comparing the two cases:

in loglog scale:
When we separate big and small particles:

And the energy:

It looks like that we can see the transition to a crystal in the MSD (at $2\times10^5s$) for the dissipative case

We see that for the thermal case the particles giggled way more than in the dissipative case. The energies are different, but even when the energy of the thermal simulation is lower than the dissipative simulation its MSD can be higher (An other simulation with equal energy is running to show it!)

Perhaps when we decrease $\alpha$ we kill the coefficient of diffusion (We observe a stronger collectif drift too which "drains" the energy given by the kicks into something useless for the self-assembly) and the particles don't have enough velocity to rearange in a solid. But it's unlikely since a S1 crystal (formed at $\alpha = 0.95$) melts with $\alpha = 0.85$.

Also, even in the thermal case, we see the increase of $E$ with $q_4$.

_____
#### Very vague ideas

In the thermal case all configurations are equiprobable as long as they are allowed by geometric arguments. In the dissipative case; configurations where particles are close to each others should (I think) be more probable than the others.

In this case:
- Can we found back the thermal case if we simulate the dissipative case with high energy? Intuitively, I don't know if increasing the temperature tends to equalize all probabilities **in this case**.
- If this is just a question of temperature, we should find back S1 even when $\alpha \sim 0.85$ which isn't the case (for the range of temperature I simulated).
- Can we link this fact to the necessary increase of pressure to form S1 in the dissipative case (what is really the meaning of this?).
- Can we find a pressure treshold at which S1 won't form that is independent of $\alpha$ (if we can even define a meaningful pressure in the dissipative case). It's probably not the case since S1 won't even form with $\alpha = 0.85$ at high temperature and the pressure should not be that small in these simulations...
-We can run simulation for Hex and see for different $\alpha$ how the pressure at the phase transition changes. But in the Hex case, the crystal will always forms --whatever $\alpha$ is-- as long as we are at a high enough packing fraction. It seems different for S1. For Hex, the effect might even be the opposite, we would observe a phase transition at a lower packing fraction than in the thermal case since the particles are more clamped in the dissipative case

### [16/05/22]

The MSD in the dissipative case presents a weird behaviour on large time scales (https://www.dropbox.com/scl/fo/dl3jom5g7w8mklw69svwh/h?dl=0&rlkey=fzk13x1cld1zyi0im32hunkub):

It might be corrolated to the crystal growth, but the duration doesn't look right everytime.

On top of that, the thermal system has indeed way more mobility than the dissipative one (the energy are close enought to compare):

and

I also managed to observe some very unclear QC8's with $\alpha = 0.85$:

All updated configurations are here: https://www.dropbox.com/scl/fo/s6dy2fqtnuv02mt7ts8nj/h?dl=0&rlkey=2y81sd2l82yrjydt2nb5drao3

### [17/05/22]

We observe S1 at $\alpha = 0.8$ for very high densities and long time (x20 compared to the thermal case. I don't know which part of it is explained by the dissipation, and which part is explained by the extremely high density):

### [02/06/22]

I was playing a bit with the pressure to see if I observed (and where?) a phase transition even in the dissipative case for a monodisperse system of hard sphere:


The main issues with the 2 last plots is that each point corresponds to a different temperature, I need to change my thermostat over time in order to equilibrate at the same temperature for each configuration.

For example, at each ~$\gamma^{-1}$ (time to equilibrate) I could change $T_b$ this way:

$$T_b\rightarrow = T_be^{E_{wanted} - E}$$

Until $T_b\simeq E_{wanted}$ and then measure the pressure.
There are a lot of way to do this, it's perhaps not the best one. If I do this, I get for the dissipative case:


The phase diagram is better. Nonetheless it seems hard to converge to $E_{wanted} (=1)$ especially at high $\phi$ (some tests indicate that the system simply didn't have enough time to equilibrate at high $\phi$).

If we compare the dissipative case with the thermal one:

It looks like that the transition happens at a higher packing fraction for the dissipative case but it's hard to tell... The thermal pressure is almost always higher, it's probably due to the fact that the energy of the dissipative system is a bit lower than the energy of the thermal system.

### [03/06/22]

The transition is less clear with a dissipative system, but we might expect to find it again at higher packing fraction.


As suggested by Andrea, the pressure might not be a good parameter to analyze a phase transition out of equilibrium. I will try to play a bit with $q_6$ next time.

### Phase diagram

Equilibrium:

### Funny things while waiting for the cluster

#### Cascade of energy


(Extremely) coarse grained velocity field
{%youtube CWTh8lpioII %}

As in the navier stokes equation we observe a cascade of energy toward the low k modes. A bit unexpected because the nature of the dissipation is not the same.

* [What about the breakdown of equipartition?](https://physics.stackexchange.com/questions/20944/the-equipartition-theorem-in-momentum-space)
* Can we add energy "at some scale" to see the eventual forward enstrophy cascade (still in $E(k)$). The constant $E(k)$ for large $|k|$ might be caused by the fact that we are free cooling the gas. What happens if we add a little bit of energy?
* Is it not just a lucky artefact of the coarse graining (hard to define a velocity field without an aggressive coarse graining since we have clusters and big holes)?
* What happens in a binary mixture?

____________

- Explore Models where the coefficient of restitution depends on the velocity of the incoming particles.
## References

### OPEN QUESTIONS

* Do we see crystallisation and *quasi-cristallization* in this *simplified* model ?
*