Pattern formation with repulsive soft-core interactions: Discrete particle dynamics and Dean-Kawasaki equation

Brownian particles interacting via repulsive soft-core potentials can spontaneously aggregate, despite repelling each other, and form periodic crystals of particle clusters. We study this phenomenon in low-dimensional situations (one and two dimensions) at two levels of description: by performing numerical simulations of the discrete particle dynamics and by linear and nonlinear analysis of the corresponding Dean-Kawasaki equation for the macroscopic particle density. Restricting to low dimensions and neglecting fluctuation effects, we gain analytical insight into the mechanisms of the instability leading to clustering which turn out to be the interplay among diffusion, the intracluster forces, and the forces between neighboring clusters. We show that the deterministic part of the Dean-Kawasaki equation provides a good description of the particle dynamics, including width and shape of the clusters and over a wide range of parameters, and analyze with weakly nonlinear techniques the nature of the pattern-forming bifurcation in one and two dimensions. Finally, we briefly discuss the case of attractive forces.

Figure 1

Dynamics of a system of 6000 particles evolving according to Eq. (3) in a periodic one-dimensional domain of size $L=3$, so that ${\rho }_0=2000$. The interaction potential is GEM-3 with parameters $R=0.1$ and $\epsilon =0.0333$. (Left) $D=0.4$ so that $\tilde{D}=D/\left(\epsilon {\rho }_{0}R\right)=0.06$. (Right) $D=0.6$ so that $\tilde{D}=0.09$. The lower panels show the spatiotemporal trajectories (in black). For clarity of the plots, only 600 trajectories are shown here. The top panels show a coarse-grained normalized density, $〈\rho \left(x\right)〉/{\rho }_0$, at the late stages. The coarse graining is done by averaging spatially over boxes of width $0.02R$ and temporally over the last 150 temporal configurations separated by ${10}^{-3}$ time units. In the top left panel, the clusters are fitted by Gaussians (black lines).

Figure 2

(Left column) Snapshots of the positions of $N=1000$ particles in 2D at large times for $R=0.1,\epsilon =0.0333,L=1$, and ${\rho }_0=1000$. (Right column) Coarse-grained densities of the same configurations. The coarse graining is done by averaging in space over distances $0.017R$ and in time over 500 configurations separated by ${10}^{-3}$ time units. (a),(b) GEM-3 potential with $D=0.02$ so that $\tilde{D}=D/\left(\epsilon {\rho }_0{R}^2\right)=0.06$. (c),(d) GEM-3 potential with $D=0.01 \left(\tilde{D}=0.03\right)$. (e),(f) GEM-1 potential with $D=0.005 \left(\tilde{D}=0.015\right)$. Figure (g) is the density $\rho$ along the white dotted lines shown on panels (b) and (d) (red and blue squares, respectively). The black curves correspond to a fit by a sum of Gaussian functions.

Figure 3

(a) Radial distribution function $g^{\left(2\right)}\left(r\right)$ and (b) structure factor $S\left(k\right)$ [as a function of the dimensionless wave vector $k=qR$, where $q$ is the Fourier variable in the definition of $S\left(q\right)$ in Eq. (6)] for several values of $\tilde{D}=D/\left(\epsilon {\rho }_0{R}^2\right)$ in 2D particle simulations for $L=3,\epsilon =0.0333$ and $R=0.1$. GEM-3 potential and $\tilde{D}=0.05$ (solid green circles), 0.06 (red squares), 0.07 (blue diamonds), 0.08 (yellow triangles), 0.09 (black stars) and GEM-1 potential with $\tilde{D}=0.05$ (open gray circles). (c) Normalized maximum value, $S_{max}/N$, of the main peak of the structure factor plotted with respect to $D$ (inset) or $\tilde{D}$ (main plot) in 2D particle simulations for $L=9,\epsilon =0.0333$, and $R=0.1$. Yellow triangles, ${\rho }_0=2000$; green squares, ${\rho }_0=1500$; red diamonds, ${\rho }_0=1000$; solid black circles, ${\rho }_0=500$. The open blue circles correspond to results obtained by integration of the 2D deterministic DK equation for the same parameters and ${\rho }_0=2000$. In that case, $S\left(k\right)$ was obtained by taking the Fourier transform of the density-density correlation function.

Figure 4

Steady solutions of Eq. (12) in 1D for $\alpha =3$. (a) $\tilde{D}=0.06$; (b) $\tilde{D}=0.09$.

Figure 5

Density functions in the steady state for the GEM-3 potential obtained by integration of Eq. (12) with the pseudospectral method. (a) $\tilde{D}=0.06$; (b) $\tilde{D}=0.03$. Panel (c) is the density $\rho$ along the white dotted lines shown in panels (a) and (b) (red and blue squares, respectively). The black curves correspond to a fit by a sum of Gaussian functions.

Figure 6

Growth rates, from Eq. (13) (a) 1D. GEM-3 potential with $\tilde{D}=0.075$ (solid green squares), 0.084 (red triangles), 0.09 (blue diamonds), and 0.1005 (yellow circles) and GEM-1 potential with $\tilde{D}=0.075$ (open black squares). (b) 2D. GEM-3 potential with $\tilde{D}=0.015$ (solid green squares), 0.03 (red triangles), 0.06 (blue diamonds), and 0.12 (yellow circles) and GEM-1 potential with $\tilde{D}=0.015$ (open black squares).

Figure 7

Effective potential $h\left(x\right)$. Solid line, GEM-3 potential, with intercluster distance $c=a/R=1.38$, which is (see Table 1) the periodicity given by the linear stability analysis. Dashed line, $h\left(x\right)$ for GEM-1. We use the same value $c=1.38$, although the peak at the origin is independent of the intercluster separation. Dotted line, $h\left(x\right)$ for the GEM-3 potential but for a larger intercluster separation $c=1.8$.

Figure 8

Clusters width as a function of $\tilde{D}$ for a GEM-3 potential: red circles, from the steady density solution of the DK equation; blue squares, from a coarse-grained density in particle simulations. In both cases, the Gaussian fit was made focusing on the top of the clusters. (a) One-dimensional case. Dashed black line is Eq. (27). (b) Two-dimensional case. Dashed black line is Eq. (33). The intercluster distance $c$ used in the analytical formula was obtained from the position of the second peak of the radial distribution function $g^{\left(2\right)}\left(r\right)$ [see Fig. 3], which gives a result slightly better than using the linear instability value $2\pi /k_c$. In the particle case the coarse graining in 1D and in 2D was done as in Figs. 1 and 2, respectively.

Figure 9

For each value of $\tilde{D}$ the maximum and the minimum value of the steady density in 1D for $\alpha =3$ are plotted. The inset gives an enlargement of the bifurcation region. Squares $(\square )$, values from direct numerical simulation of the DK equation; *, $\times$, and +, values from a coarse-grained density in particle simulations at ${\rho }_0=1000$, 1500, and 2000, respectively ($R=0.1,\epsilon =0.0333)$. Solid line, the weakly nonlinear approximation Eq. (34), including up to the second harmonic; dotted line, only the first harmonic term in Eq. (34); dashed line, the maximum of the pattern as given by the Gaussian approximation, Eq. (30). In the particle case the coarse graining was done as in Fig. (1).

Figure 10

Maximum value ${\rho }_{max}/{\rho }_0$ of the steady density solution of the DK equation in 2D for $\alpha =3$. We start from $\tilde{D}=0.077$, slowly increase it up to $\tilde{D}=0.102$ (red disks), and then slowly decrease $\tilde{D}$ back to its initial value (black squares and green diamonds). The plot clearly highlights a hysteretical behavior. The green diamonds indicate that ${\rho }_{max}/{\rho }_0$ remains close to 1 during the times accessible to our simulation, but clearly increasing, even if very slowly. For the black squares in the lower branch, ${\rho }_{max}/{\rho }_0$ is observed to decrease towards 1. Thus, the behavior is consistent with the theoretical threshold for stability of the homogeneous solution ($\tilde{D}\approx 0.0823$; see Table 1, marked with a dotted vertical line) and we expect the green diamonds to reach the upper density branch in a sufficiently long simulation. The thin undotted lines give the theoretical prediction for ${\rho }_{max}$ from Eq. (36): two branches, giving the upper one a prediction for the amplitude of the stable hexagonal density. The right dotted vertical line indicates its turning point.

Figure 11

Dynamics of 600 Brownian particles with attractive interactions in a periodic domain of size $L=3$, so that ${\rho }_0=200$. $\epsilon =-0.33,R=0.1$, and $D=3.96\times {10}^{-3}$, so that $\tilde{D}=6\times {10}^{-4}$. (a) GEM-1 attractive potential; (b) GEM-3 attractive potential.

Figure 12

Density functions at large times for the attractive GEM-3 potential and $\tilde{D}=0.0975$. (a) Coarse-grained density of particle simulations. $R=0.1,\epsilon =0.0333,N=2000,L=1$ (so that ${\rho }_0=2000$). The coarse-graining procedure is as in Fig. 2. (b) Integration of Eq. (12) with the pseudospectral method. Panel (c) is the density $\rho /{\rho }_0$ along the white dashed lines shown in panels (a) and (b) (red and blue squares, respectively), shifted for convenience. The black curves correspond to a fit by a sum of Gaussian functions.

Figure 13

Growth rate [Eq. (13)] for attractive GEM-3 and GEM-1 potentials (solid and open symbols, respectively). For GEM-3: $\tilde{D}=0.015$ (green squares), 0.9 (red circles), and 3 (yellow diamonds). For GEM-1: $\tilde{D}=0.6$ (blue squares), 4.5 (black circles), and 21 (orange diamonds).