Attraction controls the entropy of fluctuations in isosceles triangular networks
Fabio Leoni Yair Shokef
cond-mat.soft







School of Mechanical Engineering and Sackler Center for
  Computational Molecular and Materials Science, Tel-Aviv University, 
  Tel-Aviv 69978, Israel 


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.


Attraction controls the entropy of fluctuations in isosceles
  triangular networks
    Yair Shokef[]
    
==================================================================================





§ 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 <cit.> and spin ice <cit.>
to magnets <cit.>, magnetic island <cit.>, high 
transition-temperature superconductors <cit.>, elastic
beams <cit.> and colloids 
<cit.>.  
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
<cit.>, to crystals and gels <cit.>.

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 <cit.>. 
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 <cit.>, and to
glassy dynamics <cit.>.  
Multiple configurations corresponding to the same ρ_cp can be
obtained by alternating straight stripes of up and down spheres
(Figure <ref>a) or by any set of zigzagging stripes 
(Figure <ref>b). 

This ground-state degeneracy implies a subextensive residual
entropy at ρ_cp (S_0∼√(N), with N the number of
particles in the system) <cit.> 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 colloids 
touching the same wall (giving rise to a frustrated bond) and by
colloids touching the opposite wall (satisfied bond).
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 <cit.>. 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
<cit.>.  
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 <cit.>. 

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 <cit.>, 
and theoreticelly with spins starting from the Wannier
antiferromagnetic Ising model on a triangular lattice
<cit.> by allowing for the lattice to elastically deform
<cit.>.  
Zigzagging stripes patterns have been found also in the Ising model
on an anisotropic triangular lattice <cit.>. 
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 <cit.>). 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 <cit.>.

For temperature slightly larger than zero, the degeneracy is removed by
the order-by-disorder effect
<cit.>   
and in the elastic Ising model straight stripes represent the stable
phase <cit.>, while bent stripes are selected in
the colloidal system for ρ<ρ_cp when colloids are modeled
using hard or soft repulsive potentials <cit.>. 
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 large enough compared to repulsion 
<cit.>.
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 repulsivly-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=btan(α/2) (see
Figure <ref>e).  
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  

    ℋ^h=∑_m,n∑_l=1^3K2(dr_l-dr_0)^2,

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  

    dr^2=(dx+du)^2+(dy+dv)^2=dr_0^2+2(dxdu+dydv)+du^2+dv^2

where dx=x_i-x_j, dy=y_i-y_j, du=u_i-u_j, dv=v_i-v_j and
dr_0=(dx^2+dy^2)^1/2 is the length at rest of the spring linking
particles i and j, which can take the values c or 2a (see
Figure <ref>e and Tables <ref>,
<ref>).  
Since we consider the expansion around mechanical equilibrium, we
will drop terms linear in du and dv and write:
dr^2=dr_0^2+du^2+dv^2.  
Expanding to harmonic order the expression of dr that we get from
Eq.(<ref>), we obtain  

    dr=dr_0+du^22dr_0(1-dx^2dr_0^2)+dv^22dr_0(1-dy^2dr_0^2)-dxdydudvdr_0^3.

The Hamiltonian of the straight-stripes configuration, with particle
positions specified for the unit-cell in Table <ref>, is

    ℋ_s^h=K2∑_m,n[du_1^2+a^2c^2(du_2^2+du_3^2)+b^2c^2(dv_2^2+dv_3^2)-2abc^2(-du_2dv_2+du_3dv_3)].

Using the relations du_l=u_l-u_0 and dv_l=v_l-v_0 we get

    [                                                                         ℋ_s^h= K∑_m,n[u_0^2-u_0u_1+a^2c^2(2u_0^2-u_0u_2-u_0u_3)+b^2c^2(2v_0^2-v_0v_2-v_0v_3).;                                                                               ;                                                                                                                       .-2abc^2(u_0v_2+u_2v_0-u_0v_3-u_3v_0)]. ]

The Hamiltonian of the bent-stripes configuration, with particle
positions specified for the unit-cell in Table <ref>, is

    [                                                                                                      ℋ_b^h= K2∑_t,n{du_10^2+cos^2α
      du_50^2+sin^2α dv_50^2-sin(2α)
      du_50dv_50+sin^2(α2)
      (du_40^2+du_41^2+du_31^2).;                                                                                                            ;                                                                                                                                           +sin^2(3α2)du_21^2+cos^2(α2)
      (dv_40^2+dv_41^2+dv_31^2)+cos^2(3α2)
      dv_21^2;                                                                                                            ;                                                                                                                                                       .-sinα(-du_40dv_40+du_41dv_41-du_31dv_31)+sin(3α) 
      du_21dv_21}, ]

where 1≤ t≤ L/2. Indeed, the unit-cell of the bent stripe
configuration includes two particles (see Figure <ref>d):
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 du_l0=u_l-u_0, du_l1=u_l-u_1 and
dv_l0=v_l-v_0, dv_l1=v_l-v_1 we get

    [                                                                                     ℋ_b^h=           K2∑_t,n{[2-sin^2α+3sin^2(α2)+sin^2(3α2)](u_0^2+u_1^2)+(sinα-sin2α+sin3α)(u_0v_5.;                                                                                           ;                                                                                                         +u_5v_0)-sinα(u_0v_4+u_4v_0-u_1v_4-u_4v_1+u_1v_3+u_3v_1)-sin3α(u_1v_2+u_2v_1);                                                                                           ;                                                                                            +[1-cos^2α+3cos^2(α2)+cos^2(3α2)](v_0^2+v_1^2)-2sin^2α  
      v_0v_5-2cos^2(α2)(v_0v_4+v_1v_3;                                                                                           ;                                                                                                                                                        . +v_1v_4)-2cos^2(3α2)v_1v_2}. ]


For a harmonic interparticle potential, the Hamiltonian for straight
and bent stripes configurations we expanded around mechanical
equilibrium takes the quadratic form: 
ℋ=K∑_m,nA_m,nq_mq_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 <cit.>: Δ 
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 Ref.<cit.> 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
<ref>). 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 <ref>a), which includes the simultaneous
fluctuation of all particles in the system.  

In Figure <ref>a 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. 

Now we consider the same 2D triangular network model, but for
repulsive harmonic interactions, that is we consider the following
Hamiltonian 

    ℋ^r=∑_m,n∑_l=1^3K2(dr_l-dr_0)^2θ(dr_0-dr_l),

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 <ref>b 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 <ref>c,d). For n=3 we consider the particles 0,
1, 2 and 0, 1, 4 for the straight and bent configurations,
respectively (see Figure <ref>c,d).  
As for the 3D spring network of Ref.<cit.>, 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 <ref>). 
The Hamiltonian of the free particle 0 can be written as
ℋ_0^r=1/2Kρ^2∑_i=1^6g_i(α,γ) where
ρ^2g_i(α,γ) is the contribution to ℋ_0^r
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

    Z_0^r=∫_0^∞∫_0^2πexp[-βK2ρ^2∑_i=1^6g_i(α,γ)]ρ
    dρ dγ=1β
      K∫_0^2πdγ∑_i=1^6g_i(α,γ)=I(α)β
    K

where β=1/(K_BT) is the Boltzmann factor and
I(α)=∫_0^2π[∑_i=1^6g_i(α,γ)]^-1dγ=∫_0^2πf^-1(α,γ)dγ
with f(α,γ)=∑_i=1^6g_i(α,γ).  
In this case we have Δ s=ln(I_s(α)/I_b(α)).

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: 

    [                                                       f_s^r(α,γ)= cos^2γ[θ(γ)θ(π/2-γ)+θ(γ-3π/2)θ(2π-γ)]+sin^2(α/2+γ)[θ(γ)θ(π-α/2-γ);                                                                  ;                                                                               +θ(γ-2π+α/2)θ(2π-γ)]+sin^2(α/2-γ)[θ(γ-α/2)θ(π+α/2-γ)];                                                                  ;                                                                      +cos^2γ[θ(γ-π/2)θ(3π/2-γ)]+sin^2(α/2+γ)[θ(γ-π+α/2)θ(2π-α/2-γ)];                                                                  ;                                                                                      +sin^2(α/2-γ)[θ(γ-π-α/2)θ(2π-γ)+θ(γ)θ(α/2-γ)]. ]

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 <ref>c and
Table <ref>), we can write the function f_s^r 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_s^r(α,γ)=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: 

    [                                                       f_b^r(α,γ)= cos^2γ[θ(γ)θ(π/2-γ)+θ(γ-3π/2)θ(2π-γ)]+sin^2(α/2+γ)[θ(γ)θ(π-α/2-γ);                                                                  ;                                                                              +θ(γ-2π+α/2)θ(2π-γ)]+cos^2(α+γ)[θ(γ-π/2+α)θ(3π/2-α-γ)];                                                                  ;                                                                     +sin^2(3α/2+γ)[θ(γ-π+3α/2)θ(2π-3α/2-γ)]+sin^2(α/2+γ)[θ(γ-π+α/2);                                                                  ;                                                                          ·θ(2π-α/2-γ)]+sin^2(α/2-γ)[θ(γ-π-α/2)θ(2π-γ)+θ(γ)θ(α/2-γ)] ]




In Figure <ref> we show 1/f for hamonic and repulsive
harmonic interactions for straight and bent configurations as a
funciont 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 <ref>.  
From Figure <ref> we can see, particularly for big angles
α for which Δ s takes its larger values (see Figure
<ref>b), that repulsion accentuates the contribution to the 
function I for bent stripes (corresponding in Figure <ref>
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 <cit.> and Ising <cit.>
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
<cit.>. 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 <cit.>, for
which only states visited by the system at the observational time
scale are taken into account, or following the Edwards approach
<cit.>, 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 <cit.>. 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). 

This research was supported by the Israel Science
  Foundation Grant No. 968/16. 



