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

https://arxiv.org/pdf/1804.03638.pdf Phase separation through time dependent delta
https://arxiv.org/pdf/1801.00850.pdf Discontinuous to continuous liquid to solid via h/A change

### A quick recap of the results (with or without direct coexistence method)

#### Self assembly of the langevin model:

Here, dissipation shift toward higher densities the crystallization

#### Self assembly of the delta model without damping

Herre, the delta (at moderate $\alpha$) shifts the crystallization toward lower values. Note however, that at fixed $\alpha$ and changing $\Delta$ the crystallization stays the same (it is not a question of renormalization of velocity and time)!

#### Rationalization:


The $\Delta$ model has an effective repulsion while the langevin has an effective attraction!

#### Attractive well at equilibrium

Indeed a small attraction shifts the crystallization towartd higher values of the packing fraction:
But a too strong one will completely make the system crystallize:

### Direct Coexistence Method!

#### Direct Coexistence langevin noise!

Here, we perfom some direct coexistence simulation with an expanding crystal (we start from a crystal in vacuum that can expand... and then eventually create a crystal/liquid):

We change $\alpha$. We also note that for small $\alpha$ the interface is, often not as flat (slab like) than for an equilibrium system:

We do not, in a first time, care about the fact that the solid might have an additional stress due to the y confinement. We perform the simulation with different initial condition so that we can always be in the coexistence. This will also change the y spacing for the same state point and will allow us to check if the results are really that much dependent on the y spacing (in the end we should really take care of that neately...). Everything is at fixed forcing $T$ thus, decreasing $\alpha$s have a decreasing energy.

Here is a very preliminary result:

$\phi_i$ is the initial density, squares are the liquids and circles are the solids.

We see that $\phi_i = 0.844$ is in the coexistence region until more or less $\alpha = 0.95$ at which point, we lose the solid phase and every points below are around $0.844$.

$\phi_i = 0.8557$ is more or less always in the coexistence region but for some reason, the solid phase is losing density. I think this is because my algorithm computing the density of each phase is fucked...

$\phi_i = 0.867143$ is pretty good for the dissipative region but is messed up at high $\alpha$ because i think it had not had the time to equilibrate.

Anyway, the lattice spacing does not seem to mess up thing too much:

### Coexistence with $\Delta$ alone

As expected the coexistence does not change qwith fixed $\alpha$. Note also the two points for the liquid and solid for the same density. These are different latticve spacing.


(the delta axis has to be read in reverse)

(the delta axis has to be read in reverse)

Here is a first test with full $\alpha$, this is very very noisy and the liquid density for $\phi = 0.82$ is fucked up because we are a bit too close to the solid phase:
($\alpha = 0.283$ at $\phi = 0.82$)

while the solidof $\phi = 0.81$ has a strange decrease, which migh be in part explained by a lack of thermalization?

### Coexistence with $\Delta$ + damping
fixed $\alpha$ I'm varying $\Delta$.
#### $\alpha$ fixed


I think that for reasonably high $\alpha$, with $\Delta$ alone we do not change the crystallization zone with varying $\Delta$ because $\Delta$ without $\gamma$ simply renormalize the time. BUT, with $\gamma$, decreasing $\Delta$ corresponds to getting closer to the absorbing region and hence having really strong repulsion, which might explain why the system crystallize at a lower density with decreasing $\Delta$

### Temperature difference between phases

Langevin model:


Delta:

### Clean data.

Now that we have a nice idea about the behavior of the system, we'd like to have clean data.
The crystal must not be constrained in the y direction. In our simple example, we can simply try to find the case in which the x and y lattice spacing are the same. This would also corresponds to the case where the initial crystal has the same density as the crystal coexistence density.

This is a bit painful to do because it requires us, for each new simulation (new $\alpha$, $\gamma$, $\Delta$), to look at a various number of different lattice spacing and determine the one that lead to the equilibrated $x$ lattice spacing.

From a simulation:

We compute the bonds between neighbors


Which then gives us the lattice spacing on average, here for pure $\Delta$ at $\phi = 0.84$:
We then select the y lattice for which it coincides with the x lattice spacing

Unfortunately, this procedure has to be redone for every change of parameters because it will lead to a different crystal density, for example, langevin $\alpha = 0.95$ and $T = 0.5$ (the data are a bit messy because it has not had the time to thermalize).


This one requires a way denser system, look at the fluid phase emerging from the melting of a very dense crystal!!!


____________

For now, let's be lazy and use some old data:

langevin:

Increasing temperature might tend to make the coexistence smaller, as in Equilibrium.


These data were taken without care and we did not care about the x and y spacing (everrysimulation has the same spacing 2.085 with particle diameter = 2), we can nonetheless, look a posteriori if the data might make sense by looking at the hx/hy ratio.

The solid lines represent the ratio between the measured y lattice and the one imposed at the beginning. If the simulation went well, we should expect this quantity to alway be equal to 1. Here, for example, we see that for small $\alpha$ we had some issues. And indeed when looking at the snapshots, we easily had a bubble forming instead of a slab, and even some times, the whole crystal rotating. For higher $\alpha$ it is fine.

For the others point, we can focuse more on the ratio between the x and y lattice spacing (which are the points). Here is a zoomed part:


We see that they always wander around 1, not perfect perfect, but perhaps, one could say, good enough to be sure that the trend observed on the phase coexistence diagram is correct.

If we look at the distribution of lattice spacing, they are really close together:

____________



Turns out, the phase diagram obtained from this method is worse :)


___________

### NewData!



____________





__________


#### Reaching the eq limit:

$T_{thermo}$ is the temperature of the langevin bath while $E$ is the measured energy. The vertical lines are the equilibrium results:

https://arxiv.org/pdf/1003.4462.pdf




| Theory (just the general shape, not the right values) | results |
| -------- | -------- |
|  |  |

____

However, the system is still strongly in a non-equilibrium state because the particles have a different temperature:


Result by trizac:


indeed, the result is expected to not vary too much with the temperature of the thermostat (because it would only change the pair correlation function!). At least for a liquid.

______

#### Test with monodisperse
Note that this information of Equilibrium -> $T_{thermo} = 0$ might not be encoded in the g of r.

|Res varying in colormap| |  |
| -------- | -------- | -------- |
|E or E/T varying colormap|||

With a finer grid:


### Theory of Nonequilibrium Symmetry-Breaking Coexistence

https://arxiv.org/pdf/2309.10341.pdf

At equilibrium:

$$\mu_{s} = \mu_l\Rightarrow \int_{v_l}^{v_s}(p(v) -p_{coex})dv$$

$$p_s = p_l = p_{coex}$$

$$T_s = T_l$$

Out of equilibrium, there is no reason to have an equality of temperature but it will require egality of pressure. Moreover, if $\mu$ is not defined, the coexistence really requires that the number of particles going from $s$ to $l$ are equal to the number of particles doing the opposite thing.

For a system with dissipative collisions, the pressure will depend on the density and on the phase $\psi$. It will also depends on the temperature, itself dependent on the density (and probably the phase):

$$p(\phi, \psi, T)\equiv p(\phi, \psi = \{l, s\}, T(\phi, \psi))$$

Thus it might make sense to take simply these two equations and call it a day, for the coexistence of dissipative hard disks:

$$\int_{v_l}^{v_s}(p(v, T(v)) -p_{coex})dv$$

$$p_s = p_l$$

where it is understood that $T(v)$ also depends on the branch/phase we chose for the pressure.

The issue with this method is that it requires to know the point at which we switch from the solid to the liquid equation of state. Ie: the spinodal.

______________

This is too hard :) Let's do liquid gas:

With a varying $\Delta$:

$$\Delta(t_i, t_j) = \Delta_{0}(2-e^{t_j/\tau_{r}}-e^{t_i/\tau_{r}})^\beta$$

Still using the usual collision rule:

$$\Delta E_{coll}(t_i, t_j) = m\Delta(t_i, t_j)^2+\alpha\Delta(t_i, t_j)\sqrt{\pi m T}-T(1-\alpha^2)$$

Let's assume poissonian collision: $$P(t_i) = \dfrac{1}{\tau}e^{-t_i/\tau}$$
with $\tau$ the inverse collision rate: $\tau = \omega(\phi, T)^{-1} = (8\phi\chi\sqrt{T/\pi m}/\sigma)^{-1}$.

Easier, but we can do simpler:

$$\Delta E_{coll} = \Delta E_i(t_i) + \Delta E_j(t_f)-T(1-\alpha^2)$$

$\Delta E_i=\delta E + \Delta E_0(1-e^{-t_i/\tau_r})^\beta$ or $\Delta E_0\times\text{min}((t/\tau_r)^\beta, 1)$, turns out the exponential one is easier to deal with.

The steady state is given by:

$$\langle\Delta E_{coll}\rangle = 2\langle \Delta E_i(t_i)\rangle -T(1-\alpha^2) = 2\delta E + 2\Delta E\dfrac{\tau/\tau_r\Gamma(1+\beta)\Gamma(\tau/\tau_r)}{\Gamma(1+\beta + \tau/\tau_r)}-T(1-\alpha^2) = 0$$

and shows a VdW loop in the pressure:


________________________

# MIPS

### Intro:

Underdamped mips not well studied.

Here is the model I'll be using:

$$m\dfrac{d^2\boldsymbol x_i}{dt^2}=\gamma \left(v_0 \boldsymbol n_i-\dfrac{d\boldsymbol x_i}{dt}\right)-\dfrac{d U_{WCA}}{dx_i}+\sqrt{2\gamma T}\eta_i(t)$$
with $\boldsymbol n_i = (\cos(\theta_i), \sin(\theta_i))$ and:

$$\dfrac{d \theta_i}{dt}=\sqrt{2D_r}\zeta_i$$

We chose U steep enough and change it with $v_0$ so that we are effectively in the hard sphere limit.

We will mostly place ourselves in the limit $T = 0$.

In this case, the system is controled by 3 parameters: $\tau_m = m/\gamma$ some inertial time, $\tau_r = 1/D_r$ an rotational time/the time needed to change direction and finally $l_0=v_0\tau_r$ the distance a particle travels in a straight line in overdamped limit.
_________
### Quick overview of overdamped mips:
3D:
$l_0=v_o/\gamma*m$

#### Theory:

From hydrodynamics:





We obtain:
Far away from the boundary the gradient termS will vanish leaving equyality of pressure (purely mechanically).



NOTE THAT USUALLY, THE MODIFIED MAXWELL CRITERIUM WILL DEPEND ON THE PREFACTOR OF THE GRADIENT TERMS! IT DEPENDS ON THE INTERFACIAL ENERGY

_________
$Pe\sim v_o/D_r$

Here https://www2.thphy.uni-duesseldorf.de/~hlowen/doc/op/op0416.pdf it says that MIPS is lost when increasing Pe. Which is not what I find. Here https://iopscience.iop.org/article/10.1088/1367-2630/abd80a they argue the same ($f_o = \gamma v_o$):

They also argue a lot of stupid things I think...

Here is what I find:



Nontheless, the los of MIPS with increasing mass is genuine. At least probably genuine :) Since both coexistence branches converge to thge same point. It must be genuine :) :)

| small simulation time | long simulation time |
| -------- | -------- |
|  | |

_____________

#### Why do we lose MIPS at high m?

| Small m/gamma vs rotational time | Big m/gamma vs rotational time |
| -------- | -------- |
| |  |
| | |

The right video is just here to illustrate.

Here is a video where at the middle, we change the mass to a larger one:

#### Why is MIPS highly dependent on $N$ when we increase $v_o$?

##### Is it interfacial size?
Since we increase the velocity, perhaps, the size of the interface grows uncontrollably with $v_0$ because the particles, after hitting the dense phase, a gas particle might jump far away (not expected since it is damped)

Not really!

The interface size grow with m tho!
##### Is it dynamically impossible to reach mips at high $v_0$ and small $N$ but it is stable..

Perhaps, here is a phase diagram starting from a phase separated system. Everything is well phase separated!:

|$\langle \phi\rangle$| Starting phase separated | Starting homogeneous |
|----| -------- | -------- |
|$\langle\phi\rangle = 0.6$ t= large|  |  |
|t=medium|||
|t=small|||

Here is the size of the box divided by the persistence lenght: 

|$\langle \phi\rangle$| Starting phase separated | Starting homogeneous |
|----| -------- | -------- |
|$\phi = 0.7$|||

The nuclei is harder to form!

Might not be the end of the story! Here is a starting config at small vo and increasing vo suddendly:



### Toward an effective $v(\rho)$

We define $v(\phi) = \langle \boldsymbol n \cdot \boldsymbol v\rangle$ the effective velocity.

Here are the things we obtain 

the colormap is the value of $v(\phi)/v_0$

______________________________________________

Increasing m decreases the $v(\phi = 0)$. And the decay is not linear anymore!


$v(\phi = 0)/v_o$ is independent on $v_o$

$v(\phi)/v_o$, even for small mass is dependent on vo

| Column 1 | Column 2 |
| -------------------------------------------------- | -------------------------------------------------- |
|  |  |

$v(\phi = 0, m)/v_o$ decreases almost linearly with $m$ (probably like $1/(1+m)$).

The thing is that, while increasing $m$ we might decreases the $v(\phi)$ observed but it might not matter too much as long as the reorientation time is extremely large.

_________

The diffusion coefficient:

(dots are measured and lines are extrapolated from the ideal theory).

____________

### Bidisperse:
$\langle\phi\rangle = 0.65$


_______________________
Start homogeneous:


BOTTOM LEFT POINT



__________________

Not well equiliçbrated:
| homo start | Mips start |
| -------- | -------- |
|  |  |

#### Theory of free inertial AOUP

$$m \mathbf{\ddot x}=-\gamma \mathbf{\dot x} + \gamma v_0 \mathbf{n} + \mathbf F+ \sqrt{2\gamma T}\mathbf{\eta} ~~~\textrm{ and }~~~\dot \theta = \sqrt{2D_r}\zeta$$

Or equivalently ($2D_r = 1/\tau_r$ and $\gamma/m= 1/\tau_m$ ):

$$\mathbf{\ddot x}=- \mathbf{\dot x}/\tau_m + v_0/\tau_m \mathbf{n} + \mathbf F/m +\dfrac{\sqrt{2\gamma T}}{m}\mathbf{\eta} ~~~\textrm{ and }~~~\dot \theta = \sqrt{1/\tau_r}\zeta$$

Or **probably** equivalently (with $F=0$) with a change of variable $t = \tau_r\tau$:

$$\mathbf{\ddot x}=- \mathbf{\dot x}\dfrac{\tau_r}{\tau_m} + \dfrac{\tau_r^2}{\tau_m} v_0 \mathbf{n} + \tau_r^2\dfrac{\sqrt{2\gamma T/\tau_r}}{m}\mathbf{\eta} ~~~\textrm{ and }~~~\dot \theta = \zeta$$

From which (only the first two :)) we obtain:

#### Explicit formula:

$$\mathbf{v}(t)=\mathbf{v}_0e^{-t/\tau_m}+\dfrac{1}{m}\int_0^t\left[\mathbf F(t')+\gamma v_0 \mathbf n(t')+\sqrt{2\gamma T}\eta\right]e^{(t'-t)/\tau_m}dt'$$
$$\mathbf{x}(t)=\mathbf{x}_0+\tau_m\mathbf v_0\left(1-e^{-t/\tau_m}\right)+\frac{1}{\gamma}\int_0^t\left[\mathbf F(t')+\gamma v_0 \mathbf n(t')+\sqrt{2\gamma T}\eta\right]\left(1-e^{(t'-t)/\tau_m}\right)dt'$$
$$\theta(t) = \theta_0 + \sqrt{1/\tau_r}\int_0^t\zeta(t')dt'$$

#### Correlation:

\begin{split}\langle (\theta(t_1)-\theta(t_2))^2 \rangle&=\dfrac{1}{\tau_r}\left\langle\left( \int_0^{t_1}\zeta(t')dt'- \int_0^{t_2}\zeta(t^")dt^"\right)^2\right\rangle\\
&=\dfrac{1}{\tau_r}\left(-2\int_0^{t_1}\int_0^{t_2}\delta(t'-t^")dt^" dt' + \int_0^{t_2}\int_0^{t_2}\delta(t'-t^")dt^" dt' + \int_0^{t_1}\int_0^{t_1}\delta(t'-t^")dt^" dt'\right)\\
&=\dfrac{1}{\tau_r}\left(-2\int_0^{t_1}\Theta(t_2-t^")dt^" + \int_0^{t_2}\Theta(t_2-t^")dt^" + \int_0^{t_1}\Theta(t_1-t^")dt^"\right)\\
&=\dfrac{1}{\tau_r}\left(t_2+t_1-2\min(t_2, t_1)\right)=\dfrac{1}{\tau_r}|t_2-t_1|
\end{split}

\begin{split}
\langle \mathbf{n}(t_1)\cdot\mathbf{n}(t_2)\rangle&=\langle \cos(\theta(t_1)-\theta(t_2))\rangle=\langle\Re(e^{i\Delta \theta})\rangle = e^{-\langle \Delta \theta^2\rangle/2}=e^{-\dfrac{|t_1-t_2|}{2\tau_r}}=e^{-D_r|t_1-t_2|}
\end{split}

\begin{split}
\langle \mathbf{v}(t_1)\cdot\mathbf{n}(t_2)\rangle &= \dfrac{v_0}{\tau_m}\int_{-\infty}^{t_1}\langle \mathbf{n}(t')\cdot\mathbf{n}(t_2)\rangle e^{(t'-t_1)/\tau_m}dt'=\dfrac{v_0}{\tau_m}\int_{-\infty}^{t_1}e^{-D_r|t'-t_2|-(t_1-t')|/\tau_m}dt'\\
&=\dfrac{v_0}{1+\tau_m/\tilde\tau_r}e^{-|t_1-t_2|/\tilde\tau_r}
\end{split}

assuming $t_2>t_1$ and defining $\tilde \tau_r=1/D_r$:

\begin{split}
\langle \mathbf{v}(t_1)\cdot\mathbf{v}(t_2)\rangle&=\dfrac{1}{m^2}\int_{-\infty}^{t_1}\int_{-\infty}^{t_2}\left[(\gamma v_0)^2 \langle\mathbf n(t')\cdot\mathbf n(t) \rangle +2\gamma T\langle\eta(t')\cdot\eta(t)\rangle\right]e^{(t'+ t-t_1-t_2)/\tau_m}dtdt'\\
&=\dfrac{2 T}{m} e^{-(t_2-t_1)/\tau_m}+\dfrac{\gamma^2 v_o^2}{\gamma^2-(mD_r)^2} \left(e^{-D_r(t_2-t_1)}-\dfrac{mD_r}{\gamma}e^{-(t_2-t_1)/\tau_m}\right)\\
&=\dfrac{2 T}{m} e^{-(t_2-t_1)/\tau_m}+\dfrac{v_o^2}{1-(\tau_m/\tilde\tau_r)^2} \left(e^{-(t_2-t_1)/\tilde\tau_r}-\dfrac{\tau_m}{\tilde \tau_r}e^{-(t_2-t_1)/\tau_m}\right)
\end{split}

\begin{split}
\langle \mathbf{x}(t_1)\cdot\mathbf{x}(t_2)\rangle&=\int_{-\infty}^{t_1}\int_{-\infty}^{t_2}\langle \mathbf{v}(t)\cdot\mathbf{v}(t')\rangle dt dt'=2(t_2-t_1)\int_0^{t_2-t_1}\langle \mathbf{v}(t)\cdot\mathbf{v}(0)\rangle dt-2\int_0^{t_2-t_1}t\langle \mathbf{v}(0)\cdot\mathbf{v}(t)\rangle dt\\
&=2\left(\dfrac{2\tau_m^2T}{m}-v_o^2\dfrac{\tilde\tau_r\tau_m}{(\tilde\tau_r/\tau_m)^2-1}\right)\left(t/\tau_m+e^{-t/\tau_m}-1\right)+2v_o^2\dfrac{\tilde\tau_r^2}{1-(\tau_m/\tilde\tau_r)^2}\left(t/\tilde\tau_r+e^{-t/\tilde\tau_r}-1\right)\\
(t\gg\tau_m, \tilde\tau_r)&\simeq 2\left(2\dfrac{T}{\gamma}+v_0^2\tilde\tau_r\right)t
\end{split}