Collective motion and its connection to the energy landscape in 2D soft crystals

Abstract

We examine the collective motion in computational models of a two-dimensional dusty plasma crystal and a charged colloidal suspension as they approach their respective melting transitions. To unambiguously identify rearrangement events in the crystal, we map the trajectory of configurations from an equilibrium molecular dynamics simulation to the corresponding sequence of configurations of local potential energy minima (“inherent structures”). This inherent structure (IS) trajectory eliminates the ambiguity that arises from localized vibrational motion. We find that the evolution of the IS trajectory in the crystal can be split into comparatively longer-lived ground states and shorter-lived discrete excited states. These discrete excited energy levels are a consequence of discrete numbers of defect clusters in the crystal. We find that the collective rearrangement occurs through different mechanisms: (i) small closed-loop motion in the ground states without the facilitation of defects, and (ii) much larger and complex open-ended particle motions in excited states that are facilitated by clusters of defects. In both cases, clusters of displacing particles can be separated into much smaller groups of replacing particles with a loop-like structure. In contrast to glass-forming liquids, the mass of the rearranging groups grows on heating towards the melting temperature rather than cooling. We find that crystal melting in these systems can be anticipated by the merging of the average time the crystal spends in the ground state with the average time in the excited states.

---

1 Introduction

Understanding the collective dynamics of particle rearrangements in both amorphous and crystalline materials, including supercooled liquids, gels, and colloidal systems, has been an active research topic for many years.^1–9 It is observed that the constituent atoms, molecules, or colloids with the greatest displacements over a specified time window are spatially correlated, often in the form of particle exchange motions termed “strings” because of their polymeric structure.^10,11 However, the continuum range of collective particle displacements on a hierarchy of timescales often makes it challenging to unambiguously identify specific groups of collectively rearranging particles. This difficulty has led to various criteria for identifying clusters of cooperatively rearranging particles.^7,10,12–15 Typically, a fraction of particles with the greatest displacement over a specific time window are examined, where the displacement results from the accumulation of multiple small displacements that are difficult to uniquely identify. Nonetheless, a consistent picture has emerged regarding how the average size of the rearrangements varies with temperature or density/packing fraction. To improve the understanding and quantification of this collective motion and its implications, it is valuable to examine a system where rearrangement events can be unambiguously identified in both space and time.

To this end, we investigate model crystalline materials in two dimensions as minimal systems exhibiting collective particle exchange motion; we consider two variations that represent either a dusty plasma crystal (DPC) or a charged colloidal crystal^16–20 which exhibits similar dynamic features. These crystals offer several simplifying features that make them ideal model systems: crystals have well-defined lattice positions that provide an underlying template for stable locations and two-dimensional materials allow for convenient visualization and quantification of dynamic heterogeneity in both space and time. We can further simplify the identification of rearranging groups by mapping the snapshots of a trajectory of particles onto the so-called “inherent structure” (IS) trajectory through energy minimization of each snapshot.²¹ The combination of these factors allows us to directly relate the collective motion to discrete states in the potential energy landscape (PEL).

The spatial and temporal heterogeneity in particle displacements of quasi-2D DPC and colloidal crystals have been widely examined,^6,22–24 and it is known that string-like collective particle exchange motion is prevalent.^23,25,26 Large-scale particle rearrangements have also been observed in 3D.^24,27–29 van der Meer et al.^30,31 identified two mechanisms of collective motion in 2D colloids: (i) closed-loops of particles and (ii) open-ended string-like particle rearrangement facilitated by pairs of high- and low-density defects, also termed vacancy (low-density) and interstitial (high density) defects. Qi et al.³² simulated the dynamics of melting in a 2D dusty plasma crystal (DPC) and found that, by mapping the real trajectory to the IS trajectory, defects can be labeled as either stable (if they persist in the IS) or unstable (if they do not persist in IS). In other words, defects in the actual trajectory may be associated either with local vibrational motions or with genuine particle rearrangement events involving irreversible particle displacements, and thus it is important to distinguish these distinct types of motions and associated defects.

Our findings both reaffirm and extend these prior observations. We show that closed-loop replacement motions of dust particles can occur without facilitation by stable defects, but these motions contribute minimally to the overall rearrangement within the crystal; these loops span a very short time duration and occur independently of each other. Much larger string-like particle rearrangement events can be identified by longer-lived energy excitations away from the ground state and are apparently facilitated by stable defect clusters. These excitation events provide the dominant contribution to displacements within the crystal. In particular, our results show that motion facilitated by only a single pair of vacancy-interstitial defect clusters (as described by van der Meer et al.^30,31) occurs only in the first excited state. Far from melting, this is the dominant mechanism for motion, but the picture becomes increasingly complex approaching the melting transition where higher energy levels become important. We show that the number of stable defect clusters grows linearly with the energy level and that the motion cannot be described solely by the creation of pairs of vacancy-interstitial defect clusters. The defects in higher energy excited states can take a variety of structures, including dipoles and trimers with differing stability. The size of both small loops and large string-like rearrangement events grows approaching the melting transition. We also find that the temperature T dependence of the mean time in the ground state between excitations, a measure of the temporal heterogeneity, directly tracks the T dependence of traditional measures of the heterogeneity time scale. In contrast, the mean time in excited states has opposite T dependence, and the time scales for ground and excited states coincide near the melting transition. Our findings offer an avenue to better understand how earlier measures of dynamic heterogeneity relate to directly observed collective displacements within our model crystalline material. At the same time, we recognize that the origin of the collective motion in the crystalline context is not necessarily the same as in amorphous glass-forming materials, but it is our expectation that insights from one system are instructive for the other.

2 Model and simulation details

We simulate 2D crystal systems using intermolecular pair potentials that mimic the properties of (i) a dusty plasma crystal^33–35 and (ii) a colloidal crystal.³¹ To carry out our simulations, we use the LAMMPS (large-scale atomic/molecular massively parallel simulator) simulation package.³⁶ We first describe the model and simulations for the dusty plasma, followed by the colloidal system.

In the DPC case, repulsive interactions among the dust grains dominate due to an excess of electrons that accumulate on the dust grain surface, which are attenuated by positive ions in the plasma. We will refer to these grains as particles to avoid potential confusion with ordered grains in semi-crystalline materials. Despite the complexity of the dusty plasma, the interactions among the particles can often be reasonably described by a simple Yukawa potential^37–40

where Q denotes the net charge on a particle and κ = 1/λ represents the screening parameter characterizing the shielding due to the plasma between two particles separated by a distance r. The dust particles are so small compared to the lattice spacing that core interactions can be neglected. To mimic the DPC studied in ref. 35, we use parameters Q = −15 700e (where e is electron charge) and κ = 2.38 mm⁻¹. The crystal state is defined by a triangular lattice with area number density of dust particles, n_d = 3.5 × 10⁶ m⁻² and a 2D Wigner–Seitz radius of mm, which corresponds to a lattice spacing of . The dust particle mass m = 5.2 ± 0.5 × 10⁻¹³ kg. We emphasize that this form of pair potential approximates dusty plasmas in three dimensions and that our choice of potential is made to apply to quasi-two dimensional dusty plasma materials.

For our simulations, we use N = 3286 particles in a nearly square box with dimensions 30.99 × 30.59 mm, and we invoke periodic boundary conditions in both directions. We carry out our simulations using the natural units of the DPC; specifically, we measure the length in units of the Wigner–Seitz radius a, mass in units of the dust particle mass m, and energy in units of U₀ = Q²/4πε₀a, which then defines time in units of , where ω_pd = 89 s⁻¹ is the inverse dust plasma frequency. Temperature is measured in units of U₀/k_B, where k_B is Boltzmann's constant and temperature is controlled using the Nosé–Hoover method. Our simulations cover a temperature range of 0.0052 ≤ T ≤ 0.0065. We use a time step dt = 0.107, and all systems are at fixed density n_d = 1/π. To provide sufficient statistics, at each temperature, we simulate 40 independent samples for 50 million time steps.

As in the experiments, melting is driven by temperature, and the dusty plasma remains crystalline for T ≲ 0.0059 in our reduced units. In simulations, it has been shown that the crystal transitions to an intermediate hexatic phase before becoming fluid; this is a subject of some controversy, particularly in the experimental system. Our intent is to focus on the dynamics in the crystal state, and not to debate the precise nature of the melting transition.

For the charged colloidal crystal case, the colloid charge interactions are attenuated by the embedding solution so that a Yukawa potential again captures most of the essential interactions. Unlike the dusty plasma, the colloidal core interactions cannot be neglected, and hence the potential can be described by,^30,41,42

where U₀ is the overall strength of the interaction and σ is the particle diameter. Following the parameters used by van der Meer et al.,^30,31 we take U₀/k_BT = 235 and κσ = 2.25. In this case, temperature is fixed and crystal melting is driven by changing the packing fraction . We use reduced units where mass is in units of the mass of a colloidal particle and length is in units of particle diameter. Like the DPC, we use N = 3286 colloidal particles and the box size varies from 145.098 × 146.997 (ϕ = 0.121) to 138.92 × 140.74 (ϕ = 0.132). To emulate the solvent forces in a colloidal suspension we use Brownian dynamics where the LAMMPS damping parameter time scale t_damp = 1.0 and the timestep dt = 0.005; we need a smaller time step than the plasma system due to the core interactions. We also considered Nosé–Hoover dynamics and found that our qualitative findings were not affected by the choice of the integration method. To provide sufficient statistics, at each packing fraction, we simulate 20 independent samples for 65 million time steps.

Now that we have defined both systems, we next examine the change in potential energy and diffusion coefficient (with temperature for the dusty plasma and with packing fraction for the colloidal system) to gain an initial understanding of the dynamics in different phases. Fig. 1(a) shows that the approximate location of melting in the DPC is evident from the increase in potential energy or the diffusion coefficient D (determined from the mean-squared displacement (MSD) 〈r²(t)〉, see Fig. S1, ESI†). Previously, van der Meer et al.³¹ reported that the colloidal system is crystalline for ϕ > 0.128, which is consistent with our simulation results. Because melting is driven by changing ϕ, which in turn alters the interparticle spacing, the potential energy is a very strong function of ϕ; consequently, the potential energy provides a less apparent indication of melting. Thus, Fig. 1(b) only shows the diffusion coefficient, which increases due to melting. Like the DPC, the fluidization appears to occur via an intermediate hexatic phase, which has been discussed in van der Meer et al.³¹ We consider other more precise metrics to identify melting in the final discussion section of the manuscript.

Fig. 1 (a) The average potential energy as a function of temperature (black) and the diffusion coefficient D (red) for the DPC. Both quantities drop sharply upon crystallization. (b) The diffusion coefficient of the colloidal system drops sharply due to crystallization. We do not include the potential energy because it is a strongly dependent function of the packing fraction, making it difficult to identify the location of the phase change.

Even in the crystalline state, it can be difficult to distinguish large local vibrations from hopping among crystal locations. To better separate these types of motion, we map the real trajectory R(t) of particles to the “inherent structure trajectory”'.²¹ The term inherent structure (IS) refers to the configuration at the nearest potential energy minimum to a given instantaneous configuration,^43,44 and in practice the IS is found by energy minimization. In LAMMPS, we minimize energy using the conjugate gradient method. The IS trajectory R^I(t) results from energy minimization of the configurations of the original trajectory R(t), yielding two parallel time series as illustrated in Fig. 2. In the IS trajectory, the system evolves from one IS to another, eliminating local vibrations. Deep in the crystal state, the IS trajectory primarily stays in the defect-free ground state with short-lived excitations due to defects; the duration of excited states increases as melting is approached, as will be discussed in detail.

Fig. 2 Schematic representation to explain how the inherent structure trajectory R^I(t) is generated from the real trajectory R(t), adapted from Schrøder et al.²¹ All dynamics occur in R(t), which are mapped to R^I(t) via energy minimization at each time t_n.

3 Results and discussion

To explore the nature of collective motion in these crystals, we first need an unambiguous way to identify those particles that have moved from their lattice positions. To do so, we calculate the probability density of particle displacements 2πrG_s(r, Δt) where G_s(r, Δt) is the self-part of the van Hove correlation function. In Fig. 3, we plot G_s(r, Δt) of the DPC for both the real R(t) and IS R^I(t) trajectories at the lowest temperature (T = 0.0052) and a higher temperature (T = 0.0058) near the melting transition at Δt = 1750; the precise value of Δt does not matter, provided it is large enough that one can distinguish distinct peaks of 2πrG_s(r, Δt). For the real particle trajectory, the distinction between the first and second peaks in the displacement distribution is difficult to identify due to lattice vibrations. In contrast, the displacement distribution in the IS trajectory shows a broad minimum from 0.7 ≲ r ≲ 1.3 where the probability is reduced by several orders of magnitude, providing a clear distinction between a particle that is localized near its starting lattice position and one that has moved to a new lattice location. The broad minimum means that we can select a displacement cutoff to distinguish localized versus displaced particles anywhere in the range 0.7 ≲ r ≲ 1.3 without qualitatively altering our findings; we choose the lower bound r_cut = 0.7 to capture as many displaced particles as possible. A similar analysis for the colloidal systems shows that a cutoff of 1.5 clearly distinguishes displacing particles from local vibrations.

Fig. 3 The self part of the van Hove correlation function for the DPC plotted in the crystal state at T = 0.0052 and T = 0.0058 (near the melting transition). The plots for T = 0.0058 are shifted vertically by 10³ for better visualization. This plot illustrates how the vibrational effects are reduced by energy minimization and gives an unambiguous choice of distance threshold to identify the hopping motion. Similar results occur for the colloidal crystal.

Having established a simple criterion to distinguish localized from displaced particles, we simply visualize the particles in the DPC that displace from their initial lattice location as a prototypical example of collective motion observed in a much broader class of materials. In previous studies, displaced particles were identified by comparing two configurations separated by a characteristic time scale for motion; Fig. 4(a) shows an example of the displaced particles using such an approach. Unfortunately, this representation loses information about when the displacements actually occur over the time interval. Because our crystals are in 2D, it is comparatively easy to show the full space-time evolution of particle displacements, as shown in Fig. 4(b), where the horizontal axis represents time. It is apparent from the image that, in the present case, the particle displacements occur primarily in two temporally distinct periods. In this sense, it is apparent that there is heterogeneity in both the spatial distribution of particle motions and the times during which motion occurs. Fig. 4(c) shows the corresponding time series for the inherent structure energy E_IS(t), demonstrating that periods of large-scale motion correspond directly to an increase in E_IS from the ground state. The same feature occurs for the colloidal system.

Fig. 4 (a) All the displaced particles in 2D within a slice of the time window of Δt = 20, (b) space-time visualization of the displaced particles in 3D over the same time interval where the time axis runs horizontal and allows one to see both where and when displacements occur. (c) The time evolution of the IS energy E_IS(t) and (d) the configuration-to-configuration displacements ΔR^I(t). The highlighted blue regions show that excited energy states correspond to extensive particle displacements, illustrated by the example clusters of red and blue particles. In addition, there are much smaller and short-lived closed-loop rearrangements that can occur in the ground state, illustrated by the green particles. Panels (a) and (b) were generated using OVITO.⁴⁵

We can relate the displacing groups more directly to the energy state by evaluating the time series for the displacement between consecutive configurations

in the IS trajectory of the N_d displaced particles. Fig. 4(c and d) affirms that transitions of E_IS(t) from the longer-lived ground state to comparatively short-lived excited states (highlighted by light-blue shading) coincide with large displacements quantified by ΔR^I(t). The figure also shows an example illustration of those particles that moved during the excited state, consistent with the observations of Fig. 4(a). In other words, the excited states correspond to large-scale particle rearrangements. While the excited state displacements are the dominant feature, ΔR^I(t) also shows very short-lived peaks in ΔR^I(t) that do not correspond to excited states. Examination of these transient peaks reveals that they result from comparatively small and nearly simultaneous closed-loop rearrangements within the ground-state crystal which occur independently. These rapid loop exchanges can occur without surmounting any energy barrier in the IS trajectory. Hence, we identify two distinct channels for motion: (i) small simultaneous loops of particles in the ground state, and (ii) much larger scale rearrangements that occur in longer-lived excited states. We point out that the distinction between motion in the ground and excited states is different from what has been observed in glass-forming materials²¹ where motion tends to occur at the transitions between IS with little or no motion within an IS.

Using our approach to identify transient particle motion, we next seek to quantify the average size and structure of groups of rearranging particles in the DPC. Since the size and structure of displacing particle groups in the ground and excited states are qualitatively different, we analyze the motions in the ground and excited states separately. First, for each excitation from the ground state (a rearrangement event), we evaluate the number of particles that have displaced at any time during the event. If a particle moves multiple times during the event, it makes multiple contributions to the overall number of displacing particles in that event. This also means that a particle may displace from and return to the same location; in practice, such returning particles are infrequent, comprising less than 3% of all displacing particles. Typically, a rearrangement event is dominated by a single large cluster. These large clusters can be decomposed into sub-groups of particles that have replaced each other. Such clusters of geometrically exchanging particles are commonly referred to as strings¹⁰ because of the quasi-linear topology of the particle exchange clusters. Two particles belong to the same string if one particle is replaced by another within a spatial tolerance of δ ≤ 0.5 (much less than the separation between neighboring particles); in addition, the displacements must occur within Δt < 1700 (less than the time for diffusion to emerge from the MSD). We tested that the spatial thresholds for replacement have a minimal quantitative effect on our findings for 0.50 < δ < 1.60. As qualitatively indicated above, this definition confirms that displacement events in the ground state only occur in closed loops of replacing particles. Hence, all groups of displacing particles, whether in the ground or excited states, can be sorted into distinct strings.

Using these definitions for the sizes of displacing particle events and the constituent string subgroups, Fig. 5(a) shows that the average size of strings 〈s〉 in the ground state grows only weakly with heating, and the behavior of the total number of displacing particles is similar; the number of displacing particles is slightly larger than the string size since, occasionally, a displacing particle group in the ground state can be decomposed into multiple strings. The average size is modest, only exceeding 10 near melting. In the excited states, the properties of the rearranging groups are quite different from those of the ground-state rearrangements. First, since the number of displacing particles in an excited state rearrangement event can be split into many strings, the average string size 〈s〉 (Fig. 5(b)) in the excited state is smaller by an order of magnitude than the number of displacing particles (Fig. 5(c)); the string size in the excited state is larger than the string size in the ground state by a factor of ≈4. In addition, the T dependence of the size of displaced particles and strings is quite different. The average string size 〈s〉 grows only modestly on heating (≈20%), similar to the scale of growth of the strings in the ground state. In contrast, the number of displacing particles grows by nearly a factor of 3, which requires that the number of strings, rather than their size, increases rapidly as we approach the melting transition. While there are collective rearrangements in both the ground and excited states, the collective rearrangements in the excited state are by far the dominant contribution to the MSD and hence the diffusion coefficient. We include the same data for the colloidal crystal in Fig. 5(d–f), and the conclusions that we made for the DPC are generally the same for this case. The growth of the size of rearranging regions on heating contrasts with the temperature dependence of the size of collective motion in the alpha-relaxation time scale of glass-forming fluids (which grow on cooling); that said, smaller collective motion at the vibrational time scale (“stringlets”) that grow on heating have been previously observed.⁴⁶

Fig. 5 The mean size of (a) groups of displaced particles (Black) and strings (Red) in the ground state, (b) strings in the excited states, and (c) groups of displaced particles in the excited states. (a)–(c) are for the DPC and (d)–(f) are for the colloidal crystal. All data are limited to the crystalline phase. In the DPC the variation of collective motion is plotted against 1/T and in the colloidal crystal, the variation is plotted against the packing fraction ϕ. The scale of collective motion becomes larger as the crystal approaches the melting transition. The lines are intended as a guide for the eye.

Having established the growing size of rearranging particle clusters as we approach melting, we next seek to quantify the average duration of the ground and excited states and how these time scales relate to other common measures of relaxation and the time scale for dynamic heterogeneity. We first evaluate the average time the system spends in the ground states t_gr and in the excited states t_ex. As might be expected, Fig. 6 shows that, on heating toward the melting transition, t_gr decreases while t_ex increases (opposite T dependence), consistent with the growing importance of excited states for melting. The distribution of t_gr values is approximately exponential, indicating that transitions from the ground to excited states are governed by a Poisson process. It is instructive to compare the ground state time t_gr to other characteristic relaxation times, including the “diffusive time scale” (D/T)⁻¹ (obtained from 〈r²(t)〉), the relaxation time τ of the self-intermediate scattering function F_s(q₀,t) (where q₀ is the wavenumber of the first peak of the structure factor), and the characteristic timescale t* of the non-Gaussian parameter α₂(t), a commonly used measure of the timescale for dynamic heterogeneity. The data for 〈r²(t)〉, F_s(q₀,t), and α₂(t) are all provided in the ESI.† These characteristic times are also plotted in Fig. 6, from which it is apparent that these times have a very similar Arrhenius temperature dependence to t_gr. That said, there are small differences in the T dependence of these characteristic times, although it is not readily apparent from Fig. 6. By parametrically plotting each of these time scales as a function of the ground-state time t_gr (see Fig. S4, ESI†), we find that t* scales linearly with t_gr. Thus, the heterogeneous distribution of ground-state times leads to t_gr proportional to t*, consistent with the expectation that t* should capture the temporal nature of heterogeneity of molecular displacements. In contrast, τ and (D/T)⁻¹ have slightly different T dependence from that of t_gr. We observe a power-law relation (D/T)⁻¹∼t_gr^1.08, and a similar result for τ. The difference of the exponent from 1 is small, so the deviation from linearity is modest, indicating a rather strong coupling among all these dynamic time scales. Previous studies^8,9,47 found that the diffusive timescale (D/T)⁻¹ and relaxation time τ varies proportionally with temperature in the crystal state (so that there is a linear parametric relation between them), and the same holds here. These qualitative observations are also true for the colloidal crystal.

Fig. 6 The variation of the relaxation time (τ - black), the average timescale spent in the ground states (t_gr - blue), the timescale to reach maximum non-Gaussian parameter α₂ (t* - red), the average timescale spent in the excited states (t_ex - green), and the diffusive timescale ((D/T)⁻¹ - magenta) plotted against 1/T. All these timescales increase as we cool down the system except t_ex.

A natural question is whether excitations from the ground state are driven by specific structural features or if they occur more stochastically. This has been explored in glass-forming systems, where Harrowell and colleagues developed the “isoconfigurational ensemble” approach to test the correlation between transient mobile regions in a liquid and the underlying fluid structure.⁴⁸ Here, we apply a similar approach to examine the temporal correlation for the transitions from the ground state to excited states. To do so, we consider configurations where a transition to the excited state has just occurred. We then chose configurations that preceded this transition by a chosen time interval. For each such configuration, we generate an ensemble of different sets of Gaussian-distributed velocities with the same temperature, hence the term isoconfigurational ensemble. We allow each member of the ensemble to evolve and track the fraction of simulations that transition to an excited state within a generous Δt = 1000. By doing this for many different intervals preceding the original excited state, we can see how rapidly the structural correlation to reach the excited state diminishes. As shown in Fig. 7, the fraction of transitioning ensemble members decreases exponentially as we increase the time preceding the original transition. As the temperature increases toward melting, the correlation time increases due to the overall higher probability of making transitions to the excited states. That said, at all T the correlation timescale is on the order of 10¹, shorter by several orders of magnitude than any other timescale discussed earlier. Thus, we can conclude that transitions at all T are dominated by stochastic rather than structural considerations, which is also consistent with the Poisson process indicated by the distribution of ground state times.

Fig. 7 The propensity of the excited states decreases exponentially for increasing time gaps preceding the excited state. The inset shows that the characteristic time of the exponential decay increases on heating primarily due to the increased likelihood of excited states at higher temperatures.

Up to this point, we have only distinguished between the ground and excited states. Examination of the excited states shows that they are separated into well-defined discrete energy bands, rather than a simple continuum of excitations. Fig. 8(a) shows a scatter plot of all E_IS values from different T and runs in a single collected time series, from which distinct energy bands clearly emerge for both the dusty plasma and colloidal systems. The right panel of Fig. 8(a) shows the histogram of E_IS values exhibiting distinct peaks for each energy band; note that the histogram is plotted on a logarithmic scale, emphasizing the separation between the bands. In Fig. 8(b), we include the same data for the colloidal crystal at different packing fractions. Because the potential energy changes strongly with the packing fraction, we normalize to align the energy of each state across the packing fractions. It is apparent that the energy bands are more sharply defined in the colloidal system than in the dusty plasma. We attribute this difference to the difference in the softness of the potentials; the DPC potential is extremely soft, since there are no core interactions, which allows for greater spread in the ground state energies, whereas the core interactions of the colloidal system make the ground state very sharply defined. The discretization of the (purely classical) crystal energies is likely a general phenomenon, as first suggested by Stillinger and Weber.⁴⁹ Nieves and Sinno⁵⁰ also observed energy quantization associated with defect formation in 3D crystalline materials but did not consider the relation to cooperative motion.

Fig. 8 Scatter plots of the E_IS time series (left panels) indicate the banding of energy minima into discrete states for (a) dusty plasma and (b) colloidal system with normalized energies . The histograms of these energies (right panels) demonstrate the clear separation among states.

It is natural to expect that the quantization of E_IS has a structural origin. Indeed, it has already been established that structural defects play a vital role in the emergence of collective rearrangements,^8,31,51,52 and these defects must affect the system energy. In the ground state, there are no defects and each particle has six neighbors. In excited states, defects are characterized by particles with either additional or missing neighbors, most commonly leading to a particle having 5 or 7 neighbors (sometimes referred to as a Stone–Wales defect).^53,54 We identify these defected particles by a Voronoi analysis. In the IS trajectory, these defects do not form in isolation but rather emerge as clusters of defects. One might naturally expect the formation of 5–7 defect pairs but defect clusters are typically more complex.

To show the connection between defects and the excited energy states, Fig. 9 shows an example time series for E_IS in parallel to the number of defect clusters (rather than simply the number of defects). These data suggest that the first excited state consists of only a pair of defect clusters and that higher energy excited states result from increasing numbers of defect clusters and their interaction. We quantify this trend by plotting the average number of defect clusters as a function of the energy state index in Fig. 10(a) for the DPC and Fig. 10(b) for the colloidal crystal. These data show an approximately linear relationship between the number of defect clusters and the energy state. These data also show that while the first excited state emerges from the creation of a pair of defect clusters, additional defect clusters are not necessarily generated in pairs.

Fig. 9 (a) Example of the time evolution of the potential energy time series (top panel) and (b) the number of defect clusters at the corresponding times (bottom panel). This plot clearly shows that there are defects only when the system is in the excited state.

Fig. 10 The variation of the average number of defect clusters in excited states for (a) dusty plasma and (b) colloid system. The insets show the variation for the average number of defects in the excited energy states has a more complex dependence. This plot shows that the energy state of the system at any particular moment is determined by how many defect clusters are present in the system.

It is not obvious that the number of defect clusters, as opposed to the overall number of defects, is the relevant parameter determining the energy state. Indeed, our initial attempts focused on simply the number of defects, which does not show a simple relation to the energy level. The number of defect clusters, rather than the number of defects, determines the energy state because defect clusters come predominantly in specific sizes.

Generally, defect clusters can be categorized by the change in local density. Most commonly, defect clusters come in pairs with one cluster having decreased density (also termed a “vacancy”) and the other having increased density (“interstitial”). Both Zhang et al.⁸ and van der Meer et al. showed that collective motion in the excited state (a reflection of what we observe in the IS trajectory) can occur via the creation, migration, and annihilation of these defects. The structure of these defects is variable in the instantaneous configurations.³¹ The structure of clusters in the minimized configurations is far simpler to understand. For the colloidal system, nearly all defect clusters are size 6 (three 5–7 defect pairs), where vacancies feature a 6-membered ring with a hole in the middle, and interstitials have an extra (non-defect) particle at the center of the ring. For the DPC where the lack of core interactions makes the potential comparatively soft, the situation is somewhat more complex. In this case, most defects are created in pairs of clusters of size 6 and 4 (two 5–7 defect pairs); in this case, clusters of size 6 are predominantly interstitials (85%), while clusters of size 4 are predominantly vacancy clusters (60%). We note that vacancy clusters in the DPC do not exhibit an explicit hole, but are simply locally less dense. In addition, the DPC also features a non-trivial fraction of clusters of size 3 (one 5–8–5 defect string) and dipoles of a single 5–7 defect pair, which we quantify in the ESI.†

Since experiments cannot readily use energy minimization to distinguish stable and unstable defects, it is helpful to note that defect structure distinguishes most unstable defects in the instantaneous configurations. Specifically, nearly all the unstable defects in the instantaneous configurations have the form of two 5–7 defect pairs (clusters of four defected particles total) which have a nearly square shape (see ESI† for an illustration). For the colloid case where all stable defects are clusters of size six, one can immediately identify these unstable 4-clusters in the instantaneous configuration. For the dusty plasma, there are a number of such 4-clusters that are stable, so one cannot eliminate all unstable clusters just based on the structure in this case. It would be valuable to examine this criterion with experimental data.

The simple picture of defects consisting of a single interstitial–vacancy pair of clusters does not hold at higher energy levels. Beyond the first excited state, dipole clusters in particular become important. This is apparent from the fact that near melting, where higher excited states are common, the fraction of defects in the form of chains of dipoles becomes substantial; this is perhaps not surprising since the chaining of dipoles is the natural mode of self-assembly in dipolar fluids.^55,56 Such geometrically polymeric dynamic structures or “strings” have often been observed in dusty plasmas and other near two-dimensional particulate crystalline materials⁵⁷ but do not seem to be predicted by conventional models of melting of two-dimensional crystals. These dipole clusters do not contribute to the collective rearrangements via “walking” through the crystal in the way that the larger 4- or 6-member clusters do. Instead, dipoles usually occur in nearby (but distinct) pairs of dipoles that quickly annihilate each other. The correlation between dipoles can be seen in Fig. 11, which shows the pair correlation function g(r) between the center-of-mass of dipole clusters, as well as g(r) for all other clusters: for dipoles, g(r) is very strongly peaked at r ≈ 4.5. In contrast, g(r) for all other clusters is nearly featureless, showing that non-dipole defect clusters behave as a nearly ideal gas within the crystal; for r ≲ 8, g(r) vanishes because defect clusters can annihilate when they are close.

Fig. 11 The pair correlation function of defect clusters excluding the dipoles (black) and the pair correlation function of the dipoles only (red). These data indicate that the dipoles are strongly correlated in space with a characteristic length δr ≈ 5, whereas all other types of defect clusters are distributed in a nearly gas-like fashion.

4 Discussion and conclusions

We have examined the spatio-temporal features of dynamic heterogeneity in soft 2D crystals which allowed us to directly connect the dynamic heterogeneity to the states sampled in the potential energy landscape. The mapping of the equilibrium trajectory to the IS trajectory of locally stable states greatly facilitates the analysis, as it eliminates trivial defects from the equilibrium crystal leaving only “stable” defect clusters that go hand-in-hand with particle rearrangements. This analysis also reveals the existence of discrete energy levels of the defects, a phenomenon first suggested in the energy landscape analysis of crystals by Stllinger and Weber.⁴⁹ It is perhaps surprising that particle rearrangements can occur even in the ground state of the crystal; these rearrangements necessarily take the form of closed loops so that no defects (which would increase the energy) are created in the process. (Technically, there must be some very short-lived defect in the structure as the rapid hop occurs, but capturing such an effect would require very frequent sampling of configurations). Open-ended particle rearrangements occur in the excited states of the crystals where defects accompany the motion. In a chicken-versus-egg discussion, one cannot meaningfully assign causation: in the IS trajectory, defects and (open-ended) rearrangements occur together and are inexorably linked. Similar large-scale particle rearrangements have been observed in experiments and connected to changes in the relaxation and modulus in both colloidal crystals²⁸ and glass-forming materials.²⁹

Overall, the defect energies somewhat resemble a quantum anharmonic oscillator where the regular spacing of energies occurs at low energies, gradually giving rise to a nearly continuous density of energies as the melting limit is approached. The observation of quantized energies of excited energy states has been reported before in simulations of ordered nanoparticles⁵⁸ and folded proteins⁵⁹ so energy quantization in these classical dynamical systems seems to be a common phenomenon at low temperatures where the material approaches a non-ergodic integrable dynamics^60–62 state at which relaxation times become formally infinite. A similar phenomenon may also be possible in glass-forming liquids at sufficiently low temperatures.⁶³

Obviously, these particle rearrangements must ultimately lead to the melting of the crystal when the number of particles in these dynamic structures becomes sufficiently large. In this sense, the balance between the defect-free ground state and the excited states where motion occurs must have a connection to melting. Fig. 6 already showed that the average time in the ground t_gr and excited t_ex states have opposite T dependence as we approach the melting transition. By extending these data to higher T, Fig. 12(a) shows that t_gr and t_ex cross T ≈ 0.00595, coinciding with the melting of the crystal. To more carefully assess this possibility, we evaluate the translational order parameter

where G is the primary reciprocal lattice vector, as well as the associated susceptibility

χ_T = N(〈|ψ_T|²〉 − 〈|ψ_T|〉²). (5)

Fig. 12(b) shows that |ψ_T| drops suddenly and χ_T has a peak at T ≈ 0.00595 in the DPC, confirming that the crystal state is no longer present. We include the same data for the colloidal crystal in Fig. 12(c and d), where the crystallinity of the system vanishes at ϕ ≈ 0.128. In 2D crystals, this is typically understood as a transition to the hexatic phase, with an additional transition to the fluid at higher T (smaller ϕ). Our point here is not to debate whether the melting occurs via an intermediate hexatic phase or a first-order crystal–liquid transition, but rather to recognize that crystal melting occurs when the average time in ground and excited states become nearly equal. In future work, it will be valuable to assess if the observed crossing of time scales is robust in more ordinary melting of 3D crystals, and how this crossing may relate to the equilibrium melting temperature or a stability limit of the crystal on superheating.

Fig. 12 (a) The T-dependence ((c) the ϕ -dependence) of the average lifetime spent in the ground state (t_gr) and excited state (t_ex). t_ex crosses t_gr at the crystal-to-hexatic phase transition in the dusty plasma (the colloid system). (b) The location of the crystal-to-hexatic transition is identified by the drop in the magnitude of the translational order parameter |ψ_T| and the corresponding peak in the translational susceptibility χ_T at T ≈ 0.00595 for dusty plasma. (d) The same thing was found for the colloid system at ϕ ≈ 0.128. The susceptibility data was adopted from the work of Meer et al.³¹

Data availability

Data used in this study will be made available upon request to the corresponding author.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank D. Leporini for helpful discussion of this work. We acknowledge funding support from NIST award 70NANB24H025.

Notes and references

1. H. Sillescu, J. Non-Cryst. Solids, 1999, 243, 81–108.
2. M. D. Ediger, Annu. Rev. Phys. Chem., 2000, 51, 99–128.
3. R. Zangi and S. A. Rice, Phys. Rev. Lett., 2004, 92, 035502.
4. A. Widmer-Cooper and P. Harrowell, J. Phys.: Condens. Matter, 2005, 17, S4025.
5. C. J. Dibble, M. Kogan and M. J. Solomon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 050401.
6. J. Kim, C. Kim and B. J. Sung, Phys. Rev. Lett., 2013, 110, 047801.
7. F. W. Starr, J. F. Douglas and S. Sastry, J. Chem. Phys., 2013, 138, 12A541.
8. H. Zhang, M. Khalkhali, Q. Liu and J. F. Douglas, J. Chem. Phys., 2013, 138, 12A538.
9. H. Zhang, X. Wang, H.-B. Yu and J. F. Douglas, J. Chem. Phys., 2021, 154, 084505.
10. C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett., 1998, 80, 2338.
11. B. A. Pazmiño Betancourt, J. F. Douglas and F. W. Starr, J. Chem. Phys., 2014, 140, 204509.
12. W. Kob, C. Donati, S. J. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett., 1997, 79, 2827.
13. Y. Gebremichael, T. B. Schrøder, F. W. Starr and S. C. Glotzer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051503.
14. H. Zhang, D. J. Srolovitz, J. F. Douglas and J. A. Warren, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 7735–7740.
15. C. R. Berardi, K. Barros, J. F. Douglas and W. Losert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041301.
16. C.-L. Chan, W.-Y. Woon and I. Lin, Phys. Rev. Lett., 2004, 93, 220602.
17. R. Xie and X.-Y. Liu, J. Am. Chem. Soc., 2009, 131, 4976–4982.
18. K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu and A. G. Yodh, Phys. Rev. Lett., 2011, 107, 108301.
19. X. Cao, H. Zhang and Y. Han, Nat. Commun., 2017, 8, 362.
20. M. Mondal, C. K. Mishra, R. Banerjee, S. Narasimhan, A. Sood and R. Ganapathy, Sci. Adv., 2020, 6, eaay8418.
21. T. B. Schrøder, S. Sastry, J. C. Dyre and S. C. Glotzer, J. Chem. Phys., 2000, 112, 9834–9840.
22. B. Cui, B. Lin and S. A. Rice, J. Chem. Phys., 2001, 114, 9142–9155.
23. C.-L. Chan, C.-W. Io and L. I, Contrib. Plasma Phys., 2009, 49, 215–234.
24. K. V. Edmond, C. R. Nugent and E. R. Weeks, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 041401.
25. H. M. Thomas and G. E. Morfill, Nature, 1996, 379, 806–809.
26. C. Yang, C.-W. Io and L. I, Phys. Rev. Lett., 2012, 109, 225003.
27. E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield and D. A. Weitz, Science, 2000, 287, 627–631.
28. J. Sprakel, A. Zaccone, F. Spaepen, P. Schall and D. A. Weitz, Phys. Rev. Lett., 2017, 118, 088003.
29. R. Higler, J. Krausser, J. van der Gucht, A. Zaccone and J. Sprakel, Soft Matter, 2018, 14, 780–788.
30. B. Van Der Meer, W. Qi, R. G. Fokkink, J. Van Der Gucht, M. Dijkstra and J. Sprakel, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 15356–15361.
31. B. van der Meer, W. Qi, J. Sprakel, L. Filion and M. Dijkstra, Soft Matter, 2015, 11, 9385–9392.
32. W.-K. Qi, Z. Wang, Y. Han and Y. Chen, J. Chem. Phys., 2010, 133, 234508.
33. U. Konopka, G. Morfill and L. Ratke, Phys. Rev. Lett., 2000, 84, 891.
34. Z. Haralson and J. Goree, Phys. Rev. Lett., 2017, 118, 195001.
35. V. Zhuravlyov, J. Goree, J. F. Douglas, P. Elvati and A. Violi, Phys. Rev. E, 2022, 106, 055212.
36. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. Int Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
37. H. Totsuji, T. Kishimoto and C. Totsuji, Phys. Rev. Lett., 1997, 78, 3113.
38. S. Hamaguchi, R. Farouki and D. Dubin, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 4671.
39. Y. Feng, J. Goree and B. Liu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 036403.
40. V. Zhuravlyov, J. Goree, P. Elvati and A. Violi, Phys. Rev. E, 2023, 108, 035211.
41. I. Rios de Anda, A. Statt, F. Turci and C. P. Royall, Contrib. Plasma Phys., 2015, 55, 172–179.
42. N. P. Kryuchkov, F. Smallenburg, A. V. Ivlev, S. O. Yurchenko and H. Löwen, J. Chem. Phys., 2019, 150, 104903.
43. F. H. Stillinger and T. A. Weber, Phys. Rev. A: At., Mol., Opt. Phys., 1983, 28, 2408.
44. F. H. Stillinger, Science, 1995, 267, 1935–1939.
45. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
46. B. A. Pazmiño Betancourt, F. W. Starr and J. F. Douglas, J. Chem. Phys., 2018, 148, 104508.
47. H. Zhang, Y. Yang and J. F. Douglas, J. Chem. Phys., 2015, 142, 084704.
48. A. Widmer-Cooper and P. Harrowell, J. Chem. Phys., 2007, 126, 154503.
49. F. H. Stillinger and T. A. Weber, J. Chem. Phys., 1984, 81, 5095–5103.
50. A. M. Nieves and T. Sinno, J. Chem. Phys., 2011, 135, 074504.
51. A. Maciołek and S. Dietrich, Rev. Mod. Phys., 2018, 90, 045001.
52. C. Reichhardt and C. Reichhardt, Phys. Rev. B, 2023, 108, 155131.
53. F. A. Lavergne, A. Curran, D. G. Aarts and R. P. Dullens, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6922–6927.
54. D. Chen, Y. Zheng, L. Liu, G. Zhang, M. Chen, Y. Jiao and H. Zhuang, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2016862118.
55. J. Stambaugh, K. Van Workum, J. F. Douglas and W. Losert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 031301.
56. K. Van Workum and J. F. Douglas, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 031502.
57. R. Quinn and J. Goree, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051404.
58. D. J. Wales and T. V. Bogdan, J. Phys. Chem. B, 2006, 110, 20765–20776.
59. N. Nakagawa and M. Peyrard, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 5279–5284.
60. L. Galgani and G. Benettin, Lett. Nuovo Cimento Soc. Ital. Fis., 1982, 35, 93–96.
61. A. Carati and L. Galgani, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 61, 4791.
62. A. Carati, L. Galgani and A. Giorgilli, Chaos, 2005, 15, 015105.
63. X. Xu, J. F. Douglas and W.-S. Xu, Macromolecules, 2022, 55, 8699–8722.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm01405g

---