Attraction Controls the Entropy of Fluctuations in Isosceles Triangular Networks

We study two-dimensional triangular-network models, which have degenerate ground states composed of straight or randomly-zigzagging stripes and thus sub-extensive residual entropy. We show that attraction is responsible for the inversion of the stable phase by changing the entropy of fluctuations around the ground-state configurations. By using a real-space shell-expansion method, we compute the exact expression of the entropy for harmonic interactions, while for repulsive harmonic interactions we obtain the entropy arising from a limited subset of the system by numerical integration. We compare these results with a three-dimensional triangular-network model, which shows the same attraction-mediated selection mechanism of the stable phase, and conclude that this effect is general with respect to the dimensionality of the system.


Introduction
Geometrically-constrained systems may show peculiar features compared to their unconstrained counterparts. In particular, geometric constraints can lead to frustration because the system cannot simultaneously minimize all local interaction energies, or free energies. Frustrated systems usually show a degenerate ground state and then they may posses a residual entropy. Frustration is relevant in physical and biological systems that range from water [1] and spin ice [2,3] to magnets [4], magnetic island [5], high transition-temperature superconductors [6], elastic beams [7,8] and colloids [9][10][11][12]. The possibility to control colloidal interparticle interactions and to visualize and manipulate each particle and follow its motion in both space and time makes colloidal suspension a powerful tool to study phenomena in condensed-matter physics, ranging from glass formers [13], to crystals and gels [14].
A prototypical geometrically-confined system is composed of short-range repulsive colloids confined in a slit pore of two plates. Varying density and plate separation, discontinuous phase transitions between layered, buckled, rhombic and adaptive prism crystal structures occur [15][16][17]. In the case of a slit pore with a plates interdistance slightly larger than a colloid diameter, when density approaches the close-packing value ρ cp , colloids, due to their free-volume-dominated free energy, tend to touch opposite walls, giving rise to effective antiferromagnetic interactions [9,18], and to glassy dynamics [19]. Multiple configurations corresponding to the same ρ cp can be obtained by alternating straight stripes of up and down spheres (Figure 1a) or by any set of zigzagging stripes (Figure 1b).
This ground-state degeneracy implies a subextensive residual entropy at ρ cp (S 0 ∼ √ N, with N the number of particles in the system) [9] so that the residual entropy per particle tends to zero in the thermodynamic limit. At ρ = ρ cp , in the straight-or zigzagging-stripes configuration, a colloidal sphere is surrounded by two colloids touching the same wall (giving rise to frustrated bonds) and by four colloids touching the opposite wall (satisfied bonds), except for the limit case of α = 90 • , in which case every sphere is surrounded by four colloids touching the same wall and four colloids touching the opposite wall, where α is defined in Ref. [20]. The terms frustrated and satisfied bonds refer both to the packing consideration of neighboring spheres wanting to touch opposite walls, and also from the obvious connection to an antiferromagnetic spin-1/2 Ising model [18]. In this case, nearest neighbors at opposite spin states satisfy the antiferromagnetic interaction between them, while those at the same spin state represent the frustration of this Ising model on the triangular lattice [21]. The three-dimensional (3D) network of links connecting the centers of neighboring spheres is composed of tilted equilateral triangles, since the 3D distance between the centers of each pair of contacting spheres is equal to the sphere diameter. When the centers of the colloids are projected on the plane parallel to the plates, it reduces to a 2D network composed of isosceles triangles. In this network, the longer link of each triangular plaquette corresponds to a frustrated bond, while the two shorter links correspond to satisfied bonds. Straight and zigzagging stripes are the only configurations corresponding to a 2D network of isosceles triangles which can tile the plane [9]. maximally zigzagging or bent configuration. Shells, which order is indicated with n s , are denoted by increasingly darker colors for increasing n s . Thicker, gray lines correspond to springs of longer (for α > π/3) or shorter (for α < π/3) length at rest; (c,d) are the unit cell for the straight-and bent-stripes confgirations, respectively. Numbers associated to particles correspond to the particle positions in tables in Section 2; (e) Parameters associated to every plaquette: a, b, c and α.
This geometric mechanism underlying the ground state of the buckled colloidal system composed of straight or zigzagging stripes, has been realised, experimentally also by packing a granular system in a container under the effect of gravity [22,23], and theoretically with spins starting from the Wannier antiferromagnetic Ising model on a triangular lattice [21] by allowing for the lattice to elastically deform [24,25]. Zigzagging stripes patterns have been found also in the Ising model on an anisotropic triangular lattice [26]. A relevant phenomenon involving the selection of one between two competing phases related to the distorsion of equilateral into iscosceles triangles and due to angular frustrations concerns stiff, nematically ordered, polymer molecules such as DNA [27]. The Ising antiferromagnet on a deformable triangular lattice has the same degeneracy of the ground state and subextensive entropy at T = 0 as the colloidal system at ρ = ρ cp (that is S 0 ∼ √ N, equivalent to the N 2/3 scaling found in Perovskite Oxynitrides [28]). In the following we refer to T = 0 as temperature being arbitrarily close to zero. Indeed, it has been shown that the third law of thermodynamics, implying the unattainability of absolute zero temperature in a finite number of steps and within a finite time, holds for arbitrary classical or quantum systems or involving infinite-dimensional reservoirs [29].
For temperature slightly larger than zero, the degeneracy is removed by the order-by-disorder effect [30][31][32][33][34][35][36][37] and in the elastic Ising model straight stripes represent the stable phase [20,25], while bent stripes are selected in the colloidal system for ρ < ρ cp when colloids are modeled using hard or soft repulsive potentials [20]. By tuning the attractive vs repulsive components of an asymmetric power-law potential used to model colloids, we found that the stable phase in the colloidal system can be turned from bent to straight stripes for attraction strong enough compared to repulsion [20]. We established a connection between the effect of the attraction on the phase stability and the packing of hard spheres and their entropy. We showed that other parameters of these systems are irrelevant to the phase stability, as for example their dimensionalities. Indeed, the elastic Ising antiferromagnet is defined in 2D, while the colloidal system is 3D or quasi-2D due to the buckling of the monolayer.
In this paper we study a 2D isosceles triangular network which shows the same ground-state degeneracy as the Ising elastic antiferromagnet at T = 0 or the colloidal monolayer at ρ = ρ cp . Increasing temperature above zero, the degeneracy is removed through the order-by-disorder effect and we find the straight-stripes phase to be selected for particles linked with harmonic interactions, while the bent-stripes phase is more stable if only the repulsive component of the harmonic inter-particle potential is considerd (that we call repulsive harmonic potential, described in detail in the following section). This result suggests that the inversion of the stable phase by adding attraction to repulsively-interacting particles in triangular networks is a general mechanism, irrespective of the dimensionality of the system.

Isosceles Triangular Network Model
The model we consider is composed of particles in a 2D triangular network linked with springs of two different lengths at rest such that every plaquette or triangle is formed by a longer edge 2a and two shorter identical edges c, and thus has a height b = √ c 2 − a 2 and a head angle α such that a = b tan(α/2) (see Figure 1e). We will first consider particles interacting through a harmonic potential and then will consider a repulsive harmonic potential, as defined below. The Hamiltonian of the system for harmonic inter-particle interaction can be written as where K is the spring constant, which is assumed to be identical for all springs, 1 ≤ m, n ≤ L, N = L 2 is the number of particles or nodes of the network and the index l runs over three of the six neighbors each particle has to avoid double counting of bonds. The positions of the particles are described by the coordinates {x i , y i }. At T = 0 multiple degenerate states, represented by straight or zigzagging stripes, minimize the system free energy. For T > 0 the system cannot jump from one configuration to another, but at T = 0 it is at mechanical equilibrium in every state under consideration. Therefore, the present model is not ergodic by construction. In the conclusion we will discuss entropy calculation in non ergodic systems. In order to study the stability of the system at T > 0, we will consider small fluctuations about the equilibrium position described by small displacements {u i , v i } of all particles in the straight and bent configurations. dr is the distance between particles i and j, and its square is thus given by is the length at rest of the spring linking particles i and j, which can take the values c or 2a (see Figure 1e, Tables 1 and 2). Since we consider the expansion around mechanical equilibrium, we will drop terms linear in du and dv and write: dr 2 = dr 2 0 + du 2 + dv 2 . Expanding to harmonic order the expression of dr that we get from Equation (2), we obtain The Hamiltonian of the straight-stripes configuration, with particle positions specified for the unit-cell in Table 1, is Using the relations du l = u l − u 0 and dv l = v l − v 0 we get The Hamiltonian of the bent-stripes configuration, with particle positions specified for the unit-cell in Table 2, is ∑ t,n du 2 10 + cos 2 αdu 2 50 + sin 2 αdv 2 50 − sin(2α)du 50 dv 50 + sin 2 ( where 1 ≤ t ≤ L/2. Indeed, the unit-cell of the bent stripe configuration includes two particles (see Figure 1d): particle 0, which represents the particles with odd m, and particle 1, which represents the particles with even m. Therefore, we set for particle 0, m = 2t − 1, and for particle 1, m = 2t.
Using the relations For a harmonic interparticle potential, the Hamiltonian for straight and bent stripes configurations we expanded around mechanical equilibrium takes the quadratic form: H = K ∑ m,n A m,n q m q n , where {q} = {u, v} represents small displacements about the equilibrium position of every particle. In the canonical ensemble the difference between the entropy per particle of the straight and bent configurations for such Hamiltonian is [20]: ∆s = (S s − S b )/N = 1/(2N) ln( A b / A s ) where the subscript s refers to straight and b to bent, and A is the determinant of A. The dimensionless matrix A depends only on the deformation angle α and on the zigzagging-stripe realization. In [20] we used a recursive method to obtain the matrix A in the case of the elastic Ising model for any subset of the network composed by shells of particles (see Figure 1). Here we apply the same method to the 2D spring network model. The number of particles n belonging to the shells up to n s is given by n = 1 + 3n s (n s − 1). In our shell-expansion calculation, these n particles are free to move, while the other N − n particles of the network are frozen in their equilibrium position. Increasing n, ∆s should converge to the exact result (see Figure 2a), which includes the simultaneous fluctuation of all particles in the system.
In Figure 2a we show ∆s for the 2D harmonic network model for a number of shells up to n s = 20, that is n = 1141 particles free to move. From it we can see that ∆s > 0 for every deformation angle α of the network and for every order n s of the expansion (except for very small deviations at small α and small n s ). Considering only one particle free to move, n s = 1, gives a qualitative indication on the behavior of ∆s for every orders of approximation. This rapid convergence with increasing n s gives confidence in this expansion method when we will apply it for purely repulsive interactions, for which we are technically much more limited in the number of particles that we may numerically calculate the simultaneous fluctuation of. Table 1. Distances between the neighboring particles and the central particle 0 in the unit cell of the straight-stripe configuration. Particle positions are graphically shown in Figure 1c. Table 2. Distances between the neighboring particles and 0 and 1 particles in the unit cell of the maximally zigzagging-stripe configuration. Particle positions are graphically shown in Figure 1d.
Now we consider the same 2D triangular network model, but for repulsive harmonic interactions, that is we consider the following Hamiltonian where θ() is the Heaviside step function. Namely, now each spring applies a restoring force only when compressed (dr l < dr 0 ) and there is no resistance to stretching (dr l > dr 0 ). In this case we have to numerically integrate the partition function in order to get the entropy of straight-and bent-stripe configurations. In Figure 2b we show ∆s for repulsive harmonic interactions for n = 1, 2, 3. Numerical calculations for n > 3 are beyond our computational reach (numerical integration for a number of variables larger than 6, in our case suffers from fluctuating results for a limited capacity in the precision). For n = 1 particle free to move, each particle can be equivalently chosen to be free. For n = 2 we consider the particles 0 and 1 (see Figure 1c,d). For n = 3 we consider the particles 0, 1, 2 and 0, 1, 4 for the straight and bent configurations, respectively (see Figure 1c,d). As for the 3D spring network of [20], we find that also in our 2D model ∆s < 0 for repulsive interactions. For the case n = 1 we can show how the inversion of the stable phase when turning from harmonic to repulsive harmonic potential depends on the spatial configurations of straight and bent stripes, and in particular of the angular distribution of neighboring particles around each particle in the network.
For n = 1 the computation of the canonical partition function for the repulsive harmonic system can be easily reduced to the integration of single variable functions using polar coordinates (ρ, γ) (see Figure 3). The Hamiltonian of the free particle 0 can be written as H r 0 = 1/2Kρ 2 ∑ 6 i=1 g i (α, γ) where ρ 2 g i (α, γ) is the contribution to H r 0 coming from the neighbor i of the particle 0, and the function g i depends on the coordinates of the particle i, as specified below. The canonical partition function is thus where β = 1/(K B T) is the Boltzmann factor and In this case we have ∆s = ln(I s (α)/I b (α)).

Figure 3.
Example of straight (a) and bent (b) configuration for n = 1 particle free to move. The deviation of the central particle from its equilibrium position is described by polar coordinates (ρ, γ).
The difference between the calculation of the harmonic and the repulsive harmonic partition function is that in the former case the functions g i contributes to the integral I(α) for any angle 0 ≤ γ ≤ 2π, while in the latter case every g i contributes to I(α) for a specific range of angles only. For the repulsive harmonic potential in the straight-stripes configuration we need to consider the azimuthal ranges coming from each one of the neighboring particles: Due to the symmetry of the straight-stripes configuration, thanks to which the reflection about the origin of each neighbor transforms it in another neighboring particle (see Figure 1c and Table 1), we can write the function f r s by just taking the contribution of every g i without the condition imposed by the Heaviside step functions and dividing it by 2, that is f r s (α, γ) = cos 2 γ + sin 2 (α/2 + γ) + sin 2 (α/2 − γ) = cos 2 γ(2 − cos α) + cos 2 (α/2) For the repulsive harmonic potential in the bent configuration we need to consider the azimuthal ranges coming from each one of the neighboring particles: In Figure 4 we show 1/ f for hamonic and repulsive harmonic interactions for straight and bent configurations as a function of the lattice deformation angle α and the azimuthal direction in space γ, over which we numerically integrate in order to get the value of the function I for that specific angle α, which in turn sets the entropy via Equation (9). From Figure 4 we can see, particularly for big angles α for which ∆s takes its larger values (see Figure 2b), that repulsion accentuates the contribution to the function I for bent stripes (corresponding in Figure 4 to a wider region composed of brighter colors, i.e., white, yellow and red, for repulsive harmonic over harmonic interaction in the case of bent stripes). More in general, we can say that repulsion accentuates differences in the contribution to the partition function and thus to the free energy between symmetric and asymmetric distribution of neighboring particles.

Conclusions
We studied a 2D triangular-network model composed of particles interacting through harmonic or repulsive harmonic springs. At T = 0 the ground state is degenerate and composed of straight or any set of zigzagging-stripes configurations. At T > 0 we found that the stable phase is composed of straight or bent stripes depending on the harmonic or repulsive harmonic nature of particle interaction, respectively. This selection mechanism of the stable phase through the order-by-disorder effect is equivalent to that observed in the colloidal [20] and Ising [25] antiferromagnets irrespective of the dimensionality of the system. This suggests that the phase inversion of isosceles triangular networks is controlled by the attraction component of the interparticle interaction. We suggest that this is due to the fact that repulsive interactions accentuate differences in the contribution to the free energy between symmetric and asymmetric distribution of neighboring particles, as we have shown for n = 1 free particle calculation.
Both the Ising antiferromagnet on a deformable triangular lattice and the 2D isosceles triangular network model at T = 0 from one side, and the colloidal monolayer at ρ = ρ cp from the other side, have the same degeneracy and subextensive entropy S ∼ √ N and thus a vanishing residual entropy per particle in the thermodynamic limit. At T > 0 for the triangular networks and at ρ < ρ cp for the colloidal monolayer this degeneracy is removed, but they can still have a residual entropy per particle for finite system size if the ergodicity is broken. Indeed, even at T > 0 or ρ < ρ cp a system may be trapped in a local minimum of the free-energy landscape and thermal fluctuations are not large enough for a small system to overcome energy barriers. For ergodic systems the time average of observables can be computed by using ensemble averages thanks to the Birkhoff theorem [38]. From the other hand, for non-ergodic systems, the phase space is divided into disjoined sets. In this case, states can be counted either following the kinetic view [39], for which only states visited by the system at the observational time scale are taken into account, or following the Edwards approach [40], for which all possible states are considered regardless of whether they are explored or not by the system. Recently, the Edwards hypothesis has been proved to be valid at the un-jamming point [41]. In thermal ergodic systems at equilibrium, the two sampling methods give the same result. Following the Edwards approach we can conclude that an indication of the presence of residual entropy in a system is given by the ergodicity breaking (which can be checked for generic temperature or density) instead of by the degeneracy of the ground state (defined for T = 0 or ρ = ρ cp only).