Melting of a two-dimensional monodisperse cluster crystal to a cluster liquid
Wenlong Wang, Rogelio Díaz-Méndez, Mats Wallin, Jack Lidmar, Egor Babaev
MMelting of a two-dimensional monodisperse cluster crystal to a cluster liquid
Wenlong Wang, ∗ Rogelio D´ıaz-M´endez, Mats Wallin, Jack Lidmar, and Egor Babaev Department of Physics, Royal Institute of Technology, Stockholm, SE-106 91, Sweden (Dated: April 26, 2019)Monodisperse ensembles of particles that have cluster crystalline phases at low temperatures can model anumber of physical systems, such as vortices in type-1.5 superconductors, colloidal suspensions and cold atoms.In this work we study a two-dimensional cluster-forming particle system interacting via an ultrasoft potential.We present a simple mean-field characterization of the cluster-crystal ground state, corroborating with MonteCarlo simulations for a wide range of densities. The efficiency of several Monte Carlo algorithms are comparedand the challenges of thermal equilibrium sampling are identified. We demonstrate that the liquid to cluster-crystal phase transition is of first order and occurs in a single step, and the liquid phase is a cluster liquid.
I. INTRODUCTION
Melting of two-dimensional crystals is a fundamen-tally interesting problem that attracted interest for decades.Within the celebrated Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory [1, 2], the melting is a two-stageprocess occurring with an intermediate hexatic phase of quasi-long range orientational order separating the crystal from thefluid. In this scenario, the two phase transitions (crystal-to-hexatic and hexatic-to-fluid) are predicted to be of second or-der, driven by the unbinding of different types of topologi-cal defects (dislocations and disclinations) [2, 3]. However,after many decades of experimental and numerical research,evidence was collected that several 2D systems fall outsideKTHNY theory, melting in a single first order transition, orin a two-step process by one first order and one continuoustransitions [4–6]. While the characterization of melting tran-sitions of 2D crystals is a long-standing problem that has re-ceived great attention, almost all the effort has been focusedon simple potentials forming simple crystals, with only a fewnumerical studies available on the generalized exponentialmodel [7, 8].Monodisperse systems, in which particles are all of thesame type, can nonetheless exhibit very complex hierarchi-cal structure formations such as lattices of clusters and exoticphases such as glasses [9–15]. The wide array of behaviorsincluding rich self-organizing patterns, dynamics and the ther-modynamic equilibrium phases in such systems are far fromunderstood [16–19]. Particles with cluster-forming interac-tions emerge in various contexts ranging from soft matter,complex molecules, to cold atoms and vortices in supercon-ductors [20, 21]. The problem is also relevant for the physicsof multi-component superconductors where a rich variety ofvortex cluster solutions have been found [13, 22, 23].For non-cluster-forming systems, the ordered phase con-sists in a normal crystal with a single particle located at eachlattice site. On the other hand, for cluster-forming potentialsthe particles arrange in a cluster crystal, in which each latticesite is occupied by a cluster of more than one particle. Un-der the above conditions, there is a simple criterion by Likos ∗ Electronic address: [email protected] et al. relating the cluster-forming ability of the interactionsto the sign of its Fourier transform [24, 25]. Namely, if theFourier transform of the interaction potential is positive def-inite, the system is non-cluster-forming, and cluster-formingotherwise.Since the ordered states of cluster-forming particles in prin-ciple allow more types of defects, this raises the questionof the nature of melting transition of cluster-crystals in 2D,and to what extent the KTHNY theory applies. In this work,we focus on a cluster-forming ultrasoft potential that has at-tracted considerable attention as a model for colloidal suspen-sions and ultracold atoms [13, 20]. Recently a similar kindof cluster-forming interaction potential was found to arise be-tween the vortices in type-1.5 superconductors. In the bulk,intervortex forces are long-range attractive and short-range re-pulsive; see e.g. Refs. [26–29]. For thin films an additionalrepulsive interaction arises due to magnetic stray fields thatgives raise to vortex cluster crystals [13, 23]. We present asimple mean-field analysis of the ordered cluster-crystallineground state. In particular, our analysis captures a number ofinteresting properties such as the existence of an onset den-sity for the cluster formation, the density dependence of thelattice constant, and also that the triangular lattice is more sta-ble than the square lattice at all densities. We corroborate ourresults with a number of Monte Carlo simulations, finding aqualitative and also reasonable quantitative agreement.In order to ensure the proper equilibration of the system andalso for comparing the efficiency of different Monte Carlo al-gorithms, we here go beyond temperature annealing schemes,exploring a number of different Monte Carlo methods. Wealso use these methods to identify the most challenging fac-tors involved in the thermal equilibration. The algorithms wehave implemented are widely used in spin glasses [30–32],where systems are highly disordered and frustrated, and equi-libration is usually essential for progress. They include simu-lated annealing [33] for optimizations, and parallel tempering[34–36] for equilibrium sampling. One more recent but lessknown algorithm is the population annealing [37–44]. Bothpopulation annealing and parallel tempering are extended en-semble methods with similar efficiency for spin glasses. Wefind, however, that parallel tempering can be more efficientthan population annealing for our less frustrated model. Usingour equilibrium sampling to generate energy histograms for aseries of system sizes, we find a single first-order phase transi- a r X i v : . [ c ond - m a t . s o f t ] A p r tion, from cluster liquid to cluster crystal, featuring a double-peak structure in the histograms. In addition we studied thecase of particles moving in a harmonic trap. Also here weobtain a cluster crystal phase, notably with a lattice spacingapproximately independent of the distance to the trap centerand similar to the unconfined case. However, the occupationnumber of the clusters decreases with the distance to the cen-ter. Accordingly, the melting temperature becomes positiondependent with more stable highly populated clusters locatedat the center of the trap, generating an interesting inhomoge-neous melting.The paper is structured as follows. First we introduce ourmodel, observables, simulation methods, and the mean-fieldanalysis in Sec. II. The algorithms and numerical results aredescribed in Sec. III. Finally, a summary of the main findingsand discussion of their implications is presented in Sec. IV. II. MODEL, OBSERVABLES AND METHODSA. Model and observables
The two-dimensional system of monodisperse particles in-teracting via an ultrasoft potential is described by the Hamil-tonian H = (cid:88) i We have studied the system using three different MonteCarlo methods: simulated annealing (SA), population anneal-ing (PA), and parallel tempering (PT). The major and com-mon component for each algorithm is a Monte Carlo sweep.We use the Metropolis algorithm in this work. Each MonteCarlo sweep is a sequential update of all the N particles.For each update, we propose to shift a particle randomlywithin a square box of length (cid:112) /n centered on the particle.We assume readers are familiar with the single-temperatureMarkov chain, and the SA for sequential cooling through a se-ries of single-temperature Monte Carlo following an anneal-ing schedule. We additionally use PT for multiple Markovchains with exchange of replicas between different tempera-tures. The PA method is, however, relatively new, and wediscuss this algorithm in some detail here.Like PT, PA is an extended-ensemble algorithm for thermalequilibrium sampling. It is also similar to SA but with a largepopulation of R replicas. It is often the case, also in this work,that replicas are initialized randomly at β = 0 and cooled toa target temperature T min = 1 /β max . The replicas traverse anannealing schedule with N T temperatures by slowly lower-ing the temperature. The population is, however, resampled ateach temperature step to maintain thermal equilibrium, with aself-consistent reweighting process. When the temperature islowered from β to β (cid:48) , the population is resampled. The meannumber of copies of replica i with energy E i is proportionalto the appropriate reweighting factor, exp[ − ( β (cid:48) − β ) E i ] . Theconstant of proportionality is chosen such that the expectationvalue of the population size at the new temperature is R . Inpractice, the number of each replica is rounded to the flooror ceiling of the expectation value with the right probability,called nearest integer resampling. The resampling is followedby N S sweeps to each replica of the new population, and thecycle of resampling and sweeps repeats until the target tem-perature is reached.Note that PA and SA have a similar structure. Turning offresampling, PA becomes SA of R independent runs. PA hasbeen found to be similar in efficiency to PT for spin glasses,but has a number of additional useful features such as havingintrinsic equilibration measures (see below), giving direct ac-cess to free energy using the free energy perturbation method,and is massively parallel.The equilibration measure of PA in this work is based onthe family entropy S f and the entropic family size ρ s [40, 41].These quantities are, respectively, defined as S f = − (cid:88) i ν i ln ν i , (4) ρ s = lim R →∞ R/e S f , (5)where ν i is the fraction of the population that has descendedfrom replica i in the initial population. Intuitively, exp( S f ) characterizes the number of effective surviving families or thediversity of the population, and ρ s characterizes the averageeffective surviving family size. A large S f ensures accuratesampling of the equilibrium distribution. On average, S f de-creases with β with a rate that depends on the free energy land-scape and the simulation parameters. We use S f ≥ ln(100) in this work, which has been used in many systems [40, 41].If this condition is not fulfilled at the final temperature, thesimulation is restarted with a larger initial population. In ad-dition to providing an equilibrium measure, S f is also usefulfor identifying bottlenecks of the equilibration. C. Mean-field analysis of the cluster crystal ground state In this section, we develop a mean-field analysis to capturequantitative features of the cluster crystal phase. For simplic-ity, we focus on the ground state which has no thermal fluctu-ations. In particular, we are interested in how the mean clus-ter size c s and lattice constant a change as a function of theparticle density n . It is interesting to see whether the mean-field treatment can catch the onset of particle clustering. Fi-nally, we also compare the triangular lattice and the squarelattice. We emphasize that comparing average energy per par-ticle for different lattices has been a well-established practice,see e.g. recent works of [15, 17]. Nevertheless, the applicationto cluster-forming particles provides some interesting insightson the clustering features. Let the area of the unit cell be A . Then n = c s /A by defi-nition, and A = √ a / for the triangular lattice and A = a for the square lattice. Note that the optimum c s and a arehence related for a given density n . We define a single parti-cle interaction energy for a particle sitting at the origin withoutloss of generality as (cid:15) ( c s ) = ( c s − U (0) + (cid:88) i (cid:54) =0 c s U ( r i ) , (6) = (cid:88) i c s U ( r i ) − U (0) , (7)where the summation is over lattice sites and r i denotes thedistance from the lattice site i to the origin or the lattice site .The first and second terms of Eq. (6) are the intra-cluster andinter-cluster interaction energies, respectively. Minimizationof (cid:15) with respect to c s gives the optimum cluster size c s andlattice constant a . In summary, for a given density n , we aimto find the global minimum of (cid:15) ( c s ) . For each c s , we firstcompute the lattice constant a and then numerically sum thesingle particle interaction energy (cid:15) . We run this procedure forboth the triangular lattice and the square lattice.The optimum single particle interaction energy, lattice con-stant, and mean cluster size for both lattices are shown inFig. 1. Firstly, the single particle interaction energy for thetriangular lattice is always smaller than the square lattice, andthe difference increases with increasing density. Both con-verge to zero in the zero-density limit. Second, the mean clus-ter size has the following onset density of particle clusteringat n c ≈ . . When the density is lower than n c , clusteringis not favored. When n ≥ n c , clustering occurs in the groundstate. Once clustering occurs, the lattice spacing remarkablysettles down and becomes independent of the density. As aresult, the mean cluster size grows linearly with density in theclustering regime. This is compatible with the numerical re-sults of Ref. [19]. The single particle interaction energy ishence also asymptotically linear with density from Eq. (7). III. NUMERICAL RESULTS This section has two main parts. We first present the com-parison of different MC algorithms, and the bottlenecks forthermal sampling. Then we present results for equilibriumproperties of the system, including comparison with the mean-field results, the order of the phase transition, and cluster crys-tals in a parabolic trap. A. Algorithms We start by showing that procedures based on simulatedannealing cannot easily maintain full thermal equilibrationthroughout all temperatures even for the clean system. Itmight be tempting to expect otherwise, especially consider-ing that after a deep quench from a random configuration, thesystem can reasonably restore the order parameter close to thetransition temperature [13, 46]. When the system is cooled Q u a n t i t i e s n TL a TL (cid:15) TL c s SL a SL (cid:15) SL c s FIG. 1: The cluster crystal lattice constant a , single particle inter-action energy (cid:15) , and mean cluster size c s as a function of density n for both the triangular lattice (TL) and the square lattice (SL) underthe mean-field approximation. Note that the lattice constant is essen-tially independent of the density once the clustering occurs ( c s > )at about n ≈ . . Therefore, the mean cluster size grows linearly with n in the clustering regime. The single particle interaction en-ergy for the triangular lattice is always lower than the square latticeat any finite density, and is also asymptotically linear with density. slowly, one might find the crystal lattice easily, while furtherlowering the temperature would only suppress lattice fluctu-ations. However, running a simulated annealing, the systemdoes not always find a perfect lattice, but may contain defectssuch as grain boundaries near the transition temperature T C .Here, T C can be estimated as the temperature where φ de-parts from zero or from the peak of the specific heat. Oncethey have formed, it may be difficult to remove such defectsusing local updates. See Fig. 2 for an example of a reasonablycareful annealing of 1000 independent runs for N = 1000 particles at density n = 1 . . While some fraction of theruns find pretty good lattice structures, the majority of themdo not, and some fail badly with little long-range order. Themean, shown by the black curve, is systematically below thetop-most curve, which runs toward one when T = 0 . SimpleMonte Carlo running at a single temperature, as in quenchingdynamics for glass formations, would be only worse.Nevertheless, comparing the onset of the ordered phasewith the peak of the specific heat in Fig. 3, it seems reason-able that SA captures the transition temperature T C reason-ably well. Furthermore, SA remains a valuable tool for opti-mization to reach very low temperatures (including T = 0 ) asthermal equilibration down to very low temperatures could berather difficult even for small sizes, as we will discuss. SA hasbeen used previously to study phase diagrams for other clustercrystal forming potentials [47, 48].Population annealing and parallel tempering, unlike simu-lated annealing, are both designed for thermal sampling. Wefirst use the family entropy of population annealing to iden-tify bottlenecks in the simulations. As might be intuitivelyexpected, we find that the phase transition (which we will φ T Simulated annealing φ T FIG. 2: Evolution of the order parameter φ upon slow cooling usingsimulated annealing. Each thin line is an independent run and thethick black curve is the mean of these runs which is also highlightedfor clarity in the inset panel. Only 10 temperatures are shown forthe mean, although there are 1000 temperature steps. Note that notall runs go toward φ = 1 , showing that the ground state is a globalminimum but not a global attractor. These are 1000 independentruns for N = 1000 at density n = 1 . from β = 0 to T = 0 .We cool the system quite slowly with a total of 1000 steps, with100 sweeps at each temperature. Half of the temperatures are linearin β up to β = 10 , and the other half are linear in T to T = 0 .This demonstrates SA can be used for optimization, but is not fullyreliable for thermal sampling. demonstrate is first order in a later section) is a generic bot-tleneck. This is reflected as a sharp decrease of family en-tropy near the transition. Clearly, increasing the number ofsweeps near the transition can improve equilibration signifi-cantly. See Fig. 3 for an example with n = 1 . (solid lines),where we have used an annealing schedule with 100 temper-atures linear in β up to T C and 100 ones linear in T in thelow temperature part. Notice that if we spent 100 sweeps perreplica near the phase transition ( β ∈ [8 , ), the family en-tropy drops much less compared with 40 at all temperatures.Furthermore, this is even more efficient than simply doublingthe population size keeping N S = 40 , which involves moretotal work. Therefore, we conclude that the phase transitionrequires more work (either more sweeps or equivalently moretemperatures) for efficient thermal sampling.Note also that from the family entropy of n = 1 . (dashedlines) in Fig. 3 the phase transition may not be the only bot-tleneck in thermal sampling, which is accompanied by ananomalous peak in the specific heat at low temperatures. Forthis harder case, we have used an annealing schedule with 100temperatures linear in β up to T C and 200 ones linear in T inthe low temperature part. We have again spent 100 sweepsper replica near the phase transition ( β ∈ [8 , ) comparedwith 40 at other temperatures such that the transition is not asignificant bottleneck. In this example, the anomalous peak isthe major bottleneck. While this is not a generic peak for dif-ferent sizes and densities, it occasionally shows up, and sincetemperatures are very low, spending more sweeps around thispeak does not help much in contrast to the transition bottle-neck. By careful inspection, we find that the nature of thebump is due to a subtle lattice reorganization: not in the lat-tice form but in the number of clusters. There are 86 clus-ters at higher temperatures before the bump but 90 clusters atlower temperatures after the bump, and between the two thepopulation has a mixture of replicas of both kinds. Note thatboth lattices are the same triangular crystals, and there is noglassy states here as the replicas are kept in thermal equilib-rium. This is clearly a balance of energy and entropy of thetwo finite lattices. While we believe this is probably a finite-size effect, fitting a lattice in a finite space, it certainly sug-gests that equilibration below and much below T C for finitesystems may involve subtle differences. For example, if oneis running SA for a limited number of independent runs, thetrue ground state might be missed. C V / N ( N = ) S f β n = 1 . n = 1 . N S = 40 , R = 10 N S = 40 , , R = 10 N S = 40 , R = 2 × N S = 40 , , R = 2 × (a)(b)FIG. 3: Specific heat [(a)] and family entropy [(b)] as a func-tion of β in two representative scenarios, with and without the low-temperature bump. In both cases, the first-order phase transition is ageneric bottleneck, which requires more work near the transition foran efficient sampling. For example, spending 100 sweeps near thetransition (40 sweeps at other temperatures) is even more efficientthan doubling the system size keeping 40 sweeps at all temperatures.See text for more details. However, there could be an additionalbottleneck associated with the anomalous peak at low temperatures,such as for n = 1 . . Careful inspection shows that this is due tolattice reorganization: there are 86 clusters at higher temperaturesbelow the bump and 90 above. Such bumps at very low temperaturespose even greater equilibration challenges, as shown by the sharp de-crease of the family entropy. Furthermore, spending more sweepshere does not help much, in contrast to the phase transition. Theshoulder peaks at small β are not generic but occur only for smallsizes. We use a small size here as we traverse a wide temperaturerange to β = 60 which is hard for large sizes. We now focus on the efficiency of population annealing andparallel tempering. Despite the two algorithms are similarlyefficient in spin glasses that are highly disordered, it appearsthat population annealing is not as efficient as parallel temper-ing for the clean system. We have also found similar results in our recent simulations of the three-dimensional Ising model.Suppose one is interested in the phase transition regime orbelow, we think the following factor is limiting the power ofpopulation annealing for clean systems: population annealingis a sequential Monte Carlo and requires to start the popu-lation at a high temperature. But resampling is only mean-ingful when the population size R is sufficiently large, hencepopulation annealing necessarily spends lots of efforts in thehigh temperature regime before reaching the transition tem-perature. No initial equilibration is needed if the simulationstarts from β = 0 , but it takes a number of steps to get closeto T C . One can also choose to start at a high but finite tem-perature directly, but this would require population annealingto equilibrate initially a large number of replicas. In contrast,parallel tempering can equilibrate fast for clean systems (shortthermalization times and correlation times) directly near thetransition for a much smaller number of replicas, and there-fore good or reasonable statics can be obtained more quicklyfor a fixed computational effort. In summary, we think the se-quential preparation is slowing down PA, compared with PT.Our results are, however, not in disagreement with the sim-ilar efficiency in disordered systems. In that scenario, boththe thermalization times and correlation times of PT are muchlonger, driving the efficiencies of PA and PT together. Never-theless, PA has been very useful in our studies, e.g., in iden-tifying simulation bottlenecks. It provides both sharper andmore quantitative information than the exchange probabilitiesof parallel tempering. For equilibrium sampling of monodis-perse particles in disordered potentials where equilibration isslow, we expect that the massively parallel PA would againbecome relevant.Finally, there is an intrinsic complexity in the model itself:the long-range interactions and the associated O ( N ) com-plexity for one sweep. While it is an option to introduce a cut-off to the potential, the efficiency necessarily decreases withincreasing density and is less relevant when the system size isnot exceptionally large. Therefore, this is a significant limit-ing factor on the sizes accessible to equilibrium simulations.We have managed to fully equilibrate about − parti-cles around T C in this work, using state-of-the-art computingresources. In contrast, for clean systems with short-range in-teractions such as the Ising ferromagnets, one can reach ordersof magnitude larger sizes.In summary, the long-range interactions and the phase tran-sition are generic bottlenecks, and extra bumps at low temper-atures when present are often a significant bottleneck if oneis interested in the low-temperature regime. Increasing thenumber of Monte Carlo sweeps near the phase transition usu-ally improves efficiency. We find that SA is not sufficient formaintaining full equilibration throughout the annealing, andparallel tempering is more efficient than population anneal-ing for clean systems, although for disordered and frustratedsystems the efficiency becomes similar. TABLE I: Simulation parameters of some representative examplesof the Monte Carlo methods. Here n is the density of particles, N isthe number of particles, R is the number of replicas, N T is the num-ber of temperatures, N S is the number of sweeps for each replica ateach temperature, β min and β max are the minimum and maximuminverse temperatures, respectively, and M is the number of indepen-dent runs.Method n N R N T N S β min β max M SA 1.1 1000 ∞ 201 40 0 60 6PA 1.25 200 × 301 40 0 60 2PA 1.5 200 201 60 0 60 6PA 1.75 200 201 40 0 60 2PA 2 200 201 60 0 60 6PA 2.5 200 201 40 0 60 2PA 3 200 201 40 0 40 2PT 1.1 200 − . × − . × − . × − . × 11 11.5 2PT 1 1000 − . × 10 16 2PT 1.5 1000 − . × − . × − . × − . × B. Equilibrium properties We first check the mean-field results and present the phasediagram, followed by a demonstration that the liquid to clustercrystal phase transition is first order, using energy histograms.Our results for thermal sampling are studied mainly using PT,but some PA results are also presented when relevant. Wefind ground states using simulated annealing in the same wayas discussed in the previous section, down to T = 0 with atleast 1000 independent runs. The simulation parameters aresummarized in Table I.The T C is estimated from the peak of the specific heat. Weextract from the ground state the lattice constant and the meancluster size to compare with the predictions of the mean-fieldtheory. Ideally we should find well ordered ground states with φ ≥ . . While this is readily achievable for large densities,it turns out to be quite a difficult task for low densities wherethe transition temperatures are very low and particles interactonly very weakly. We therefore require this criterion only fordensities of n ≥ . Nevertheless, we find that the states stillhave reasonably good local order, despite of the poor globalorder. Remarkably we find that the location of the first non-trivial peak of the correlation function does not significantlydepend on whether the system is in the true ground state orjust a low-energy state, and often even a glassy state.Our results are summarized in the top panel of Fig. 4. Wefind the mean-field analysis is in reasonably good agreementwith the numeric results. The system has an onset density, and Q u a n t i t i e s n MFT a MC a MFT c s MC c s (a) T C n PAMC N = 200 C V PTMC N = 1000 C V SAMD N = 600 − φ (b)FIG. 4: Top panel: Comparison of the mean-field and simulated an-nealing results characterizing the ground states [(a)]. Bottom panel:Phase diagram estimated from the specific heat peaks of the equi-librium MC simulations, and the estimation of simulated annealingusing molecular dynamics and φ [13] [(b)]. Errorbars of MC resultsare smaller than the symbols. the lattice constant is independent of density in the cluster-crystal region. Below the onset density, the MFT is triviallyexact as c s = 1 . The MFT, however, does not find the criticalonset density exactly. The numerical onset density appearsto be between . and . , while the MFT predicts . . Theasymptotic lattice constant a = 1 . found in the simulationsat the two largest densities shows a very good agreement withthe MFT result a = 1 . . The small discrepancy is dueto neglecting of correlations in our MFT approximation. Forexample, clusters of different sizes are likely not randomlydistributed [13]. There is also place for finite-size effects onthe Monte Carlo side, as we cannot exclude that for large sizesthe MC results will get closer or even converge to the MFTprediction. Nevertheless, the results are reasonably close, andfrom our MFT analysis, we can clearly see that the onset ofparticle clustering is a process of minimization of the latticeenergy at low temperatures.The phase diagram is summarized in the bottom panel ofFig. 4. We have identified T C as the location of the maximumpeak of the specific heat. This again works well for high den-sities n (cid:38) . There are small finite-size effects as the peaks of N = 200 and 1000 are not at the same temperature, but theyare sufficiently close to get an estimate of the phase bound-aries. Our results regarding the melting temperature are ingood agreement with those obtained from the evolution of φ using simulated annealing [13, 19].A region of phase coexistence can be present at a first or-der melting of an ordinary particle crystal, where the densityjumps discontinuously at the transition when simulated in anNVT ensemble. For a cluster crystal, on the other hand, thelattice spacing can readjust without changing the total densityby changing the cluster occupation numbers, and for this rea-son it is not clear that a phase coexistence region necessarilyexists. In fact, our simulation results do not clearly indicate acoexistence regime in the phase diagram, suggesting that sucha region might be very narrow or absent in the model investi-gated here. c s β n = 1 . n = 1 . FIG. 5: The corresponding mean cluster size c s as a function oftemperature for Fig. 3. The density n = 1 . shows clearly a changeof c s at the anomalous peak. Other than this finite size effect, thecluster size does not depend on the temperature in the crystal phasefor this ultrasoft model. The errorbar here is for the spread of thedistribution with each point estimated from configurations, notthe errorbar of the mean of c s which is smaller by a factor of √ 100 =10 . Note, however, the clustering procedure or counting of clusters ismost reliable only at low temperatures. See the text for more details. It is interesting to know to what extent the lattice parame-ters apply to finite temperatures. It is clear that the identityand hence counting of clusters is less precise at finite tem-peratures, due to particle hopping among clusters. However,we can still look for a trend and study if there is an obvi-ous change in the size of clusters, following our clusteringmethod. Such results are particularly relevant in the low-temperature limit. Therefore, we go back to the examplesin Fig. 3 and the temperature evolution of the mean clus-ter size c s is shown in Fig. 5. We see clearly the fluctua-tions of the mean cluster size from configuration to config-uration is essentially zero in the T → limit, and other thanthe finite-size effect at the anomalous peak, the mean cluster size approximately does not depend on the temperature in thelow temperature limit. This even holds close to the transi-tion temperature, though one should be cautious on the defini-tion and counting of clusters in this regime. Nevertheless, thisclearly shows the ground state lattice constant remains a rel-evant length scale throughout the solid phase. We show thatthis is actually true even beyond the phase transition in theliquid phase, using pair correlation functions in Sec. III C.The order of this phase transition between the fluid and thecluster-crystal region has not been conclusively investigatedpreviously. It has been reported that this transition occurs ina single step for the generalized exponential model of index 4(GEM4), where the transition was encountered to be of first-order nature [7]. Using density functional theories and simu-lations a transition in the ordering of clusters but not for parti-cles was found. Here we use equilibrium sampling to generateenergy histograms of our model (similar to GEM4) to directlydemonstrate the discontinuous nature of the transition. Ourconclusion is therefore in agreement with Ref. [7].The energy histograms for several sizes near T C at density n = 1 . are shown in Fig. 6. We have refined our temperatureset to find symmetric distributions, which require a very finetemperature spacing. While there is still a bit of asymmetry,a clear trend can be observed in the energy histograms, wherea double-peak structure appears with a growing free energybarrier for increasing system size. This is a strong signatureof a first-order phase transition. In addition, we have mea-sured the maximum value of C V /N , shown in Fig. 7. We seethat after an initial transient at small sizes, the function devel-ops an approximately linear behavior, compatible with a first-order phase transition. Again, C V /N does not include anycluster parameter, so we conclude the liquid to cluster-crystalphase transition is first order, in contrast to the KTHNY the-ory. While there is no theoretical basis to exclude the possi-bility of two transitions (either two second order or two firstorder, or one first order and one second order), we do not seeany evidence that this occurs in our simulations of the ultrasoftpotential. C. Cluster liquids The clustering of particles in the cluster crystals remark-ably persists to the liquid phase, and hence we here proposethe concept of a cluster liquids associated with melted clustercrystals. It is relatively straightforward to recognize clusterliquids from some typical configurations shown in Fig. 8. De-spite this clearness in our present case, a generic definitionturns out to be difficult due to the rich phase diagrams of var-ious potentials. Here, we discuss this only qualitatively andrestrict our discussions to simple cases where the associatedcrystal phases are simple crystals like the triangular lattice,excluding more exotic ones such as stripes. The counterpartof corresponding concepts such as stripe liquids can be real-ized.To have a well defined cluster liquid, the liquid phaseshould preserve some of the cluster order of the solid phase.The configurations in Fig. 8 motivate us to look at the pair cor- H i s t og r a m ( a ) U = 1 ( c ) M = 20 H i s t og r a m E/N ( b ) U = 4 E/N ( d ) M = 40 N = 200 N = 400 N = 800 N = 1200 (c)(a) (b)(d) H i s t og r a m ( a ) U = 1 ( c ) M = 20 H i s t og r a m E/N ( b ) U = 4 E/N ( d ) M = 40 β < β c β (cid:46) β c β (cid:38) β c β > β c (g)(e) (f)(h)FIG. 6: Top panel: Energy histogram at the phase transition for aseries of sizes N = 200 , , and , respectively [(a-d)] atdensity n = 1 . . There is a double-peak structure emerging withan increasing free energy barrier with increasing size N , suggestingthat the liquid to cluster-crystal phase transition is first order. Bot-tom panel: Energy histogram for four selected temperatures beforeand after the phase transition of N = 1200 . The selected temper-atures are β = 11 . , . , . , . [(e-h)], wherethe critical temperature β c = 11 . . relation function. If the particles are indeed clustering in theliquid phase, the first few relevant crystal peaks and in partic-ular the cluster lattice constant peak should appear in the cor-relation function due to the short-range order of liquids. Thisis indeed the case, as shown in Fig. 9 where the correlationfunctions for density n = 1 . are shown near the phase tran-sition. The second peak from left is the lattice constant peakwhich is almost temperature independent. This peak persistsdeep into the liquid phase but becomes less pronounced. Thisindicates that particles are clustering before solidification tocluster crystals, and clustering does not appear suddenly atthe phase transition. Since the lattice constant is independentof density for cluster crystals, we expect that this peak in theliquid phase is also density independent, which is the case asshown in Fig. 9. Therefore, the cluster liquid is well defined. C V m a x / N ( n = . ) N FIG. 7: Maximum of the specific heat per particle max( C V ) /N asa function of N at the density n = 1 . . We see after an initial tran-sient at small sizes, the function develops approximately into a linearfunction, compatible with a first-order phase transition. Errorbars aresmaller than the symbols.(a) (b) (c) (d)(e) (f) (g) (h)FIG. 8: Typical cluster liquid dynamics for the ultrasoft potential atdensity n = 1 . and β = 10 [(a-h)]. The clusters can oscillate anddiffuse for a relatively long time before a particle hops between clus-ters. In this example, five particles are highlighted. From t = 0 to t = 300 , they form two clusters. At about t = 400 , the blue clusterloses a particle and forms a single-particle cluster. This cluster lateron at about t = 600 joins the green cluster. The times are in units ofMonte Carlo sweeps. See the relevant movie and comparison with aparticle liquid in Ref. [49]. Having defined cluster liquids, we discuss the major dif-ferences of such liquids and particle liquids in terms of theirtypical dynamical processes. In particle liquids, the domi-nant processes are particle oscillations in the potential formedby the neighbors and the eventual diffusion of the particles.These two processes also occur in cluster liquids for the cen-ters of mass of clusters. However, due to the clustering, the g ( r ) × g ( r ) × r T > T C ( n = 1 . T (cid:38) T C ( n = 1 . T = T C ( n = 1 . T (cid:46) T C ( n = 1 . T (cid:38) T C ( n = 1) T (cid:38) T C ( n = 1 . T (cid:38) T C ( n = 1 . T (cid:38) T C ( n = 2) (a)(b)FIG. 9: Top panel: The first few peaks of the pair correlation func-tion occur at the same locations for cluster liquids and cluster crys-tals [(a)]. Bottom panel: These peaks in cluster liquids are density-independent like the density-independent lattice constant of clustercrystals [(b)]. Note the second peak from left is the cluster latticeconstant peak, its height grows when T is lowered [(a)] or when n is increased at fixed N [(b)]. The figure demonstrates that clusteringoccurs already in the liquid phase and not at the phase transition. Inall cases, β = 4 , , . , , , , , , respectively. potential well is much stronger and remarkably particles in acluster can undergo many oscillations within the cluster diam-eter before the cluster occupation number changes. The dif-fusion of clusters is also much harder than single particles asthe diffusion of a cluster is now the center of mass of severalparticles. This might be related to our observation of a singlefirst-order phase transition. Note that many such oscillationsis necessary for well-defined cluster liquids. This is becauseparticles may also “cluster” together by chance for particleliquids. However, such “clusters” break almost immediatelyafter forming and as such they are different from the clustersin cluster liquids. In addition to the two processes above, clus-ter liquids have also particle exchanges among clusters or theconstant dynamical breaking and forming of clusters. Thereare two typical elementary processes, a cluster may emit orabsorb a particle, or two emitted single particles form a newcluster. The combination rate of the latter depends strongly onthe density of the system. For low densities such as n = 1 . , asingle particle may be relatively long-lived as a single-particlecluster. Nonetheless, such single particles are the most “reac-tive” in the liquid, having higher mobilities and strong ten-dencies to form clusters. For higher densities such as n = 2 ,an emitted particle joins another cluster much faster, whichappears effectively as a particle hopping among clusters.Next we present typical dynamics of a particle liquid and acluster liquid. We use here a Gaussian potential (GEM2) andthe ultrasoft potential for the two cases, respectively. We haveused N = 1000 particles at the same density n = 1 . and β = 10 for comparison, and an updating length scale for a par- ticle ± . for reasonably time scales of the dynamics. Notethat this length scale is still relatively small compared withthe lattice constant. In Fig. 8, cluster oscillations along with aparticle emission and a subsequent combination are shown forcluster liquids. For more details of both liquids we refer to amovie comparing the two cases [49]. In both cases, the timesare in the unit of sweeps. Finally, it should be emphasizedthat for MC dynamics to be realistic, our dynamics is mostlyrelevant to experimental settings where inertia effects are neg-ligible, i.e., the dynamics is heavily overdamped similar to themolecular dynamics simulations of Refs. [13, 19]. D. Confinement effects Finally, we have also studied the cluster-forming particlesin presence of a global parabolic trap, motivated by magnetictrap experiments of the atomic Bose-Einstein condensates [50,51]. More specifically, the particles interact with each otherand also with an external trap of the form V ( r ) = Ω r / . Wehave explored values from Ω = 0 . to , and obtain triangularcluster crystals for Ω (cid:38) . at low temperatures. Note thatdensity as well as periodic boundary conditions are no longerrelevant in this setting, as particles are confined by the externaltrap. The trapping frequency Ω controls the density and sizeof the system.Interestingly, the lattice spacing is again robust in this set-ting for the range of trapping frequencies studied here. As Ω increases, the lattice constant does not change significantlyas in non-cluster-forming systems [52, 53], but rather parti-cles form larger, more populated clusters, and hence a smallerlattice for the same number of particles. See Fig. 10 for typ-ical equilibrium states at four representative temperatures for Ω = 1 . The lattice constant of about 1.465 is very similarin all cases. The cluster sizes are no longer uniform due tothe trap. Instead they are larger at the center of the trap andsmaller at the edge. The problem is also related to the inter-vortex potentials with long-range attraction and multiple inter-mediate repulsive length scales, which is expected in layeredsystems [21].As a consequence the system has many energy scales. Upona slow cooling the particles order first at the center of thetrap, and then gradually in smaller clusters at the edge; c.f.with a discussion of melting of ordinary vortex lattices in atrap [54, 55]. This can be understood from the phase dia-gram shown in Fig. 4. For intermediate temperatures such as T = 0 . , the system can remarkably have a solid cluster-crystalline core surrounded by a cloud of fluid. This inho-mogeneous melting is an interesting phenomenon that can beexperimentally realized. IV. CONCLUSIONS In this work we studied a system of monodisperse particlesinteracting via an ultrasoft potential in 2D. We find that sim-ulated annealing is not effective to reach thermal equilibriumeven for this clean system without disorder, while it is still a0 -8-4048 -8 -4 0 4 8 T = 1 . -8-4048 -8 -4 0 4 8 T = 0 . -8-4048 -8 -4 0 4 8 T = 0 . -8-4048 -8 -4 0 4 8 T = 0 . (a) (b)(c) (d)FIG. 10: Monodisperse particles in an external harmonic trap ofstrength Ω = 1 with 400 particles [(a-d)]. Note the same latticestructure as the free case, but with very uneven cluster sizes. Theclusters are larger at the center of the trap, and smaller at the edgeof the lattice. The broad range of cluster sizes leads to nonuniformcrystallization or melting, as larger clusters order first followed bysmaller clusters upon further cooling. In particular, at intermediatetemperatures such as T = 0 . , the system remarkably has a solidcore with a shell of liquid fluid around it. useful tool for optimizations. Under these conditions paral-lel tempering is shown to be more efficient than populationannealing. We have presented a simple mean-field characteri-zation of the ground states and compared its predictions withMonte Carlo results, finding reasonably good qualitative andquantitative agreement.With the resolutions that we can achieve, the cluster liquidto cluster crystal phase transition is first order and occurs ina single step. We do not find indications for formation of anintermediate hexatic phase, which in contrast is the case fornon-cluster-forming systems. Our analysis is based on equi-librium energy histograms near the transition, without usingany cluster-related parameter. The results are relevant for theproblem of melting transitions in systems of soft particles andthe problem of vortex melting phase transition in supercon-ductors with multi-scale cluster-forming inter-vortex forces[13, 21, 29]. Acknowledgments W.W. and E.B. acknowledge support from the SwedishResearch Council Grant No. 642-2013-7837 and the GoranGustafsson Foundation for Research in Natural Sciences andMedicine. M.W. and R.D.M. acknowledge support from theSwedish Research Council Grant No. 621-2012-3984. Thecomputations were performed on resources provided by theSwedish National Infrastructure for Computing (SNIC) at theNational Supercomputer Centre (NSC) and the High Perfor-mance Computing Center North (HPC2N). [1] B. I. Halperin and D. R. Nelson, Theory of Two-DimensionalMelting , Phys. Rev. Lett. , 121 (1978).[2] D. R. Nelson and B. I. Halperin, Dislocation-mediated meltingin two dimensions , Phys. Rev. B , 2457 (1979).[3] J. M. Kosterlitz, Kosterlitz–Thouless physics: a review of keyissues , Rep. Prog. Phys. , 026001 (2016).[4] A. J. Armstrong, R. C. Mockler, and W. J. O’Sullivan, Isothermal-Expansion Melting of Two-Dimensional ColloidalMonolayers on the Surface of Water , J. Phys. Cond. Mat. ,1707 (1989).[5] E. P. Bernard and W. Krauth, Two-step melting in two dimen-sions: first-order liquid-hexatic transition , Phys. Rev. Lett. ,155704 (2011).[6] S. C. Kapfer and W. Krauth, Two-Dimensional Melting: FromLiquid-Hexatic Coexistence to Continuous Transitions , Phys.Rev. Lett. , 035702 (2015).[7] S. Prestipino and F. Saija, Hexatic phase and cluster crystalsof two-dimensional GEM4 spheres , The Journal of ChemicalPhysics , 184502 (2014).[8] M. Montes-Saralegui and G. Kahl, Obtaining equilibrium statesin ultrasoft cluster forming systems using a combined thermo-and barostat , J. Phys. Cond. Mat. , 325102 (2015).[9] G. Malescio and G. Pellicane, Stripe phases from isotropic re-pulsive interactions , Nature Materials , 97 (2003).[10] M. A. Glaser, G. M. Grason, R. D. Kamien, A. Koˇsmrlj, C. D. Santangelo, and P. Ziherl, Soft spheres make more mesophases ,EPL (Europhysics Letters) , 46004 (2007).[11] C. J. Olson Reichhardt, C. Reichhardt, and A. R. Bishop, Struc-tural transitions, melting, and intermediate phases for stripe-and clump-forming systems , Phys. Rev. E , 041502 (2010).[12] A. Ikeda and K. Miyazaki, Glass Transition of the Monodis-perse Gaussian Core Model , Phys. Rev. Lett. , 015701(2011).[13] R. D´ıaz-M´endez, F. Mezzacapo, W. Lechner, F. Cinti,E. Babaev, and G. Pupillo, Glass Transitions in MonodisperseCluster-Forming Ensembles: Vortex Matter in Type-1.5 Super-conductors , Phys. Rev. Lett. , 067001 (2017).[14] M. Zu, P. Tan, and N. Xu, Forming quasicrystals by monodis-perse soft core particles , Nature Communications , 2041(2017).[15] L. Caprini, E. Hern´andez-Garc´ıa, and C. L´opez, Cluster crys-tals with combined soft- and hard-core repulsive interactions ,Phys. Rev. E , 052607 (2018).[16] J. C. P`amies, A. Cacciuto, and D. Frenkel, Phase diagramof Hertzian spheres , The Journal of Chemical Physics ,044514 (2009).[17] W. L. Miller and A. Cacciuto, Two-dimensional packing of softparticles and the soft generalized Thomson problem , Soft Mat-ter , 7552 (2011).[18] K. Barkan, M. Engel, and R. Lifshitz, Controlled Self-Assembly of Periodic and Aperiodic Cluster Crystals , Phys. Rev. Lett. , 098304 (2014).[19] R. D´ıaz-M´endez, F. Mezzacapo, F. Cinti, W. Lechner, andG. Pupillo, Monodisperse cluster crystals: Classical and quan-tum dynamics , Phys. Rev. E , 052307 (2015).[20] F. Cinti, T. Macr`ı, W. Lechner, G. Pupillo, and T. Pohl, Defect-induced supersolidity with soft-core bosons , Nature Communi-cations , 3235 (2014).[21] Q. Meng, C. N. Varney, H. Fangohr, and E. Babaev, Phasediagrams of vortex matter with multi-scale inter-vortex inter-actions in layered superconductors , J. Phys. Cond. Mat. ,035602 (2016).[22] Q. Meng, C. N. Varney, H. Fangohr, and E. Babaev, Honey-comb, square, and kagome vortex lattices in superconductingsystems with multiscale intervortex interactions , Phys. Rev. B , 020509 (2014).[23] C. N. Varney, K. A. Sellin, Q.-Z. Wang, H. Fangohr, andE. Babaev, Hierarchical structure formation in layered super-conducting systems with multi-scale inter-vortex interactions ,J. Phys. Cond. Mat. , 415702 (2013).[24] C. N. Likos, A. Lang, M. Watzlawek, and H. L¨owen, Crite-rion for determining clustering versus reentrant melting behav-ior for bounded interaction potentials , Phys. Rev. E , 031206(2001).[25] C. N. Likos, Why do ultrasoft repulsive particles cluster andcrystallize? Analytical results from density-functional theory ,The Journal of Chemical Physics , 224502 (2007).[26] E. Babaev and M. Speight, Semi-Meissner state and neithertype-I nor type-II superconductivity in multicomponent super-conductors , Phys. Rev. B , 180502 (2005).[27] E. Babaev, J. Carlstr¨om, and M. Speight, Type-1.5 Supercon-ducting State from an Intrinsic Proximity Effect in Two-BandSuperconductors , Phys. Rev. Lett. , 067003 (2010).[28] M. Silaev and E. Babaev, Microscopic theory of type-1.5 su-perconductivity in multiband systems , Phys. Rev. B , 094515(2011).[29] E. Babaev, J. Carlstr¨om, M. Silaev, and J. Speight, Type-1.5superconductivity in multicomponent systems , Physica C: Su-perconductivity and its Applications , 20 (2017).[30] S. F. Edwards and P. W. Anderson, Theory of spin glasses , J.Phys. F: Met. Phys. , 965 (1975).[31] D. Sherrington and S. Kirkpatrick, Solvable model of a spinglass , Phys. Rev. Lett. , 1792 (1975).[32] A. P. Young, ed., Spin Glasses and Random Fields (World Sci-entific, Singapore, 1998).[33] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization bysimulated annealing , Science , 671 (1983).[34] R. H. Swendsen and J.-S. Wang, Replica Monte Carlo simula-tions of spin glasses , Phys. Rev. Lett. , 2607 (1986).[35] C. Geyer, in Computing Science and Statistics: 23rd Sympo-sium on the Interface , edited by E. M. Keramidas (InterfaceFoundation, Fairfax Station, 1991), p. 156.[36] K. Hukushima and K. Nemoto, Exchange Monte Carlo methodand application to spin glass simulations , J. Phys. Soc. Jpn. ,1604 (1996).[37] K. Hukushima and Y. Iba, in The Monte Carlo method inthe physical sciences: celebrating the 50th anniversary of theMetropolis algorithm , edited by J. E. Gubernatis (AIP, 2003), vol. 690, p. 200.[38] E. Zhou and X. Chen, in Proceedings of the 2010 Winter Sim-ulation Conference (WSC) (Springer, Baltimore MD, 2010), p.1211.[39] J. Machta, Population annealing with weighted averages: AMonte Carlo method for rough free-energy landscapes , Phys.Rev. E , 026704 (2010).[40] W. Wang, J. Machta, and H. G. Katzgraber, Population anneal-ing: Theory and application in spin glasses , Phys. Rev. E ,063307 (2015).[41] W. Wang, J. Machta, and H. G. Katzgraber, Evidence againsta mean-field description of short-range spin glasses revealedthrough thermal boundary conditions , Phys. Rev. B , 184412(2014).[42] L. Y. Barash, M. Weigel, M. Borovsk´y, W. Janke, and L. N.Shchur, GPU accelerated population annealing algorithm ,Computer Physics Communications , 341 (2017).[43] C. Amey and J. Machta, Analysis and optimization of popula-tion annealing , Phys. Rev. E , 033301 (2018).[44] A. Barzegar, C. Pattison, W. Wang, and H. G. Katzgraber, Op-timization of population annealing Monte Carlo for large-scalespin-glass simulations , Phys. Rev. E , 053308 (2018).[45] P. Tan, M. Steinbach, and V. Kumar, Introduction to Data Min-ing (Addison-Wesley Longman Publishing Co., Inc., Addison-Wesley, Reading, 2006).[46] R. D´ıaz-M´endez, G. Pupillo, F. Mezzacapo, M. Wallin, J. Lid-mar, and E. Babaev, Phase-change switching in 2d via soft in-teractions , Soft Matter , 355 (2019).[47] B. M. Mladek, P. Charbonneau, and D. Frenkel, Phase Coexis-tence of Cluster Crystals: Beyond the Gibbs Phase Rule , Phys.Rev. Lett. , 235702 (2007).[48] K. Zhang, P. Charbonneau, and B. M. Mladek, Reentrant andIsostructural Transitions in a Cluster-Crystal Former , Phys.Rev. Lett. , 245701 (2010).[49] For a more detailed comparison of particle liquids and clus-ter liquids, please see the movie: .[50] N. Henkel, F. Cinti, P. Jain, G. Pupillo, and T. Pohl, Super-solid Vortex Crystals in Rydberg-Dressed Bose-Einstein Con-densates , Phys. Rev. Lett. , 265301 (2012).[51] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez, The Defocusing Nonlinear Schr¨odinger Equation:From Dark Solitons to Vortices and Vortex Rings (SIAM,Philadelphia, 2015).[52] M. Cerkaski, R. G. Nazmitdinov, and A. Puente, Thomson ringsin a disk , Phys. Rev. E , 032312 (2015).[53] B. Ash, J. Chakrabarti, and A. Ghosal, Static and dynamic prop-erties of two-dimensional Coulomb clusters , Phys. Rev. E ,042105 (2017).[54] S. Kragset, E. Babaev, and A. Sudbø, Thermal Fluctuationsof Vortex Matter in Trapped Bose-Einstein Condensates , Phys.Rev. Lett. , 170403 (2006).[55] S. Kragset, E. Babaev, and A. Sudbø, Effects of boundariesand density inhomogeneity on states of vortex matter in Bose-Einstein condensates at finite temperature , Phys. Rev. A77