Order–disorder transition in soft and deformable particle assembly with dynamic size-dispersity in two dimensions

Abstract

Soft and deformable objects are widespread in natural and synthetic systems, including micellar domains, microgel particles, foams, and biological cells. Understanding their phase behavior at high concentrations is crucial for controlling long-range order. Here, we employ a Voronoi-based model to study the packing of deformable particles in two dimensions under thermal fluctuations. Particles are represented as interconnected polygons, with the system energy comprising penalties for deviations in area and perimeter from preferred values. The strengths of these penalties capture two key features of packing: dynamic size dispersity, mimicking chain exchange in block copolymer micelles or solvent exchange in microgels, and particle line tension, reflecting the energy cost of shape changes. The model exhibits an order–disorder transition (ODT): low perimeter penalties yield disordered states, while higher penalties produce a hexagonal crystal lattice. Large dynamic size dispersity shifts the ODT to higher perimeter penalties. We explain this by analyzing particle sizes, defect formation barriers, and Voronoi entropy, which show that defect formation is easier when the area penalty term is smaller, providing a mechanistic basis for the ODT trends. In regimes far from ODT, deviations from the hexagonal lattice are accurately described by normal mode displacement fields, confirming that thermal fluctuations rather than defects govern the structure.

---

1 Introduction

Materials composed of soft and deformable constituent particles are ubiquitous in both natural and synthetic systems. Examples of such deformable particles include grafted core–shell particles,^1–4 dendrimers,^5,6 star polymers,^7,8 block copolymer micelles,^9–11 vesicles,^12,13 microgels,^14,15 emulsions,^16,17 foams,¹⁸ and biological cells.^19,20 A key characteristic of deformable particles is their ability to shrink and/or deform in response to interactions with neighboring particles under crowded conditions. For example, densely packed microgels exhibit faceting and interpenetration of polymer chains,^15,21 and can nearly eliminate interstitial voids, thereby approaching space-filling configurations. Collectively, these physical systems can be modeled as a packing of soft, deformable particles. By contrast, particles such as silica microspheres and metallic nanoparticles are significantly less deformable. The maximum packing efficiency of monodisperse hard spheres is limited to ∼74% in three dimensions, forming close-packed structures such as a face-centered cubic structure, and to ∼91% for hard discs in a hexagonal close-packed arrangement in two dimensions. The packing of hard particles has been the subject of extensive theoretical, computational, and experimental investigation over several decades.^22–24 In comparison, the packing behavior of deformable particles remains far less understood, despite its relevance to diverse soft-matter and biological systems.^15,16,21

Owing to their deformability, such particles can adapt their shape in response to external stimuli, enabling the design of adaptive materials for stretchable electronics, soft robotics, and responsive surfaces.²⁵ Thin film self-assembly studies have direct applications in advanced coatings and membranes to functional interfaces.²⁶ Moreover, quasi-2D studies, particularly through simulations and model experiments, have provided fundamental insights into jamming, collective dynamics, and melting in two dimensions,^27–29 which often differ from their three-dimensional counterparts. Two-dimensional melting, in particular, has attracted considerable recent attention, especially with respect to the possible existence of intermediate hexatic phases and the nature of the melting transition, including continuous versus weak first-order scenarios.^28,30 Recent efforts have also focused on understanding the packing behavior of hard disks with static polydispersity.^31,32 In contrast, the impact of dynamic size dispersity on the packing structure remains relatively underexplored, despite its widespread relevance in experimental systems where particle sizes can fluctuate. For instance, the size of block copolymer micelles can fluctuate due to the exchange of copolymer chains between micelles,³³ and microgel particles can swell or de-swell in a solvent, thereby altering their size within dense packings.³⁴

To understand packing behavior, it is instructive to distinguish between packing and tiling, a distinction that becomes essential in nearly confluent assemblies of deformable particles. Packing refers to the arrangement of objects within a space, whereas tiling describes the partitioning of that space into distinct, non-overlapping compartments. In nearly confluent systems of densely packed deformable particles, the distinction between packing and tiling vanishes, as particle deformation allows them to collectively fill space without significant voids. In this regime, many-body interactions naturally emerge,^1,35 breaking down the common simplifying assumption of pairwise additivity that is often employed in computational simulations. To capture these collective effects, several modeling frameworks have been developed, including machine learning approaches for polymer-grafted nanoparticles,^1,2,35 deformable particle models representing particles as 2D spring rings,³⁶ and Vertex^19,37–39 and Voronoi-based^40,41 models commonly used to describe epithelial monolayers with polygonal cells.

In this study, we adopt a Voronoi-based model to investigate the packing behavior of soft and deformable particles in 2D where Voronoi cells represent the faceted polygonal particles under dense packing. We sample the system configurations from Boltzmann distributions at a fixed temperature, quantifying the effect of thermal fluctuations on the self-assembly. We further examine the impact of dynamic size dispersity on phase behavior, mimicking experimental scenarios in which the size of constituent particles fluctuates. We quantify topological defects in such systems and show their effect on the order of the system. Finally, we perform normal mode analysis of the ordered structures and discuss the resemblance between our model and the harmonic spring network.

2 Materials and methods

2.1 Voronoi model

We adopt a Voronoi-based model to simulate the 2D assembly of densely packed (no gaps between particles) soft and deformable particles. The phase behavior of such dense packings has been explored previously in three dimensions.^42,43 In the model, the Voronoi diagram represents the system configuration with Voronoi cells describing the faceted deformed particles. The Voronoi and Vertex models are closely related. While polygon centers serve as the degrees of freedom in the Voronoi model, polygon vertices are the degrees of freedom in the Vertex model. The energy of the system in our Voronoi model is given by the following expression:where N is the total number of particles. The first term in the energy equation captures the area elasticity and imposes an energy penalty when the instantaneous area A_i of the particle i deviates from its preferred area A₀ and k controls the strength of the penalty, thereby controlling the degree of dynamic area-dispersity (size-dispersity) of the assembly in the system. The second term quantifies the perimeter penalty on the particle whose instantaneous perimeter P_i deviates from the preferred perimeter P_ci which is taken as the perimeter of an equivalent circle of area A_i, i.e., . The strength of the perimeter penalty is controlled by the line tension, γ. The choice of P_ci reflects that an isotropic circle geometry has the minimum perimeter density for the lowest line tension energy. The parameters k and γ capture two key aspects of packing behavior in soft, deformable particles. In experiments, k scales the energy penalty by the particle size fluctuations—accounting for chain exchange in block copolymer micelles and the resulting size dispersity—while γ quantifies particle deformability via the energy cost of shape changes. For example, in emulsion systems, γ can be mapped onto the effective interfacial tension of emulsion droplets, while in polymeric microgels or block copolymer micelles, γ reflects the energetic cost of particle deformation arising from polymer stretching.

We nondimensionalize the quantities in eqn (1) by scaling energies with thermal energy k_BT, where k_B is Boltzmann's constant and T is absolute temperature, and by scaling lengths with A₀^1/2. In nondimensional form, eqn (1) becomes

where Ẽ = E/k_BT, k̃ = kA₀²/k_BT, Ã_i = A_i/A₀, [small gamma, Greek, tilde] = γA₀^1/2/k_BT, P̃_i = P_i/A₀^1/2, and P̃_ci = P_ci/A₀^1/2. Quantities denoted with tildes are dimensionless, and this convention is used consistently throughout the paper.

2.2 Monte Carlo simulations

We performed Monte Carlo (MC) simulations in the canonical ensemble for a two-dimensional system of N = 3600 particles. The simulation box dimensions were set to 55.836 × 64.476, ensuring a preferred area of one in dimensionless units per particle. Also, the box dimensions are commensurate with a hexagonal crystal lattice. Periodic boundary conditions were applied along both the x- and y-directions. Each Monte Carlo move consisted of selecting a particle at random and displacing it by a random vector with magnitude up to a prescribed maximum displacement. After every move, we constructed the Voronoi tessellation of the particle configuration, using particle positions as the centers. This was implemented with MATLAB's “delaunayTriangulation” library.⁴⁴ The system energy was evaluated from the Voronoi diagram using eqn (2), and equilibrium configurations were sampled via the standard Metropolis algorithm. We refer to these as “thermal systems,” reflecting the inclusion of thermal energy in the simulations. Based on eqn (2), increasing k̃ or [small gamma, Greek, tilde] reduces the relative thermal energy in the system and vice versa. The maximum displacement was chosen such that the acceptance ratio lies in the range of 30–60%. We equilibrated our system for at least 40 000 Monte Carlo cycles where one cycle consists of N displacement moves. At least 100 simulation snapshots, evenly sampled from 20 000 Monte Carlo cycles in the production run, were used to calculate the ensemble average of thermodynamic properties. Errors were estimated as the standard deviation of five block averages from the mean.

We investigate the phase behavior of particles interacting with the energy functional expressed in eqn (2) in three different regimes of dynamic size-dispersity. To do so, we choose three different values of k̃ to obtain different degrees of dynamic size-dispersity in the system, k̃ = 1000, k̃ = 100, k̃ = 30, to ensure low, moderate, and high dynamic size-dispersity, respectively. The initial configurations for Monte Carlo simulations were prepared in a perfect hexagonal lattice, oriented such that the two farthest-separated vertices of any hexagon lie along the x-axis. Three series of simulation runs were carried out with k̃ = 1000, k̃ = 100, k̃ = 30, each over a range of [small gamma, Greek, tilde] values. Since energy is scaled by the thermal energy, a decrease in [small gamma, Greek, tilde] corresponds to an increase in the relative thermal energy of the system. To examine the system-size effect, we also simulated systems with N = 900 and N = 2500 in the boxes with dimensions commensurate with a hexagonal crystal lattice for the case of k̃ = 100 at varying values of [small gamma, Greek, tilde] (see Fig. S1 in the SI). The system behavior appears to be independent of system size. We chose a moderate system size of N = 3600 to focus on the two dominant phases of interest, ordered and disordered, that are robustly captured within this regime. We also assessed the influence of initial conditions on the equilibrium phases (see Fig. S2 in the SI). Simulations initialized from both a perfect hexagonal lattice and a configuration obtained after equilibrating at nearby values of [small gamma, Greek, tilde] converge to the same equilibrium state, demonstrating robustness to initial conditions enabled by sufficient Monte Carlo sampling.

2.3 Structure factor

We characterize the system structure by computing the structure factor S(q⃑) of the particle assembly bywhere is the position vector of the Voronoi center of the ith particle, N is the total number of particles in the system, is a wave vector, and 〈⋯〉 denotes an average over at least 100 evenly spaced simulation snapshots. The allowed wave vectors were determined by discretizing the simulation box spatially into 300 × 300 grids. Then, we obtain azimuthally averaged S(q̃) from where q̃ is the magnitude of . For this, q̃ was assigned to discrete radial bins spanning the desired q̃-range. Within each bin, the corresponding values were averaged to yield the powder-averaged S(q̃).

2.4 Orientational order parameter

We characterized the system structure by quantifying the hexatic order in the system since the honeycomb structure of uniform hexagons has the least line tension energy.⁴⁵ We compute the local bond-orientational order parameter ψ₆ for each particle using eqn (4).In the eqn, i is the particle identity, z_i is the coordination number of the ith particle (same as the number of vertices of the particle), θ_ij is the angle between the positive x-axis and the vector joining the Voronoi center of the particle to the vertex j. We refer to θ_ij as the vertex angle. In the expression, 6 denotes the six-fold (hexatic) symmetry of the regular hexagon. By construction, ψ_6,i is a complex number, where the magnitude |ψ_6,i| quantifies the local six-fold order of the particle, and the phase arg(ψ_6,i) indicates the orientation of the particle relative to the positive x-axis. After computing ψ₆ for all the particles, we obtain the mean global bond orientational order parameter Ψ₆ of a system configuration by using eqn (5) where 〈…〉 denotes an average over at least 100 evenly spaced simulation snapshots. For our initial configuration where two farthest-separated vertices of hexagons lie along the x-axis, Im(Ψ₆), the imaginary component of Ψ₆, is zero, and Re(Ψ₆), the real component of Ψ₆, is one. Moreover, for a completely disordered arrangement of polygons under the constraint where the mean coordination number of polygons is 6, Ψ₆ is essentially zero.

2.5 Voronoi entropy

We also compute the Voronoi entropy S_vor which is the Shannon entropy of the Voronoi tessellations to compute the structural order in the system defined as follows:

S_vor = 〈−zP_z ln P_z〉 (6)

where P_z denotes the fraction of Voronoi polygons with coordination number z, and 〈…〉 denotes an average over at least 100 evenly spaced simulation snapshots. This entropy measures the deviation of the observed Voronoi cell distribution from that of an ideal single crystal, and therefore vanishes for a perfect hexagonal arrangement. While orientational order parameters and structure factor calculations primarily quantify long-range structural order, Voronoi entropy provides a complementary measure that captures the degree of disorder in the system.

2.6 Athermal defect formation simulations

We map the energy landscape of defect formation in a system of N = 400 particles with no thermal fluctuations. To introduce 5–7–7–5 bound dislocation pairs, we select a 6–6–6–6 four-particle cluster composed of two neighboring particles along the y-axis and their two immediate neighbors along the x-axis. The system is initialized in a perfect hexagonal configuration. We then actively modulate the interparticle distances l̃_x and l̃_y, defined as the separations between opposite particles in the x and y directions of the chosen four-particle cluster, while allowing all surrounding particles to relax in order to minimize the total energy. Specifically, we vary l̃_y first while keeping l̃_x fixed, and then vary l̃_x for each chosen value of l̃_y. We compute the energy of the relaxed system in the l̃_x − l̃_y plane.

2.7 Normal mode analysis

We perform normal mode analysis to characterize the collective vibrational (phonon) modes of a system of N = 3600 particles arranged in a perfect hexagonal configuration. To do so, we construct the Hessian matrix from the second derivatives of the total energy with respect to particle displacements using the following expression:where the indices i/j are particle identity, μ/ν are the Cartesian coordinate directions, ũ is particle displacement, and Ẽ is the total energy computed using eqn (2). We use a finite difference method to quantify the second derivatives. The Hessian matrix is of the size 2N × 2N for two-dimensional system. The Hessian eigenvectors V_n correspond to the normal modes denoted by n, describing the relative directions and amplitudes of motion for each particle, while we associate the eigenvalues [small lambda, Greek, tilde]_n with the squared frequencies [small omega, Greek, tilde]_n² of these collective vibrations i.e., . Since our physical system with thermal fluctuations is overdamped, this analysis instead provides the eigenvectors and eigenfrequencies of a corresponding “shadow system” that retains the same interparticle interactions but neglects damping.

We also characterize normal modes of an analogous system of a 2D elastic spring network connecting point masses forming the hexagonal crystal lattice. The springs are Hookean and only connect the nearest neighbors. The total energy of the spring network Ẽ_spring is computed using

where i′ is the spring identity, 2N is the total number of springs connecting N point masses, k̃_spring is spring constant and set to one, is the instantaneous length of the spring i′, and is the equilibrium spring length corresponding to the nearest-neighbor distance in the ground state of the hexagonal crystal lattice. We follow the same procedure as outlined above to obtain information about different collective modes.

For each mode n, the displacement amplitude was drawn from a normal distribution with zero mean and variance 〈ã_n²〉 = 1/[small lambda, Greek, tilde]_n, consistent with equipartition. Excluding near-zero eigenvalues associated with translational modes, the total displacement field was constructed as a linear combination . Particle positions were then updated by displacing particles by ũ from the perfect hexagonal crystal lattice. Repeated sampling produced an ensemble of imperfect hexagonal crystal configurations, which allowed us to compute Ψ₆.

In a nearly perfect hexagonal crystal, particle displacements give rise to minor distortions of the Voronoi cells from regular hexagons. To quantify these distortions, we computed the probability distribution of the vertex angles θ_ij. The resulting distribution is hexamodal, with each peak centered around (j − 1)π/3, where j ∈ [1, 6] indexes the six-fold symmetry directions. A representative probability distribution is shown in Fig. S3 in the SI. The standard deviations of the six modes are nearly identical, and their average value, denoted as δ, is used as a measure of the overall angular distortion in the system configurations. We compute δ for the imperfect hexagonal crystal configurations obtained from normal mode analysis and for k̃ = 1000, k̃ = 100, and k̃ = 30 thermal systems in the large-γ limit. Building on this angular-deviation framework, we next examine its implications for the bond-orientational order parameter. Using eqn (4), we calculate the local bond-orientational order parameter of a single distorted Voronoi cell by sampling vertex angle, θ_j, from a normal distribution with mean (j − 1)π/3, where j ∈ [1, 6] and variance δ². The global bond-orientational order parameter Ψ₆ is then evaluated by averaging local bond-orientational order parameters over at least 1000 individual Voronoi polygon configurations.

3 Results and discussion

3.1 Order–disorder transition

Our model system shows an order–disorder transition (ODT) when [small gamma, Greek, tilde] is varied at a fixed k̃. Fig. 1a and b shows the representative simulation snapshots of the k̃ = 1000 system equilibrated at [small gamma, Greek, tilde] = 8 and [small gamma, Greek, tilde] = 6, above and below the ODT [small gamma, Greek, tilde] point, respectively. The 2D scattering profile for each state is also shown in the bottom right corner of each snapshot. In the case of [small gamma, Greek, tilde] = 8, the presence of six Bragg peaks exhibiting a six-fold symmetry indicates that the system is in a hexagonally ordered phase, with particle Voronoi centers arranged on a hexagonal crystal lattice. On the other hand, at [small gamma, Greek, tilde] = 6, we see an isotropic ring in the scattering profile which indicates a disordered phase with no crystal symmetry associated with the arrangement of particles.

Fig. 1 Representative simulation snapshots of the k̃ = 1000 system at (a) [small gamma, Greek, tilde] = 8 and (b) [small gamma, Greek, tilde] = 6. The simulated scattering profile is also shown in the bottom right corner of each figure. Particles are color-coded based on their coordination numbers z.

To obtain an estimate of the location of the ODT in the [small gamma, Greek, tilde] parameter space at a fixed k̃, we compute the azimuthally averaged structure factor S(q̃) from the 2D scattering pattern (see Section 2.3 in Materials and methods). Fig. 2a shows the qualitative comparison of S(q̃) obtained at three values of [small gamma, Greek, tilde] near the transition. At large [small gamma, Greek, tilde], for example [small gamma, Greek, tilde] = 10, S(q̃) has the clear primary peak at q̃₁ ≈ 6.8, secondary peak at , and so on, confirming the existence of the underlying hexagonal arrangement of particles. As we decrease [small gamma, Greek, tilde] (increasing relative thermal energy) of the system, the peak heights in S(q̃) decrease gradually indicating the weakening of the hexagonal order in the system. When [small gamma, Greek, tilde] is small enough (relative thermal energy is high enough), we see a jump in the peak heights with respect to [small gamma, Greek, tilde] at the ODT as the system transitions into the disordered phase from the ordered phase. The persistence of finite peak heights even when the system is globally disordered indicates a short-range liquid-like packing, which is due to the excluded size interaction between particles (i.e., k̃ is non-zero). This trend can be seen more clearly if we track the primary peak height extracted from S(q̃) with [small gamma, Greek, tilde] (Fig. 2b). Here we have normalized the S(q̃) of primary peak heights with that of the primary peak height of a perfect hexagonal crystal lattice. We also characterize the ODT in the system using the mean global orientational order parameter, Ψ₆. Since Im(Ψ₆) is negligible in both the disordered phase and the ordered phase, due to the uniform distribution of ψ₆ and the chosen crystal orientation of the initial configuration, respectively, we focus our analysis on the behavior of Re(Ψ₆) (Fig. 2b). When [small gamma, Greek, tilde] is decreased from a large [small gamma, Greek, tilde], say, [small gamma, Greek, tilde] = 100, Re(Ψ₆) monotonically decreases and corroborates the observation of weakening of the hexagonal order in the system from above analysis. As the system disorders from the ordered state, Re(Ψ₆) jumps from 0.25 to approximately zero in the disordered phase, indicating a complete loss of hexagonal symmetry. We estimate the location of the ODT by taking the average of the largest [small gamma, Greek, tilde] at which the system is disordered and the smallest [small gamma, Greek, tilde] at which the system is ordered obtaining [small gamma, Greek, tilde]_ODT = 6.5 for the k̃ = 1000 system. We further characterized the order–disorder transition using the bond-orientational correlation function (see the SI for details). The disordered phase exhibits finite bond-orientational correlation lengths, whereas the ordered phase displays long-range correlations (refer to Fig. S4 in the SI). In conclusion, our simple model shows an ODT behavior. We note that studies on such 2D systems have shown the existence of a hexatic phase – sandwiched between disordered and hexagonally ordered phases – with a narrow stability region.^28,46,47 However, we do not observe a hexatic phase in the present study likely due to the moderate system size considered. In two-dimensional systems, intermediate hexatic phases are known to be particularly sensitive to the system size effect,^30,48 and we therefore cannot rule out the possibility that a hexatic phase may emerge in larger systems. With this caveat in mind, [small gamma, Greek, tilde]_ODT is an estimate of location of ODT to a first approximation.

Fig. 2 (a) Azimuthally averaged structure factor S(q̃) as a function of q̃ for the k̃ = 1000 system at three values of [small gamma, Greek, tilde]. (b) Normalized primary peak heights of S(q̃) relative to that of a perfect hexagonal lattice (solid line with closed markers) and mean global orientational order parameter (solid line with open markers) with [small gamma, Greek, tilde]. Both quantities jump at the ODT, which is estimated to be at [small gamma, Greek, tilde]_ODT = 6.5.

In our model, k̃ sets the degree of dynamic size-dispersity in the system at a given [small gamma, Greek, tilde] with smaller k̃ leading to greater size-dispersity among particles. To investigate the impact of dynamic size-dispersity on the ODT, we take k̃ = 30 (high size-dispersity) and k̃ = 100 (moderate size-dispersity) along with k̃ = 1000 (low size-dispersity) systems. The behavior of Re(Ψ₆) with [small gamma, Greek, tilde] for these systems is shown in Fig. 3a. We estimate size-dispersity in our system by the standard deviation in particle area from its mean preferred value (which is one in dimensionless units in our simulations), denoted here as [small sigma, Greek, tilde]_A. Fig. 3b shows the variation of [small sigma, Greek, tilde]_A with [small gamma, Greek, tilde] for three systems. Similar to the k̃ = 1000 system, the other two systems with k̃ = 100 and k̃ = 30 also show the ODT when [small gamma, Greek, tilde] is modulated (Fig. 3a). Apart from that, we also observe a monotonic decrease in [small gamma, Greek, tilde]_ODT with k̃. The estimate of [small gamma, Greek, tilde]_ODT is 17.5 for k̃ = 30 and 10.5 for k̃ = 100 systems along with 6.5 for the k̃ = 1000 system. The bond-orientational correlation function analysis for the three systems is also shown in Fig. S4. Qualitatively, the trend can be explained by considering two effects. First, increasing dynamic size-dispersity in our model enhances the configurational entropy of the disordered phase, as the disordered system can access an enormous number of distinct Voronoi tessellations. This entropic gain stabilizes the disordered phase relative to the ordered phase shifting the ODT. A similar stabilization mechanism has been observed in experiments⁴⁹ and simulations⁵⁰ of three-dimensional colloidal systems with static polydispersity, as well as in computational studies of hard disks with static polydsipersity,³² where permutations of particle sizes introduce additional degrees of freedom that increase the configurational entropy of the fluid phase.^51,52 Second, for our system at a given k̃, increasing [small gamma, Greek, tilde] reduces size-dispersity (Fig. 2b). This result is a direct consequence of our energy expression (eqn (2)) where the ground state is the perfect hexagonal crystal lattice. At a given k̃, increasing [small gamma, Greek, tilde] moves the system closer to the ground state, while simultaneously driving the particle size and shape closer to a regular hexagon of unit area (in dimensionless units). Thus, in our model, the ordered phase is energetically favored while the disordered phase is entropically favored. Taken together, a higher [small gamma, Greek, tilde]_ODT is observed for systems with lower k̃.

Fig. 3 (a) Mean global orientational order parameter and (b) dynamic size-dispersity between particles for k̃ = 1000, k̃ = 100, and k̃ = 30 systems with [small gamma, Greek, tilde].

At a given k̃, [small gamma, Greek, tilde] not only affects the size-dispersity but also the shape of the particles. To gain a deeper understanding of the microscopic mechanisms underlying the transition, we computed 〈|ψ₆|〉 which is a measure of the mean degree of hexagonality of particles in the assembly and show its dependence on [small gamma, Greek, tilde] in Fig. 4a. We also computed 〈|ψ₆|〉 for random Voronoi tessellation where Voronoi centers were generated randomly (Poisson distribution) in the simulation box. In line with Fig. 3, 〈|ψ₆|〉 shows a jump with respect to [small gamma, Greek, tilde] at the ODT; however, the jump is less sharp. Below the ODT in the disordered phase, these three systems merge onto the common branch of 〈|ψ₆|〉 with [small gamma, Greek, tilde]. As [small gamma, Greek, tilde] is decreased further and approaches zero, 〈|ψ₆|〉 for all the system saturates to about 0.36 which is close to 0.37 computed for random Voronoi tessellation implying that the arrangement of particles is close to random. The difference is due to the fact that k̃ is non-zero for three systems leading to excluded-area interactions among particles. The arrangement of particles being approximately random at low [small gamma, Greek, tilde] is also verified by comparing their |ψ₆| distribution against that of random Voronoi tessellation shown in Fig. 4b. The signature of ODT is also evident from the pronounced changes in the |ψ₆| distribution across the transition, reflecting the collective rearrangement of particle shapes as the system approaches ODT. Due to the deformable nature of particles and thermal fluctuations, even deep into the ordered phase, a non-negligible population of particles adopt irregular hexagonal shapes as indicated by the shape distribution for the k̃ = 30 system at [small gamma, Greek, tilde] = 90.

Fig. 4 (a) Mean degree of hexagonality of particles 〈|ψ₆|〉 with γ for k̃ = 1000, k̃ = 100, k̃ = 30, and random Voronoi tessellation. Error bars are smaller than the marker size. (b) Probability density function (pdf) of |ψ₆| for the k̃ = 30 system at [small gamma, Greek, tilde] = 2, at [small gamma, Greek, tilde] = 17 in the disordered phase, and at [small gamma, Greek, tilde] = 18 and [small gamma, Greek, tilde] = 90 in the ordered phase along with the distribution of |ψ₆| for random Voronoi tessellation.

Particle shape and size distributions both influence the ODT in deformable particle systems. Notably, the 〈|ψ₆|〉 for the k̃ = 30 system at [small gamma, Greek, tilde] = 17 in the disordered phase exceeds that of the k̃ = 1000 system at [small gamma, Greek, tilde] = 7 in the ordered phase, shown in Fig. 4a. This suggests that large size dispersity can suppress long-range order even when many particles exhibit a high |ψ₆|. An open question is: what is the upper limit of 〈|ψ₆|〉 that still permits disorder due to size dispersity? The coupled effects of particle size and shape on phase behavior will be examined in future work.

We also characterize the structural order in the system by measuring the Voronoi entropy of the particle assembly. This analysis has proven useful in analyzing 2D point patterns.^53,54 We emphasize that S_vor is an intensive quantity, whereas thermodynamic entropy is extensive. In Fig. 5, we report S_vor as a function of [small gamma, Greek, tilde] for three systems, along with results for a random Voronoi tessellation. For the random Voronoi tessellation, we compute S_vor ≈ 1.69, consistent with values reported in the literature.^54,55 Consistent with the trends discussed above, S_vor exhibits discontinuous changes at the ODT: the disordered phase is characterized by a higher entropy, while the hexagonal ordered phase exhibits a lower entropy. At the same [small gamma, Greek, tilde], systems with a greater dynamic size dispersity (or lower k̃) display a higher entropy, indicating that size dispersity promotes structural disorder. Random Voronoi tessellation serves as the limiting behavior of our system indicating that lowering k̃ and [small gamma, Greek, tilde] drive the system closer to the random Voronoi tessellation.

Fig. 5 Voronoi entropy S_vor with [small gamma, Greek, tilde] for k̃ = 1000, k̃ = 100, and k̃ = 30 systems, and random Voronoi tessellation.

3.2 Defect formation

The structural disorder in the ordered phase in our system arises from the presence of topological defects. A defect is defined as a particle whose number of sides (or neighbors) is different from 6. Since particles in 2D assembly are less correlated compared to their 3D analog, defects spontaneously form and annihilate in the assembly, and there exists a non-zero equilibrium concentration of defects in the system.^27,56 These defects are generated by structural rearrangement of particles in succession such that two adjacent particles swap their neighbors, called T1 topological events/transitions. A T1 event transforms a four-particle cluster of six-sided particles (6–6–6–6) into a cluster of 5–7–7–5-sided particles, forming what are known as bound dislocation pairs. The introduction of such dislocation pairs into a crystal creates singularities in the underlying local order parameter field, thereby disrupting the local structural order. We estimate the concentration of topological defects in the ordered phase, f_def, by calculating the fraction of non-hexagonal particles, which is the ratio of non-hexagonal particles to the total number of particles. Topological defects are not defined in the static disordered configurations due to the absence of an underlying crystal lattice; however, the fraction of non-hexagonal particles can still be quantified. Fig. 6 shows the variation of f_def with [small gamma, Greek, tilde]. Far from the ODT in the ordered phase, we see a negligible concentration of defects in all the systems; however, the defect concentration exponentially increases (inset of Fig. 6) as the system approaches the ODT followed by melting of the crystalline order.

Fig. 6 Concentration of topological defects f_def with [small gamma, Greek, tilde] for k̃ = 1000, k̃ = 100, and k̃ = 30 systems. The inset shows the linear fit on ln f_def with [small gamma, Greek, tilde] for three systems. The solid markers are the data points, and dashed lines are the best fit.

This observation is reminiscent of the defect-mediated melting in 2D described by the Kosterlitz–Thouless–Halperin–Nelson–Young theory.^27,57,58 The exponential dependence of defect concentration on [small gamma, Greek, tilde] indicates that defect formation and annihilation follow Arrhenius-type behavior, with the slope of the best-fit line providing an estimate of the activation barrier. Similar behavior has been reported in neighbor exchange events in a disordered cell monolayer studied using the Vertex model.⁵⁹ From the inset, the slopes follow the order k̃ = 1000 > k̃ = 100 > k̃ = 30, implying that defect formation becomes easier as k̃ decreases.

To gain insights into defect formation in our thermal system, we take a simpler system of a perfect hexagonal crystal with no thermal fluctuation. We then perform a T1 transition to create 5–7–7–5 dislocation pairs. To do so, we modulate l̃_x and l̃_y, interparticle distances of opposite particles in the x and y directions in the four-particle configuration shown in the top panel of Fig. 7c, respectively, and minimize the energy of the surrounding particles (see Section 2.6 in Materials and methods). This procedure allows us to create defects in a controlled manner while looking at the underlying energy landscape of defect formation. We compute contour maps of the system energies relative to the ground state, i.e., the total energy minus the energy of the ground state of the hexagonal crystal lattice, in the l̃_x − l̃_y plane for the three cases of k̃ at [small gamma, Greek, tilde] = 20. These systems are ordered under these conditions when thermalized (see Fig. 3). Both k̃ = 30 and k̃ = 100 systems show a metastable defect-state with respect to the perfect hexagonal state (global energy minimum); however, we do not see any signature of metastable defect-state for the k̃ = 1000 case (refer to Fig. S5 in the SI for energy contour maps for k̃ = 100 and k̃ = 1000). We extract the energy barrier heights at the T1 transition as a function of l̃_x along the l̃_x = l̃_y line for the three cases of k̃. Note that l̃_x is equal to l̃_y at the T1 transition. We show the comparison of energy barriers in Fig. 6b. The minima of the energy barrier curves denote the intersection of the minimum-energy paths (MEPs) taken by the system to form defects starting from the perfect hexagonal state. We observe that the energy barrier height along the MEP increases with k̃ suggesting that it becomes increasingly improbable to form defects when k̃ is increased. Moreover, the curvature of the energy barrier around its minima is also affected by the value of k̃. A large value of k̃ sets a higher curvature of energy barrier around the MEP indicating that other energy pathways taken by the system on the energy landscape to form defects through T1 transition become less probable with increase in k̃. This analysis likely suggests that a larger k̃ suppresses the defect formation in our thermal systems at a given [small gamma, Greek, tilde], giving a potential explanation why the ODT for larger k̃ values occurs at smaller [small gamma, Greek, tilde] values.

Fig. 7 (a) Contour maps of the system energies relative to the ground state, i.e., the total energy minus the energy of the ground state of the hexagonal crystal lattice, in the l̃_x − l̃_y plane at [small gamma, Greek, tilde] = 20 and k̃ = 30. The dashed line is drawn at l̃_x = l̃_y which is at the T1 transition. (b) Total energy minus the energy of the hexagonal state at the T1 transition as a function of l̃_x where l̃_x = l̃_y at [small gamma, Greek, tilde] = 20 for k̃ = 30, k̃ = 100, and k̃ = 1000. The closed circles indicate the minimum of the curve. (c) Top panel A represents the perfect hexagonal state, middle panel B is the representative configuration at the T1 transition, and bottom panel C is a representative configuration of a defect-state.

3.3 Far from the ODT

In a crystalline solid at finite temperature, atoms are not stationary but undergo collective vibrations, known as phonons, with frequencies determined by the interatomic potential and thermal energy. We expect analogous phonon-like behavior in our system, particularly far from the ODT in the ordered phase. Because our system is overdamped, this should manifest as similar displacement fields in the particle configurations. To probe these modes, we perform normal mode analysis on our system (see Section 2.7 in Materials and methods). For a two-dimensional crystal, the Debye model predicts that the cumulative number of normal modes scales with the square of the frequency in the low-frequency regime.⁶⁰ Our model conforms with the Debye model in large-[small gamma, Greek, tilde] limit as shown in Fig. S6 in the SI. In this limit, we also expect that the many-body interactions will be negligible since particles vibrate around the lattice positions with small displacements, and the coupled interparticle interactions in our model will effectively converge to harmonic interactions. In other words, the Voronoi model should behave as a harmonic spring network model in the large-[small gamma, Greek, tilde] regime. To test this hypothesis, we compare the phonon density of states, D([small omega, Greek, tilde]), of our crystal system with that of a harmonic spring network, as shown in Fig. 8. We observe two divergent peaks in D([small omega, Greek, tilde]) for the Voronoi model, consistent with the prediction of two van Hove singularities in D([small omega, Greek, tilde]) for a two-dimensional hexagonal crystal, corresponding to the longitudinal and transverse modes.⁶¹ A similar feature is observed in D([small omega, Greek, tilde]) for the harmonic spring network model.

Fig. 8 Phonon density of states, D([small omega, Greek, tilde]), with normal mode frequency, [small omega, Greek, tilde] for (a) the Voronoi model with k̃ = 100 and [small gamma, Greek, tilde] = 500 and (b) a harmonic spring network with pair-wise interactions with a spring constant of 1 (in dimensionless units).

Furthermore, we estimated the displacements of all particles in the Voronoi system from their corresponding lattice positions using the Hessian eigenvalues and eigenvectors (see Section 2.7 in Materials and methods). Applying these displacement fields to the perfect hexagonal lattice generated ensembles of perturbed configurations, from which Ψ₆ was computed. In the absence of defects, thermal energy in an undamped system drives collective oscillations of particles, and similarly, we expect collective particle displacements in our thermal systems far from the ODT, i.e., [small gamma, Greek, tilde] > 50—where the defect concentration is negligible—to be governed purely by thermal energy. To test this, we compared Re(Ψ₆) obtained from normal mode analysis with the results from thermal simulations at k̃ = 1000, k̃ = 100, and k̃ = 30 (Fig. 9), where the x-axis represents the overall angular distortion in the system configurations, denoted as δ, defined as the average of the six standard deviations corresponding to the six modes in the hexamodal vertex angle probability distribution (see Section 2.7 in Materials and methods). The close agreement between the thermal systems and normal mode analysis indicates that all k̃-systems in the large-[small gamma, Greek, tilde] limit exhibit collective particle displacements consistent with the predictions of normal mode theory. We further examined a case of a single Voronoi cell in which vertex angles are sampled from normal distributions with mean (j − 1)π/3 where j ∈ [1, 6] and variance δ², which also reproduced the observed behavior, supporting the conclusion that vertex angular deviations in thermal systems are well described by a normal distribution in the large-[small gamma, Greek, tilde] limit.

Fig. 9 Mean global orientational order parameter with δ, defined as the average of the six standard deviations corresponding to the six modes in the hexamodal vertex angle probability distribution, for the system configurations of k̃ = 1000, k̃ = 100, and k̃ = 30, for configurations generated by applying displacement fields to a perfect hexagonal lattice from normal mode analysis. Note that [small gamma, Greek, tilde] is variable along the curve for thermal systems with different k̃ values. Error bars are smaller than the marker size. The dashed line represents the results obtained for a single distorted Voronoi cell, where the vertex angles, θ_j, were sampled from a normal distribution with mean (j − 1)π/3, for j ∈ [1, 6] and variance δ².

4 Conclusions

In this study, we use a Voronoi model to investigate the packing behavior of soft, deformable particles in two dimensions. Two model parameters, [small gamma, Greek, tilde] and k̃, control particle line tension and dynamic size-dispersity, respectively. Monte Carlo simulations at constant temperature reveal an order–disorder transition: the system forms a hexagonal crystal lattice at high [small gamma, Greek, tilde] and melts into a disordered phase at low [small gamma, Greek, tilde]. The ODT shifts with k̃, with lower k̃ pushing the transition to higher [small gamma, Greek, tilde]. Consistent with previous studies, large size-dispersity (low k̃) suppresses long-range order. Analysis of defect formation in athermal systems through T1 transitions shows that defects form more easily at low k̃, with the activation barrier increasing with k̃. This is supported by Voronoi entropy and defect concentration measurements, which are higher for low-k̃ systems, reflecting a greater defect content and providing a mechanistic explanation for the observed ODT trends.

We also studied the phonon modes of our system far from the ODT and found that it behaves like a harmonic spring network with pairwise interactions. The vibrational density of states exhibits two van Hove peaks, in agreement with theoretical predictions. The similarity between the vibrational density of states of the Voronoi model and the harmonic spring network suggests that the coupled interactions in the model effectively become pairwise far from the ODT. Furthermore, deviations in the orientational order parameter from that of the hexagonal crystal lattice in this regime, where the defect concentration is negligible, can be well captured by constructing displacement fields from normal mode analysis.

We believe that our model provides a useful conceptual framework for investigating soft-matter systems in which effective particle interactions can be approximated using the energy functional employed in this study. Representative experimental realizations could include emulsions, microgels, or micellar assemblies, where particle deformability plays a central role. In practice, the material parameters such as line tension and resistance to mass exchange between particles along with average particle size would need to be tuned to map onto the dimensionless parameters quantified here. The energy functional considered is also simplified, and realizing an exact quantitative experimental analog could be a limitation of the present simulations. Further deviations between model predictions and experimental observations may arise from non-space-filling particle configurations and from departures from the idealized two-dimensional geometry assumed in our simulations.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI) is available that includes system size analysis, the influence of initial conditions, and bond-orientational correlations across the order–disorder transition, along with other supporting results. See DOI: https://doi.org/10.1039/d5sm01097g.

Relevant raw data are available on Zenodo. See DOI: https://doi.org/10.5281/zenodo.18133950.

Acknowledgements

We gratefully thank the support from the National Science Foundation (Grant No. NSF 2230946).

References

1. Y. Zhou, S. L. Bore, A. R. Tao, F. Paesani and G. Arya, npj Comput. Mater., 2023, 9, 224.
2. D. Chintha, S. K. Veesam, E. Boattini, L. Filion and S. N. Punnathanam, J. Chem. Phys., 2021, 155, 244901.
3. G. Carrot, A. E. Harrak, J. Oberdisse, J. Jestin and F. Boué, Soft Matter, 2006, 2, 1043–1047.
4. A. Nabiyan, A. Muttathukattil, F. Tomazic, D. Pretzel, U. S. Schubert, M. Engel and F. H. Schacher, ACS Nano, 2023, 17, 21216–21226.
5. F. Zeng and S. C. Zimmerman, Chem. Rev., 1997, 97, 1681–1712.
6. Z. Lyu, L. Ding, A. Tintaru and L. Peng, Acc. Chem. Res., 2020, 53, 2936–2949.
7. O. Altintas, G. Hizal and U. Tunca, J. Polym. Sci., Part A:Polym. Chem., 2006, 44, 5699–5707.
8. K. Hayashida, T. Dotera, A. Takano and Y. Matsushita, Phys. Rev. Lett., 2007, 98, 195502.
9. J. Buitenhuis and S. Förster, J. Chem. Phys., 1997, 107, 262–272.
10. R. A. Segalman, A. Hexemer, R. C. Hayward and E. J. Kramer, Macromolecules, 2003, 36, 3272–3288.
11. I. W. Hamley, C. Daniel, W. Mingvanish, S.-M. Mai, C. Booth, L. Messe and A. J. Ryan, Langmuir, 2000, 16, 2508–2514.
12. V. Vitkova, M. Mader and T. Podgorski, Europhys. Lett., 2004, 68, 398.
13. E. L. Romero and M. J. Morilla, Int. J. Nanomed., 2013, 8, 3171–3186.
14. Y. Peng, Z. Wang, A. M. Alsayed, A. G. Yodh and Y. Han, Phys. Rev. Lett., 2010, 104, 205703.
15. M. R. Islam, R. Nguyen and L. A. Lyon, Macromol. Rapid Commun., 2021, 42, 2100372.
16. H.-Y. Chang, Y.-J. Sheng and H.-K. Tsao, J. Mol. Liq., 2022, 364, 120025.
17. Y.-H. Tsao, S.-K. Li, H.-K. Tsao and Y.-J. Sheng, Colloids Surf., A, 2024, 680, 132656.
18. I. M. Van Meerbeek, B. C. Mac Murray, J. W. Kim, S. S. Robinson, P. X. Zou, M. N. Silberstein and R. F. Shepherd, Adv. Mater., 2016, 28, 2801–2806.
19. D. L. Barton, S. Henkes, C. J. Weijer and R. Sknepnek, PLoS Comput. Biol., 2017, 13, e1005569.
20. D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, Soft Matter, 2014, 10, 1885–1890.
21. A. Scotti, M. Brugnoni, A. A. Rudov, J. E. Houston, I. I. Potemkin and W. Richtering, J. Chem. Phys., 2018, 148, 174903.
22. Y. Jin and H. A. Makse, Phys. A, 2010, 389, 5362–5379.
23. A. Donev, S. Torquato, F. H. Stillinger and R. Connelly, J. Appl. Phys., 2004, 95, 989–999.
24. S.-C. Mau and D. A. Huse, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 1999, 59, 4396–4401.
25. L. Sun, Z. Li, Y. Zhang, Y. Lu and S. Zhang, Prog. Mater. Sci., 2026, 155, 101531.
26. S. Hou, L. Bai, D. Lu and H. Duan, Acc. Chem. Res., 2023, 56, 740–751.
27. B. I. Halperin and D. R. Nelson, Phys. Rev. Lett., 1978, 41, 121–124.
28. S. C. Kapfer and W. Krauth, Phys. Rev. Lett., 2015, 114, 035702.
29. R. Guo, J. Li and B. Ai, Phys. A, 2023, 623, 128833.
30. E. P. Bernard and W. Krauth, Phys. Rev. Lett., 2011, 107, 155704.
31. P. Sampedro Ruiz and R. Ni, J. Chem. Phys., 2020, 153, 174501.
32. P. Sampedro Ruiz, Q. Lei and R. Ni, Commun. Phys., 2019, 2, 70.
33. S.-H. Choi, F. S. Bates and T. P. Lodge, Macromolecules, 2011, 44, 3594–3604.
34. I. Bouhid de Aguiar, T. van de Laar, M. Meireles, A. Bouchoux, J. Sprakel and K. Schroën, Sci. Rep., 2017, 7, 10223.
35. E. Boattini, N. Bezem, S. N. Punnathanam, F. Smallenburg and L. Filion, J. Chem. Phys., 2020, 153, 064902.
36. A. Boromand, A. Signoriello, F. Ye, C. S. O’Hern and M. D. Shattuck, Phys. Rev. Lett., 2018, 121(24), 248003.
37. M. A. Spencer, Z. Jabeen and D. K. Lubensky, Eur. Phys. J. E:Soft Matter Biol. Phys., 2017, 40, 2.
38. D. M. Sussman, Comput. Phys. Commun., 2017, 219, 400–406.
39. I. Tah, D. Haertter, J. M. Crawford, D. P. Kiehart, C. F. Schmidt and A. J. Liu, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2322732121.
40. Y.-W. Li, L. L. Y. Wei, M. Paoluzzi and M. P. Ciamarra, Phys. Rev. E, 2021, 103, 022607.
41. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
42. K. Brakke, R. Phelan and D. Weaire, Exp. Math., 1995, 4, 181–192.
43. T. O. Bello, S. Lee and P. T. Underhill, Soft Matter, 2022, 18, 5106–5113.
44. delaunayTriangulation – Delaunay triangulation in 2-D and 3-D – MATLAB, https://www.mathworks.com/help/matlab/ref/delaunaytriangulation.html, (accessed 29 September 2025).
45. T. C. Hales, Discrete Comput. Geom., 2001, 25, 1–22.
46. M. Durand and J. Heu, Phys. Rev. Lett., 2019, 123, 188001.
47. S. Prestipino, F. Saija and P. V. Giaquinta, Phys. Rev. Lett., 2011, 106, 235701.
48. N. Gribova, A. Arnold, T. Schilling and C. Holm, J. Chem. Phys., 2011, 135, 054514.
49. H. J. Schöpe, G. Bryant and W. van Megen, Phys. Rev. E, 2006, 74(6), 060401.
50. M. Cerbelaud, F. Mortier, H. Semaan, J. Gerhards, B. Crespin, R. Ferrando and A. Videcoq, Mater. Today Commun., 2024, 38, 107973.
51. M. Ozawa, G. Parisi and L. Berthier, J. Chem. Phys., 2018, 149, 154501.
52. M. Ozawa and L. Berthier, J. Chem. Phys., 2017, 146, 014502.
53. E. Bormashenko, M. Frenkel, A. Vilk, I. Legchenkova, A. A. Fedorets, N. E. Aktaev, L. A. Dombrovsky and M. Nosonovsky, Entropy, 2018, 20, 956.
54. D. Paulovics, E. Bormashenko, C. Raufaste and F. Celestini, Phys. Rev. E, 2024, 110, 024302.
55. E. Bormashenko, I. Legchenkova and M. Frenkel, Entropy, 2019, 21, 452.
56. A. Pertsinidis and X. S. Ling, Phys. Rev. Lett., 2001, 87, 098303.
57. J. M. Kosterlitz and D. J. Thouless, J. Phys. C-Solid State Phys., 1973, 6, 1181.
58. A. P. Young, Phys. Rev. B:Condens. Matter Mater. Phys., 1979, 19, 1855–1866.
59. D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, Soft Matter, 2014, 10, 1885–1890.
60. A. Ghosh, R. Mari, V. K. Chikkadi, P. Schall, A. C. Maggs and D. Bonn, Phys. A, 2011, 390, 3061–3068.
61. L. Van Hove, Phys. Rev., 1953, 89, 1189–1193.

---