Competing structures in a minimal double-well potential model of condensed matter (2025)

Julyan H. E. Cartwright^1,2julyan.cartwright@csic.es  Bruno Escribano³bruno.escribano.salazar@gmail.com  Sándalo Roldán-Vargas⁴sandalo@ugr.es  C. Ignacio Sainz-Díaz¹ignacio.sainz@iact.ugr-csic.es¹Instituto Andaluz de Ciencias de la Tierra, IACT-CSIC, 18100 Armilla, Granada, Spain
²Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, 18071 Granada, Spain
³Centro de Astrobiología (CAB), CSIC-INTA, 28850 Torrejón de Ardoz, Madrid, Spain,
⁴Department of Applied Physics, Faculty of Sciences, Universidad de Granada, 18071 Granada, Spain

Much research, both laboratory experiments and numerical simulations, is taking place on substances that show multiple liquid and amorphous solid phases. These include water Kim etal. (2020); Rosu-Finsen etal. (2023), silicon Morishita (2004); Ganesh and Widom (2009); Sastry and Angell (2003); Treacy and Borisenko (2012), carbon Togaya (1997); Glosli and Ree (1999); Sachan etal. (2019), phosphorus Katayama etal. (2000); Morishita (2001); Monaco etal. (2003),silica Lacks (2000); Trave etal. (2002), triphenyl phosphite Kurita and Tanaka (2004), carbon dioxide Bai etal. (2018), and calcium carbonate Cartwrightetal. (2012a). Noted phenomena are long- or medium-range order despite the system being amorphous, spatial heterogeneity, both dynamic and static, and local orientation Lan etal. (2021). To rationalize these features, which are common to many amorphous solid materials, it has been proposed that a particle pair potential having the form of a double well may be the generic mechanism for this rich phenomenology, termed in general polyamorphism Buldyrev etal. (2002); Wilding etal. (2006). In particular, confined and interfacial two-dimensional materials are being investigated in theory, simulations,and experiment not only for their wide range of applications but also because of their fundamental role in clarifying different surface phenomena Li etal. (2021). In this Letter we focus on two dimensions and demonstrate that long-range order that extends further than the typical liquid order, static spatial heterogeneity, and local orientation emerge from a minimal model with a double well and isotropic interaction.

More complex approaches with a double-well potential, e.g., models based on molecular thermodynamics, have shown water-like anomaliesand/or a polyamorphic transition Franzese etal. (2001); Malescio (2007); Franzese (2007); Rosu-Finsen etal. (2023); Mollica etal. (2022).Core-softened potentials such as the Hemmer-Stell-Jagla modeldisplay two phase transitions and anomalous-waterlike propertiesHemmer and Stell (1970); Stell and Hemmer (1972); Jagla (1999). Likewise, previous work with a minimal model of water, the so-called Mercedes-Benz model that incorporates directional bonding by representing water molecules as two-dimensional Lennard-Jones disks with three equivalent hydrogen-bonding arms disposed at 120 degrees, displays a first-order phase transition between a crystalline phase and a high-density amorphous phase as well as a reversible transformation between two amorphous structures of high and low density Cartwrightetal. (2012b).Here we investigate an even simpler model in two dimensions containing only a double well with no directionality in the particle bonding mechanism and without thermodynamical considerations. The choice of a double well potential differs from previous mechanisms introducing amorphous structure such as the use of particles with different (fixed) sizes Brüning etal. (2009). Here a given particle can create bonds at different distances.

With that objective in mind, we simulate the growth of a two-dimensional amorphous solid material by depositing particles on a plane surface. Particles are deposited one at a time randomly with uniform probability on the simulated plane, and can interact with each other forming permanent bonds by proximity which are not directionally constrained. In between depositions, we allow the system to relax and stabilize using a second-order Euler integrator for the equations of motion, for which we assume underdamping to avoid numerical instabilities.Typically, we allowed 10 time-steps of the Euler integrator between depositions of new particles. We did not study the effect of longer relaxations which would simulate slower depositions.

We use a double-well potential for the interaction between a pair of bound particles separated a distance $r$,

where $r_{a}$ and $r_{b}$ are the equilibrium distances for the two wells; Fig.1. For simplicity, we choose $k=1$ and model the particles as point masses with mass $m=1$, defining the unit of length $r_{0}={(r_{a}+r_{b})/2=1}$ as the distance of the potential barrier, which is symmetrically placed between the two wells. In addition, we introduce $r_{max}$ as the maximum distance for bond formation, which is necessary for avoiding the divergent asymptotic behavior of the potential, and set $\epsilon=\frac{1}{2}k(r_{max}-r_{a})^{2}(r_{max}-r_{b})^{2}$. Hence we have particle bond lengths not very far from $r_{a}$ and $r_{b}$. In this respect, we also impose bonds to be irreversible, in the sense that, once two particles are bonded, they could interact by $V(r)$ even for distances greater than $r_{max}$, although in the explored systems this situation is marginal. We control the number of bonds per particle by imposing a maximum coordination number, $n$; each particle can bond with at most $n$ neighbors.This is imposed during the neighbor search algorithm, which is executed at every time-step of the Euler method for the integration of the equations of motion. By altering the maximum coordination number, we obtain different competing structures without appealing to thermodynamic considerations based on the interplay between entropically and energetically favorable structures Roldan-Vargas etal. (2013).

Once a particle reaches $n$ bonds it becomes saturated and does not create additional bonds with other particles in its vicinity. In this respect, we also allow non-bonded repulsive interactions for $r\leq r_{min}$, given by the potential $V_{unbonded}(r)=V(r)$ if $r\leq r_{min}$ and $V_{unbonded}$ = 0 if $r>r_{min}$ (red line in Fig.1). This repulsive interaction represents a form of excluded volume, keeping non-bonded particles from a significant overlap between each other for $r\leq r_{min}$ while moving during the system relaxation after each deposition. In this respect, we also note that when a new particle is randomly deposited at a distance $r\leq r_{min}$ from a previously deposited particle, we remove that new particle and try another random deposition.

We set the distances in our model as $r_{min}=0.7$, $r_{a}=0.8$, $r_{b}=1.2$ and $r_{max}=1.3$. With these choices we: i) place the two wells, as well as $r_{min}$ and $r_{max}$, symmetrically with respect to the potential barrier (thus effectively reducing the complexity of our model), ii) impose a ratio ${r_{min}/r_{a}}=0.875$ similar to the one associated to a Lennard-Jones interaction between the Lennard-Jones particle diameter and the corresponding Lennard-Jones minimum Lennard-Jones (1931), and iii) allow any two bonded particles to cross the barrier in either sense, and jump from one well to the other during the relaxation after depositions. To investigate the effect of the occupied surface, we introduce a packing-fraction parameter,$\phi={N\pi(r_{min}/2)^{2}/L^{2}},$where $L=100$ is the fixed side length of our square simulation box, with periodic boundary conditions, and $N$ the total number of particles, which we alter to explore different packing fractions. The results presented here have been tested for large enough values of $L$ and $N$ to remove finite-size effects. We present results for $\phi=0.2$, $0.3$, and $0.4$ — $N=5197$, $7795$, and $10394$ — that is, from a dilute system to a packing fraction close to the random packing fraction of equal non-overlapping disks Hinrichsen etal. (1990). To isolate the structural effect of the maximum coordination number, we cover, for any given $\phi$, a wide range of $n\in\{2,3,4,5,6\}$. We analyze 40 independent simulations for each pair $(\phi,n)$.

We begin analyzing the system structure by looking at the radial distribution function $g(r)$Hansen and McDonald (1990) as a function of $\phi$ and $n$, Fig.2. Notice the emergence of two main peaks placed at $r_{a}$ and $r_{b}$, the two potential wells, whose height increases with $\phi$, with independence of $n$. Second, for any given $\phi$, $g(r)$ shows a more obvious (amorphous) order when increasing $n$. This is again manifested through the height of the two main peaks. More interestingly, we observe how for low $n$ ($=2,3$), and for any given $\phi$, the first of the two main peaks, placed at $r_{a}$, is more likely than the second peak placed at $r_{b}$; see the yellow and black dashed lines in Fig.2. However for high coordination numbers, $n=5$ and $6$, the second peak is more likely than the first; see the green and blue dashed lines in Fig.2. In simple terms, for high coordination numbers particles tend to bond to their neighbors at larger distances. Presumably, this more frequent bonding with the second well has a clear effect on the medium-range order for $n=5$ and $6$, especially at $\phi=0.3$ and $0.4$, where the medium-range system structure is manifested through well developed peaks for $r>2$; see insets in Fig.2. The phenomenology shown by $g(r)$ already displays the versatility of our model, which produces a medium-range order that extends further than the typical liquid structure as well as a means to control the first-neighbor structure by controlling $n$. This non-trivial structure can be connected with previous — and more complex — models and experiments. For instance: i) Koga et al.Koga etal. (2000) analyzed a quasi-2D amorphous water model by molecular-dynamics simulations using a TIP4P force field which incorporates a refined description of the internal geometry of the water molecules. They also found two main peaks and a medium-range order similar to those shown by our model, ii) similar profiles for $g(r)$ have been reported recently by Negi et al.Negi etal. (2022) for a two-dimensional ice model using a sophisticated density-functional theory approach, iii) Treacy et al.Treacy and Borisenko (2012) studied the amorphous phase of silicon by electron-diffraction experiments, finding a $g(r)$ profile consistent with our results: two main peaks and additional small peaks at longer distances; in particular, they found a coordination number of 3.8 for a narrow 3D system confined between two walls.

This non-trivial order is manifested more clearly when looking at the structure factor $S(q)$Hansen and McDonald (1990), Fig.3. First, we observe two main peaks around $q$ $\approx$ 6 and 11, in accordance with the two main peaks shown by $g(r)$ for the two potential wells. As in $g(r)$, these two peaks are more obvious upon increasing $\phi$ and $n$. Long-range order is now seen at low $q$, at distances even longer than those manisfested by $g(r)$. Again, this long-range order is obvious for high coordination numbers and, especially, at intermediate $\phi$. Thus, for $q<$ 4 and $\phi=0.3$, the system develops a plateau at high $n$ ($=5,6$) which tends to disappear at low coordination numbers, $n=2$, $3$; see insets in Fig.3. This plateau, particularly apparent around $q$ $\approx$ 2, corresponds to distances of the order of $r\geq 3$, and can be interpreted as the emergence of long-range spatial heterogeneities in the particle density. The long-range order indeed extends quite far ($r\approx 6$, $q\leq 1$), where the limit of $S(q)$ when $q\rightarrow{0}$ produces higher values upon increasing $n$, making the system more compressible, in analogy with systems controlled by temperatureHansen and McDonald (1990).The plateau disappears at $\phi=0.4$ for any $n$ and the system recovers a liquid-like spatial homogeneity, despite maintaining an obvious local medium-range order. At these high packing fractions ($\phi=0.4$) the long-range spatial heterogeneities shown by $S(q)$ for $\phi=0.3$ disappears presumably due to the absence of empty space. $S(q)$ is directly measurable in scattering experimentsHansen and McDonald (1990) and, in particular, the structure factor plots we obtained are consistent with experimental investigations based on X-ray diffraction and neutron scattering on different tetrahedral glasses and amorphous materials such as silicon and glassy selenium Wilding etal. (2006); Dyre (2006). In these works, double-peak structures are observed, the two main peaks being found at distances consistent with our model. Daisenberger et al.also found a similar $S(q)$ profile for silicon with X-ray diffraction experiments Daisenberger etal. (2007); Loerting etal. (2009). In their work, they control the system pressure in analogy with our control of the packing fraction. In addition, Mollica et al.Mollica etal. (2022) used MD simulations of significant complexity to study amorphous ice in silico using a TIP4P/ice force field where the simulated sample reduces to significantly narrow widths. In spite of the vast differences between their model and ours, the results are consistent: they show a $S(q)$ for the oxygen-‐oxygen correlations with a double-peak structure and a behavior at low $q$ that resembles the one we observe at intermediate packing fractions $\phi=0.3$.

To have a clearer intuition of the structure manifested through $S(q)$ and $g(r)$ at intermediate $\phi$, in Fig.4 (above) we show snapshots of a large control surface of two stable configurations with the same density of particles for $\phi=0.3$ with $n=3$ (left) and $n=6$ (right), where bonds in the first (second) well are depicted as black (red) lines. Although in both cases the system is not fully connected into a single cluster, the clusters we see for $n=6$ spread to larger areas and involve more particles. This is expected since at high $n$, particles can bind to more neighbors and therefore propagate their connections further. Insets in Fig.3 show for all the $\phi$ values investigated the mean number of bonds per particle, as well as the potential energy per particle, $V$/particle, as a function of $n$. There we see how, for $\phi=0.3$, particles with $n=3$ present around 2.5 bonds per particle on average whereas this number reaches a value around 4.5 with $n=6$. The long-range structure we see for $S(q)$ at low $q$ ($\approx 2$) and intermediate $\phi$ can also be observed through the particle configurations. This is what we show in Fig.4(bottom), where the previous snapshots have been coarse-grained by dividing the control surface into squares of side 3.58 ($\approx$ 2$\pi$/$q$, where $q\approx 2$). We assign a color to each square according to the number of particles it contains. Whereas for $n=3$ (left) we see that colors are mostly uniform, pointing to a uniform local density, for $n=6$ (right) we clearly see strong local deviations from the average density with high density (light yellow) and low density (violet) squares. This observation points to the emergence of local spatial heterogeneity at high $n$ for the length scale where the plateau in $S(q)$ is observed. Both observations, the emergence of large clusters and heterogeneous local density, are also manifested through the contour lines following constant densities that we show in Fig.4 (below). Since our control surface has a $10\times 10$ square resolution, we applied a B-spline smoothing method to the contour lines to improve their visibility MacCallum and Zhang (1986). Whereas for $n=3$ (left) contour lines with a given density typically extend to small areas, for $n=6$ (right) we see some contour lines that prolong to large distances. In particular, we have highlighted two contour lines, corresponding to low density, that cross the entire control surface. Similar snapshots to Fig.4 have been observed in the deposition of a monolayer of water on hydrophilic surfaces Brovchenko and Oleinikova (2008).

We lastly analyze the local orientation in our system at intermediate $\phi$, having in mind the fact that interaction in our model is isotropic and, therefore, does not contain explicit directionality. Fig.5 presents the local angular distribution according to the construction we show by the sketch. For any given particle having at least two neighbors (see red particle in the sketch), we construct all the triangles formed by the particle and any pair of its adjacent neighbors (see blue triangles in the sketch). We then collect all the angles of the triangles associated with all the particles in the system having at least two neighbors and study the different normalized histograms corresponding to each $n$ value. We see that all the histograms present four well developed peaks placed, from left to right, around $\pi/4$, $\pi/3$, $2\pi/5$, and $\pi/2$; as well as a rather uniform spectrum at large angles in the form of a long tail. These distributions show the interplay present in our system between local polycrystallinity, expressed by well defined angles, and amorphous structure, given by a long tail at large angles.

The well-defined angles we observe can be typically associated with specific triangular structures (see Fig.4 (top)): i) $\pi/4$ and $\pi/2$ could be naturally explained by the presence of square-like structures (with the red particle at a vertex) with a right angle ($\pi/2$) and two angles ($\pi/4$) defining the diagonal; ii) $\pi/3$ would be associated with equilateral triangles that appear dispersed or as part of a hexagonal structure (with the red particle in the center); and iii) $2\pi/5$ being associated with pentagonal structures. However this idealized interpretation does not take into account the two distances, potential wells, defining the system interaction. These two distances redound to the width of the peaks we observe as well as in their precise location. Although no directional interaction is imposed in our model, the revealed peak structure becomes more apparent when increasing $n$, whereas the corresponding tail at large angles decreases, leading to better defined local polycrystalline structures. This is in clear contrast to the broad and centered distribution observed for non-overlapping disks randomly deposited when analyzing their Voronoi polygons Hinrichsen etal. (1990).

We have introduced a minimal athermal model in two dimensions with isotropic interactions and three basic ingredients: a double-well bonding mechanism, an excluded volume, and a maximum coordination number. The model shows that even without any thermodynamical considerations, it is possible toproduce all the relevant features observed in several amorphous substances such as long-range order, spatial heterogeneity, and local polycrystallinity. Our results are applicable to a variety of systems including 2D liquid water and ice Meyer and Stanley (1999); Koga etal. (2000); Johnston etal. (2010); Chen etal. (2016); Ma etal. (2022); Negi etal. (2022); Mollica etal. (2022); Brovchenko and Oleinikova (2008); Bartels-Rauschetal. (2012), silicon Sastry and Angell (2003); Treacy and Borisenko (2012); Daisenberger etal. (2007); Loerting etal. (2009), and different glassy and amorphous solids Dyre (2006); Loerting etal. (2009); Wilding etal. (2006). These results support the hypothesis that a double-well potential may be sufficient to give rise to the coexistence of different types of local competing structures emerging in materials displaying polyamorphism Buldyrev etal. (2002); Wilding etal. (2006); Rosu-Finsen etal. (2023); Mollica etal. (2022); Shi and Tanaka (2020).

B.E. acknowledges support from grants PCIN-2017-098, PTA2020-018247-I and PID2020-118974GB-C21 from the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033. S.R.-V. acknowledges support from the European Commission through the Marie Skłodowska-Curie Individual Fellowship 840195-ARIADNE.