Spatio-temporal correlations in Coulomb clusters

Dynamics of particles in a complex medium is interesting both from fundamental [1,2] and technological [3–5] viewpoints. For example, heterogeneous dynamics [6,7] in glasses show an exponential tail [8] in the distribution of the particle displacements, and belong to the emerging class of "Fickian yet non-Gaussian" motion. A more exotic decay with a "stretched Gaussian" tail ($\sim \exp\left[-(x/x_0)^{\beta}\right]$, $1< \beta < 2$) has been recently found for nanoparticles confined in periodic arrays of nanoposts [2,9]. While the role of hydrodynamic drag [1] and the "core" size of interactions [10] on such slow relaxation are issues of active debate, the local heterogeneity, inherent in many systems, is also a prime suspect. Exploring these inhomogeneities is central to our study.

One way to induce a heterogeneous response is by confining the system. Motion of trapped particles changes non-trivially not only near the boundary [11], but also deep inside the system, in particular, if the inter-particle interactions are long-ranged. Rich interplay of interactions and the trap geometry on the dynamics of particles at finite temperature $(T)$, have been observed in the presence of long-range interactions like for dipolar colloids [12], dusty plasma [13] and Yukawa liquids [14]. With increasing experimental control, clusters of few electrons interacting via Coulomb forces, are fabricated having sizes comparable with the range of interactions [15]. While random impurities describe disorder in bulk systems, a finite system experiences disorder from the irregularities on its "soft" boundary created by electrostatic and magnetic means [16,17]. Such conjecture on the origin of disorder finds support from the measurement of mesoscopic fluctuations [18–20] in chaotic dots. The solid phase of Coulomb clusters at low T is called Wigner Molecule (WM), because it mimics the Wigner crystal [21]. A WM in a circular or irregular trap (defined as CWM and IWM, respectively) melts into a liquid upon increasing temperature [22–25]. The question we ask is: How do a finite number of Coulomb particles reflect the signature of the irregularity in the confining potential and also help in characterizing the thermal melting? Here, we address this question by probing the subtleties of dynamical responses in terms of the Van Hove correlation function (VHCF) [11] for two-dimensional (2D) IWM and CWM, using Molecular-Dynamics (MD) [26] simulations. The VHCF has been routinely used in unfolding the intricacies of spatio-temporal correlations at different length and time scales, yielding crucial information on static and dynamic properties [8,27].

We show that the probability distributions of the displacement of particles, given by the self-part of the VHCF, exhibit a pronounced departure from the standard Gaussian behavior [11]. At low T, this shows a nearly exponential behavior at intermediate times, whereas its long-time profile features multiple peaks. At higher T, the differences between the two traps become more pronounced: While the distribution of the displacement of particles in CWM shows stretched Gaussian relaxation near melting, the IWM shows a stretched exponential trend up to a rather large T. Such distinct dynamical responses are consistent with the mean square displacement of the particles, the spatial displacement pattern and the distinct part of the VHCF. We extend the phenomenological model proposed to explain the non-Gaussian behavior [28] to understand the subtleties of the complex motion in our system.

We model a WM using the Hamiltonian: ${\cal H} = \phi \sum_{i<j} |\vec{r}_i - \vec{r}_j |^{-1} + \sum_i V_{\text{conf}}(\vec{r}_i)$. The first term represents Coulomb repulsion between particles of charge q in the medium with dielectric constant , and $\phi=q^2/\epsilon$. We consider long-range (unscreened) Coulomb repulsion because the interaction is poorly screened in finite clusters with a small number of particles. Here, $V_{\text{conf}}$ is the confinement potential. We compare results from two traps: a) irregular trap, $V_{\text{conf}}^{\text{Ir}}(\vec{r})=a\{x^4/b+by^4-2\lambda {x^2}{y^2} +\gamma (x-y)xyr\}$ [25,29], and b) circular trap, $V_{\text{conf}} ^{\text{Cr}}(\vec{r})= \alpha r^2$, where $\alpha = m \omega_0^2/2$. We rescale the length $(r \rightarrow \phi^{1/3}\alpha^{-1/3} r)$ and energy $(E \rightarrow \phi^{2/3}\alpha ^{1/3} E)$, that transform [22] the CWM Hamiltonian to ${\cal H}^ {\text{Cr}}=\sum_i r_i^2 + \sum_{i<j} |\vec{r}_i-\vec{r}_j|^{-1}$. The time scale is renormalized by $t\rightarrow \hbar~\phi^{-2/3}\alpha^{-1/3} t$ (see footnote ¹ ). The overall factors a and α are chosen to keep the mean density in the two traps same for a fixed total number of particles, N = 150. Our model of irregular confinement generates dynamics ranging from periodic to chaotic single-particle motion by tuning λ [29]. Such motion along with broken spatial symmetries are taken as footprints of disorder. The irregularity parameters $b=\pi/4$, $\lambda \in [0.565, 0.635]$ (which controls chaoticity), and $\gamma \in [0.10,0.20]$ (which breaks reflection symmetry) [31] are adjusted to access the standard features of disordered systems [32]. Distances are measured in the unit of r₀, the mean inter-particle spacing between neighboring particles at the lowest T. The statistics of results from $V_{\text{conf}}^{\text{Ir}}$ are improved by "disorder averaging" [32] over ten configurations, each identified by a specific $\{\lambda,\gamma\}$. Results from CWM are averaged on ten independent MD simulations.

Our calculations are set up in NVT-ensemble using velocity rescaling [26]. The circular confinement conserves the angular momentum to zero (rigid rotation is subtracted from CWM results), and IWM naturally violates it. After equilibration, we study various dynamical quantities, as follows, in order to recognize the intricate response of our system. The density correlation function (obtained from VHCF) [11] is defined as $G(r,t)= \langle \sum^{N}_{i,j=1}\delta \left[ r - |\vec{r}_i(t) - \vec{r}_j(0)|\right] \rangle$, where i and j are particle indices. Here $\langle.\rangle$ denotes the average over time origins. The self-part $(i=j)$ of the VHCF, $G_s(r,t)$, gives the probability that a particle has traversed a distance r in a time interval t. Its distinct part $(i \neq j)$ of VHCF, $G_d(r,t)$, represents the probability of finding at t a different particle at a distance r from an initial particle. The mean-square displacement (MSD) is defined as ${\langle \Delta r^2 (t) \rangle = N^{-1} \sum_i \langle\left[ \Delta \vec{r}_i(t) \right]^2\rangle}$, where $\Delta \vec{r}_i(t) = \{\vec{r}_i(t) - \vec{r}_i(0)\}$.

A thermal crossover from a "solid" to "liquid" in WM asserts [22,25] a "solid" at $T=0.006$ and "liquid" at $T=0.03$ for the chosen parameters. Here we consider the temperature range $T=0.002\text{--}0.25$ which encompasses the equilibrium crossover. In order to study the spatio-temporal correlations in our system, we first investigate the time dependence of $G_s(r,t)$ [27], represented in a semi-logarithmic plot in fig. 1. Thus, at any given time t, a Gaussian r-dependence of $G_s(r,t)$ appears as an inverted parabola in this plot. Depicting the low $T~(=0.006)$ results of $G_s(r,t)$, fig. 1(a) illustrates that the dynamics in IWM can be broadly classified into three temporal regimes² : i) Ballistic motion occurs for short time $(t \leq 1)$ yielding Gaussian $G_s(r,t)$. Particles feel each other's presence only for t > 1 implying that our unit time represents the time scale of interaction. ii) Particles get "arrested" in the cage formed by neighbors during the intermediate times $(100 \leq t \leq 1000)$, reflected in the bunched up traces of $G_s(r,t)$ in the inset of fig. 1(a). iii) The dynamics for large time $(t \geq 1000)$ develops multiple peaks in $G_s(r,t)$. The motion in the CWM (fig. 1(b)) is broadly similar to IWM insofar as the temporal regimes are concerned, albeit the caging regime is narrow $(10 \leq t \leq 100)$ and untenable, causing weaker agglomeration of $G_s(r,t)$ (inset, fig. 1(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At high temperatures, in a normal liquid, particles undergo diffusive motion resulting in a Gaussian $G_s(r,t)$. We find that for both $V_{\text{conf}}^{\text{Ir}}$ and $V_{\text{conf}}^{\text{Cr}}$, $G_s(r,t)$ deviates from Gaussian even beyond the thermal crossover point $(T=0.02)$. The deviation is more pronounced in $V_{\text{conf}}^{\text{Ir}}$, where it undergoes a progressively slower decline. Such differences in the behavior of $G_s(r,t)$ persist in the two confinements (fig. 1(c), (d)) at temperatures beyond "melting".

Next, we quantify the r-dependence of $G_s(r,t)$, for a given t, by fitting it to two functional dependences: $G_{s}^{\text{small}}(r,t) \sim e^{-{r^2}/c}$ for $r \leq r_c$, and $G_{s}^{\text{large}}(r,t) \sim e^{-l r^k}$ for r > r_c. This is implemented in the parameter sub-space of $\{t,T\}$, in which the r-dependence of $G_s(r,t)$ has a monotonic evolution (this turns out to be $t \sim 400$ for $T \le 0.006$ and the entire t range for T > 0.006). For each region (r smaller or larger than r_c), our fitting procedure yields a value for $\chi^2$, measuring the goodness of the fit in that region. The final accepted value of r_c is the one that minimizes the total $\chi^2$ of the fit on the entire r-axis. The same prescription also obtains the other optimal fit parameters, e.g., c, l, k, etc. The inset (a1) of fig. 2(a) illustrates the procedure described above.

Figure 2((a) CWM; (b) IWM) presents the time evolution of the parameter k that quantifies the nature of the decay of $G_s(r,t)$. Gaussian behavior corresponds to k = 2. Henceforth, we denote the corresponding k for irregular and circular traps as $k_{\text{Ir}}$ and $k_{\text{Cr}}$, respectively. For low $T~(\leq 0.006)$, both $k_{\text{Cr}}$ and $k_{\text{Ir}}$ degrade close to unity (insets (a2) and (b1)) establishing a near-exponential decay [8] in the "solid" for large times. With increasing T, the long-time limit of $k_{\text{Cr}}$ rises; however, it broadly remains within $1 < k_{\text{Cr}}< 2$ (main panel of fig. 2(a)) confirming a stretched Gaussian trend, as observed in recent experiments [2,9]. In striking contrast, $k_{\text{Ir}}$ decreases with t causing a stretched exponential (i.e., $k_{\text{Ir}} < 1$) spatial decay (fig. 2(b)). Interestingly, $k_{\text{Ir}}$ reaches its minimum for $T=0.02$, a temperature that marks the crossover between "solid" and "liquid" [25,33]. Beyond this crossover temperature, the traces of $k_{\text{Ir}}(t)$ start rising after an initial decrease up to $t \sim 10^2$. The irregular system begins to resemble a conventional liquid (i.e., $k_{\text{Ir}} \approx 2$ for large t) only for $T \sim 0.25$. In contrast, the circular confinement yields $k_{\text{Cr}}\rightarrow 2$ at large t for $T \geq 0.05$ only. We emphasize that our result of k = 2 implies that the whole $G_s(r,t)$ curve has a single Gaussian fall. Thus, the nature of VHCF identifies and distinguishes a solid-liquid crossover in different trap geometries. While a multiple-Gaussian [34] fall (all centred at r = 0) with different widths in different ranges of r can also have k = 2, our results do not correspond to such behavior.

The classification of temporal regimes at low T is also realized in MSD $(\langle\Delta r^2 (t) \rangle)$ as shown in fig. 3((a) IWM; (b) CWM) in a double-logarithmic plot. At all temperatures, for both the confinements, we see for $t<1, \langle\Delta r^2 (t) \rangle \sim t^2$. Thus, it represents the ballistic regime. After this initial temporal regime, the low-T evolution of MSD is surprisingly flat for IWM from strong caging, whereas, CWM shows a weak Brownian behavior. Such a wide plateau in MSD (spread over several decades) for an IWM, as well as its occasional jumps, are reminiscent of glassy dynamics [35]. The plateau at intermediate t shrinks gradually with T causing a smooth crossover between ballistic and Fickian motion. The "jumps" also fade away with T. In contrast, neither a clear plateau, nor such jumps are discernible in CWM even at low T. At long time, MSD saturates when displacements approach the system size. MSD is also estimated from the VHCF as $\langle \Delta r^2 (t)\rangle = \int r^2 G_s(r,t)\text{d}^2r$, and its evaluation is shown as dots in fig. 3(a), (b), confirming the Fickian regime for $T \geq 0.01$. Interestingly, these regimes are rife with non-Gaussian motion, cf. fig. 3(a), (b). Thus, our results corroborate the notion of "Fickian yet non-Gaussian" dynamics [1,36] at long times. Overall signatures of MSD indicate that IWM resembles a supercooled system [35] more than CWM does.

Zoom In Zoom Out Reset image size

Does the spatial reorganization of the moving particles gel with the above unfolding of $G_s(r,t)$? In order to develop a deeper understanding, we analyze the particle displacements, $\Delta \vec{r}_i(t)$ in fig. 4. In the figure, thick dots represent the initial position of the particles and the connecting line signifies the displacement of the same particle at the specified time t. Different colors represent particles with different ranges of magnitude of the displacement, $\Delta r(t) = |\Delta \vec{r}(t)|$. In fig. 4, we see that most particles only rattle around their equilibrium positions for $t \ll 1000$. Motion that begins by breaking a cage causes avalanches of such events correlated in space, leading to a string-like path of delocalization [37,38]. Avalanche of such cage breaking events, whenever occurring, sets up in a very short time interval. Figure 4(a) presents $\Delta \vec{r}_i(t=1000)$ for an IWM at $T=0.006$, showing string-like alignment of neighboring particles yielding a slow relaxation producing a long tail in $G_s(r,t)$. Once the collective motion subsumes most of the system, particles undergo relatively small displacements, but an overall drift, producing secondary peaks. With further increase in t, the string-like path degenerates into sub-clusters, each performing independent collective motion (fig. 4(b)), and contribute to the multi-peak $G_s(r,t)$ of fig. 1.

Zoom In Zoom Out Reset image size

The differences between the two traps stand out in the nature of motional coherence. The strings at low T in a CWM (fig. 4(c), (d)) become approximately circular respecting azimuthal symmetry. This symmetry is well maintained at intermediate and long times, particularly near the boundary. Intriguingly, the mobile particles travel by several "lattice spacings" riding on the coherence even in the "solid" phase³! Signatures of such spatially correlated motion of particles, have been realized in thermal melting of circular plasma crystals [13], and also in a flocking transition of driven granular classical matter in confinements [39].

A phenomenological proposition to access the non-Gaussian behavior of $G_s(r,t)$ is to express it as a convolution of Gaussian particle displacements and their distribution of diffusivities P(D) [28,40]. We generalize this idea to accommodate its multi-peak nature at low T and large t:

$$\begin{equation} G_s(r,t) = \sum_{p=1}^{N_p} \int\text{d}D_p \frac{P(D_p)}{(4\pi D_p t)}{\exp}\left[ -\frac{(r-r_p) ^2}{4D_p t}\right], \end{equation} \tag{ 1 }$$

assuming $G_s(r,t)$ has N_p peaks. Naturally, a monotonic decay results for p = 1. A departure from the Gaussian behavior is thus crucially attributed by P(D) which is now a superposition of the individual distributions $P(D_p)$ for each p. We access $P(D_p)$ by first obtaining r_p and the corresponding peak width by fitting the individual peaks of $G_s(r,t)$ (fig. 5(a)). We then determine $P(D_p)$ from the statistics of individual particles whose $\langle {r_i}^2(t)\rangle$ lies within the peak width at a given r_p (see footnote ⁴). The resulting P(D) shows a long tail indicating striking inhomogeneities at all different temperatures. Such broad distribution is a footprint of dynamic heterogeneity [35], i.e., multiple time scale relaxations in a single system. P(D) shows exponential, and stretched exponential tails for up to intermediate times at low $T~(\leq 0.006)$ and becomes unusually broad with near-flat sections at large t. The long tail (mostly stretched exponential) in P(D) persists even at high $T~(\geq 0.030)$. The different nature of P(D) in separate ranges of D at long t and low T, is consistent with the physical picture presented in fig. 3. Interestingly, the $G_s(r,t)$ obtained from eq. (1) using the aforementioned extraction of P(D) matches almost perfectly (insets of fig. 5(b)–(d)) with direct simulations. While in fig. 5(b)–(d) data are shown only for irregular confinement, a broad matching holds even for circular confinement.

Zoom In Zoom Out Reset image size

The multiple relaxations are also probed by $G_d(r,t)$, the distinct part of VHCF. In the small-time limit $G_d(r,t\rightarrow 0)$ represents pair correlation function, g(r) [11,25], and thus consists of sharp peaks, like crystals, at low T. Figure 6(a) shows depleted peaks of $G_d(r,t\rightarrow 0)$ for IWM indicating a lack of positional order even at low $T~(=0.006)$ [25], while a CWM displays sharper peaks (fig. 6(b)) due to the azimuthal ordering in the symmetry direction. It is interesting to note the distinct role of the "fast" and "slow" particles towards the time evolution of $G_d(r,t)$. A particle is qualified as "fast" ("slow") if its mobility lies within the top (bottom) 7% (10%) – the mobility of the i-th particle being defined as $\langle |\Delta \vec{r}_i(t)| \rangle$. The long-time accumulation of $G_d(r \rightarrow 0)$ (fig. 6(a), (b)) arises from the fast particles occupying the locations from where another particle moved out – thereby developing the spatially correlated motion. In contrast, the slow particles trace out the sharp peaks of $G_d(r,t)$. The inset of fig. 6(b) illustrates such distinct behavior for $V_{\text{conf}}^{\text{Cr}}$. A similar picture emerges from the self-part, too – while the slow particles produce the small-r (Gaussian) behavior, the fast particles contribute to the long tail, see the inset of fig. 5(a). Figure 6(c), (d) shows $G_d(r,t)$ at $T=0.030$ for the IWM and CWM, respectively. As expected in a liquid, the large-r peaks are already washed out at smallest t. Its relaxation with time smoothens out even the weak modulations leaving only the strong peak at r = 0. The long-time filling of the gap at $G_d(r \rightarrow 0, t\rightarrow 0)$ seems stronger in the irregular trap than in circular confinement at $T=0.03$, unlike for low $T~(=0.006)$ where such filling is similar, irrespectively of the geometry. This is indicative of the fact that the spatially correlated motion causing such filling has more pronounced effects on the irregular trap, which, in turn, leads to ultra-slow relaxation in them.

Zoom In Zoom Out Reset image size

In conclusion, we show that the interplay of interactions and trap geometry leads to dynamic heterogeneity, producing unconventional dynamical behavior, like an ultra-slow spatial relaxation in Wigner molecules. We ensure that our qualitative results survive for ${35 \le N \le 500}$, and we expect these to persist for other disordered traps [15]. VHCF and its variants are directly probed in experiments [34]. While the exponential tail in $G_s(r,t)$ is common in experiments and simulations on colloids [36], gels [34], glasses [7,8,41,42], and bio-molecule suspensions [43], our finding of its stretched exponential decay opens a new direction, to the best of our knowledge. Our results will motivate further experiments probing the role of local heterogeneity on complex motion providing guidance in designing nanotechnology devices exploiting the geometry of the confinement.

We thank Deepak Dhar and Chandan Dasgupta for valuable discussions. BA acknowledges his doctoral fellowship from the University Grant Commission (UGC), India.