Soft-particle lattice-gas in 1d: one- and two-component cases

TL;DR

This paper investigates a one-dimensional lattice-gas model of soft particles, analyzing one- and two-component systems using transfer matrix methods, and identifies criteria for different behavioral regimes and thermodynamic catastrophe.

Contribution

It extends the transfer matrix analysis to a 1D soft-particle lattice-gas, including two-component interactions, and confirms a criterion for behavioral classification and catastrophe conditions.

Findings

Existence of a criterion separating behavioral regimes in one-component systems

Application of the criterion to two-component systems for catastrophe prediction

Validation of the criterion within the transfer matrix framework

Abstract

The object of the present article is a 1d lattice-gas system comprised of soft-particles, wherein particles interact only if they occupy the same or a neighboring site, as a simple representation of penetrable particles of soft condensed matter. To represent different scenarios, two different realizations of the lattice model are considered, a one-component and a two-component system, where in the two-component case particles of the same species repel and those of opposite species attract each other. The systems are analyzed entirely within the transfer matrix framework. Special attention is paid to the criterion devised in Ref. [Phys. Rev. E 63, 031206 (2001)], which serves to separate two classes of behavior encountered in a one-component penetrable particle systems. In addition to confirm the existence of a similar criterion for the one-component lattice-gas model, we find that the same criterion can be applied to the two-component system to provide conditions for thermodynamic catastrophe.

Equations

$$Z=\frac{L!}{N!(L-N)!}.$$

$$Z=\frac{(N+L-1)!}{N!(L-1)!},$$

$$\beta P=\log\big{(}1+\rho\big{)},$$

$$Z=\frac{L^N}{N!}.$$

$$\beta P=\rho .$$

$$\Xi (\beta \mu ,L)=\sum _{N=0}^{\infty }{e}^{\beta \mu N}Z(N,L),$$

$$\Xi=\sum_{n_{1}=0}^{\infty}\dots\sum_{n_{L}=0}^{\infty}e^{\beta\mu n_{1}}\dots e^{\beta\mu n_{L}}=\bigg{(}\frac{1}{1-e^{\beta\mu}}\bigg{)}^{L}.$$

$$\Xi =\sum _{n_1=0}^{\infty }\cdots \sum _{n_L=0}^{\infty }\frac{e^{\beta \mu n_1}}{n_1!}\dots \frac{e^{\beta \mu n_L}}{n_L!}=e^{L{e}^{\beta \mu }}.$$

$$p(n)=\frac{1}{1+\rho}\bigg{(}\frac{\rho}{1+\rho}\bigg{)}^{n}.$$

$$p(n)=\frac{N!}{n!(N-n)!}\frac{(L-1{)}^{N-n}}{L^N},$$

$$p(n)=\frac{{\rho }^n{e}^{-\rho }}{n!}.$$

$$H(n_1,\dots ,n_L)=\frac{K}{2}\sum _{i=1}^L{n}_i(n_i-1)+\alpha K\sum _{i=1}^L{n}_i{n}_{i+1},$$

$${\Xi }_a=\sum _{n_1=0}^{\infty }\cdots \sum _{n_L=0}^{\infty }{e}^{\beta \mu n_1}\dots e^{\beta \mu n_L}{e}^{-\beta H(n_1,\dots ,n_L)},$$

$${\Xi }_b=\sum _{n_1=0}^{\infty }\cdots \sum _{n_L=0}^{\infty }\frac{e^{\beta \mu n_1}\dots e^{\beta \mu n_L}}{n_1!\dots n_L!}{e}^{-\beta H(n_1,\dots ,n_L)},$$

$$\Xi =\sum _{n_1=0}^{\infty }\cdots \sum _{n_L=0}^{\infty }T(n_1,n_2)T(n_2,n_3)\dots T(n_L,n_1),$$

$$\Xi =\mathrm{Tr}{\mathbf{T}}^L.$$

$$\Xi =\sum _{i=1}^{\infty }{\lambda }_i^L,$$

$$\Xi ={\lambda }_1^L,$$

$$\beta P=\log {\lambda }_1.$$

$$p(n)=\frac{1}{\Xi }\sum _{n_2=0}^{\infty }\cdots \sum _{n_L=0}^{\infty }T(n,n_2)T(n_2,n_3)\dots T(n_L,n),$$

$$p(n)={\varphi }_1^2(n),$$

$$T_a(n,n^{\prime })=e^{-\frac{\beta K}{4}(n^2+n^{\prime 2})}{e}^{-\beta \alpha Kn{n}^{\prime }}{e}^{\frac{\beta {\mu }^{\prime }}{2}(n+n^{\prime })},$$

$$T_b(n,n^{\prime })=e^{-\frac{\beta K}{4}(n^2+n^{\prime 2})}{e}^{-\beta \alpha Kn{n}^{\prime }}\frac{e^{\frac{\beta {\mu }^{\prime }}{2}(n+n^{\prime })}}{\sqrt{n!n^{\prime }!}}.$$

$$p_m(n,n^{\prime })$$

$$\displaystyle\frac{p_{m}(n,n^{\prime})}{p(n)p(n^{\prime})}=1+\sum_{k=2}^{\infty}\bigg{(}\frac{\lambda_{k}}{\lambda_{1}}\bigg{)}^{m}\frac{\phi_{k}(n)\phi_{k}(n^{\prime})}{\phi_{1}(n)\phi_{1}(n^{\prime})},$$

$$\sum _{n=0}^{\infty }\sum _{n^{\prime }=0}^{\infty }{p}_m(n,n^{\prime })$$

$${\Gamma }_m=\frac{p_m(0,0)}{p^2(0)}-1.$$

$$\begin{array}[]{r l}\Gamma_{m}>0,&\quad\text{for $m$-even}\\ \Gamma_{m}<0,&\quad\text{for $m$-odd}.\end{array}$$

$$\begin{array}[]{r l}\lambda_{n}>0,&\quad\text{for $n$-odd}\\ \lambda_{n}<0,&\quad\text{for $n$-even},\end{array}$$

$$\Gamma_{m}\approx\bigg{(}\frac{\lambda_{2}}{\lambda_{1}}\bigg{)}^{m}\bigg{(}\frac{\phi_{2}(0)}{\phi_{1}(0)}\bigg{)}^{2}.$$

Full text

I Introduction

The present work investigates a 1d lattice model of soft particles. Soft interactions imply the possibility of multiple occupations of a single site, and the model is intended to be a simple representation of penetrable particles, such as the Gaussian core or the Penetrable-Sphere model.

In soft condensed matter, penetrable pair potentials have been recognized as realistic representations of effective interactions between a number of large macromolecules. For example, the Penetrable-Sphere model represents micelles in a solvent Witten89 , while the Gaussian core model accurately captures entropic repulsions between self-avoiding polymer coils in a good solvent Hansen00 ; Hansen01 . The the generalized exponential model, $\exp(-r^{4})$, accurately describes the effective repulsion between flexible dendrimers in a solution Mladek07 ; Mladek08a ; Mladek08b ; Mladek10 . Other soft potentials have been suggested to include an ever wider class of macromolecules Likos01a . Penetrability has been further extended to include charged macromolecules, represented by a divergence-free Coulomb potential Hansen12 ; Frydel13 ; Frydel16 .

Due to the absence of hard-core interactions, the density of penetrable systems may acquire arbitrarily large values. This in turn leads to new behaviors not seen in standard models. For example, a solid phase of Gaussian particles melts in two ways: when pressure is reduced, as in normal melting, and when pressure is increased, known as reentrant melting Saija05 ; Saija05b ; Ikeda13 .

Penetrable-spheres, on the other hand, do not undergo reentrant melting. A solid phase is preserved all the way into infinite densities. A lattice structure, instead of being shrunk or deformed upon compression in order to create additional sites, remains intact and excess of particles is accommodated by allowing multiple occupation of existing lattice sites. Such “stacks” of several particles sharing the same space can form already in a dense liquid phase prior to freezing. Advantage of a “stacked” over a more uniformly distributed system is that the “stacking” arrangement reduces a number of overlaps by minimizing interactions between particles in different “stacks”. This, in turn, lowers the system energy.

Between the Gaussian core and the Penetrable-Sphere model there exists a continuum of functional forms of penetrable potentials that fall into one of the two classes of behavior. Using the mean-field analysis, Likos et al Likos01 determined criterion for predicting a type of behavior for any given pair potential. If the Fourier transform of the pair interaction is everywhere positive, then one expects the Gaussian model like behavior. If, on the other hand, it is not positive everywhere, the system exhibits “stacking” formations as in the Penetrable-Sphere model.

The leading motivation for the present work is to shed light on the two classes of behavior and to better understand the Likos-Lang-Watzlawek-Löwen (LLWL) criterion in the context of a simple lattice-gas system. The model consists of a 1d array of discrete sites. There is no bound on how many particles may occupy a single site. Particles interact only if they are on the same or a neighboring site. The interaction strength is regulated with two parameters, $K$ for interactions between particles on the same site, and $K^{\prime}=\alpha K$, for interactions between particles occupying neighboring sites.

To analyze the model and its properties, we use the transfer matrix method. Within this methodology a given system is characterized by a transfer matrix. All the thermodynamic quantities can then be expressed in terms of the transfer matrix eigenvalues and eigenvectors. As the occupation number is unlimited, the transfer matrix has infinite size. In practice, however, a $20\times 20$ matrix suffices for most situations. Eigenvalues and eigenvectors of the transfer matrix are then calculated numerically using any of the standard software packages, such as Mathematica or Matlab.

In addition to a one-component scenario, the work considers a two-component system, where particles of the same species repel and those of opposite species attract each other. Such systems have been considered for penetrable-spheres Frydel17 and the Gaussian core model Frydel18 . As the two-component Gaussian core model is well behaved, the two-component Penetrable-Sphere model is thermodynamically unstable. We show that the LLWL criterion of classification Likos01 can be extended and applies to a two-component system, where it serves as a criterion of thermodynamic stability.

The work is organized as follows. In Sec. (II) we consider a 1d lattice model for non-interactive particles. Different ways of counting configurations are considered, leading to different partition functions. In Sec. (III) we consider a one-component lattice model with soft repulsive interactions. Here we determine the existence of two types of behavior, in agreement with penetrable particles of soft condensed matter. In Sec. (IV) we consider a two-component system with particles of the same species repelling and of different species attracting each other. We show that the LLWL criterion of a one-component case apply to the two-component system as a criterion of thermodynamic stability. Finally, in Sec. (V) we conclude the work.

II non-interactive particles

Given a 1d array of $L$ lattice sites and $N$ indistinguishable and non-interactive particles, the canonical partition function, for the case where at most one particle can occupy a lattice site, is

[TABLE]

and corresponds to a binomial coefficient $C(L,N)$. In this work, however, we are interested in a lattice model with multiply occupied sites. The partition function for this situation is

[TABLE]

and corresponds to the binomial coefficient $C(N+L-1,N)$ and represents the permutation formula for $L-1$ items of type one and $N$ items of type two. If $L-1$ items are assumed to represent bars, then these bars segregate $N$ items into $L$ sets, where a single set represents a lattice site. Interpreted in this way, Eq. (2) makes perfect sense.

Note, however, that the pressure per lattice site for the above system, defined as $P=\frac{\log Z}{\partial L}$, is

[TABLE]

and does not correspond to the ideal-gas behavior.

To construct the partition function that reproduces ideal-gas properties, we must proceed from the assumption that particles are distinguishable, in which case there are $L^{N}$ distinct configurations. The standard trick to obtain the partition function is to use the Gibbs correction, $1/N!$, yielding

[TABLE]

The resulting partition function now yields the correct ideal-gas pressure per lattice site,

[TABLE]

Difference between the partition function in Eq. (2) and that in Eq. (4) is well illustrated with different simulation algorithms. One algorithm generates configurations by randomly selecting a lattice site. This is followed by either adding or subtracting a particle. At the end of an update cycle, consisting of $L$ random picks of a lattice site, one ensures that the total number of particles in a system is conserved. This algorithm corresponds to $Z$ in Eq. (2).

In an alternative algorithm, configurations are generated by randomly selecting a particle (not a site), hence, particles are labeled. A selected particle is then moved to a randomly selected site. This algorithm corresponds to $Z$ in Eq. (4) and is more suitable for representing liquids.

As it is more convenient to work with grand partition functions, below we obtain the corresponding expressions. The formal relation between the grand and the canonical partition function is

[TABLE]

where $\mu$ is the chemical potential and $\beta=1/k_{B}T$. For indistinguishable particles, using Eq. (1), $\Xi=\sum_{N=0}^{\infty}e^{\beta\mu N}\frac{(N+L-1)!}{N!(L-1)!}=(1-e^{\beta\mu})^{-L}$, where only $\mu<0$ is physically meaningful. At $\mu=0$ the partition function diverges and for $\mu<0$ it becomes negative. The same result is obtained from an alternative formulation

[TABLE]

The physical interpretation of the above expression is clear. $L$ summations correspond to $L$ sites and $n_{i}$ designates the number of particles at a site $i$. In the grand ensemble an average number of particles at each site is controlled with $\mu$.

For distinguishable particles the grand partition function we use Eq. (2), leading to $\Xi=\sum_{N=0}^{\infty}e^{\beta\mu N}L^{N}/N!=(e^{e^{\beta\mu}})^{L}$. The same result follows from an alternative formulation

[TABLE]

Later in this work we consider probabilities $p(n)$, the probability that any given site is occupied by $n$ particles. Consequently, we derive these probabilities for non-interactive particles. The procedures to obtain $p(n)$ will furthermore clarify the difference between distinguishable and indistinguishable particles.

We start by recalling that $Z(N,L)$ in Eq. (2) counts the number of configurations for $N$ indistinguishable particles distributed over $L$ sites. If one site is occupied by $n$ particles, the number of configurations of the remaining $N-n$ particles distributed over $L-1$ sites is $Z(N-n,L-1)$, and $p(n)$ is given by the ratio of the two numbers, $p(n)=\frac{Z(N-n,L-1)}{Z(N,L)}$. In the thermodynamic limit, $L\to\infty$ and $N\to\infty$ such that $N/L=\rho$ (and using the Sterling formula $N!\approx N^{N}e^{-N}$ and the limiting representation of the exponential function $e^{x}=\lim_{N\to\infty}(1+x/N)^{N}$), that expression reduces to

[TABLE]

For distinguishable particles the total number of configurations is $L^{N}$. If one site is occupied by $n$ labeled particles, the number of configurations of the remaining $N-n$ particles distributed over $L-1$ sites becomes $(L-1)^{N-n}$. However, the ratio $(L-1)^{N-n}/L^{N}$ does not yield the probability $p(n)$. Since there are $\frac{N!}{n!(N-n)!}$ different ways to draw $n$ labeled particles, the ratio $(L-1)^{N-n}/L^{N}$ needs to be multiplied by that number. The correct distribution becomes

[TABLE]

which in the thermodynamic limit recovers the Poisson distribution

[TABLE]

III Interactions: one-component system

We next consider an interactive 1d lattice system represented by the Hamiltonian

[TABLE]

where the interactions between particles on the same site are given by the first term, and the interactions between particles on neighboring sites by the second term and are regulated by the dimensionless parameter $\alpha$. For $\alpha=0$ interactions between particles on neighboring sites are turned off, and for $\alpha=1$ these interactions are the same as those for particles on the same site. The case $\alpha=1$ can be regarded as analogous to the Penetrable-Sphere model, and the case $0<\alpha<1$ to the Gaussian core model. We are not interested in the scenario $\alpha>1$ which has no correspondence in actual penetrable particle systems and implies that interactions between particles on neighboring sites are greater than those for particles on the same site. Finally, the scenario $\alpha<0$ implies attraction between particles on neighboring sites, in possible analogy to the van der Waals type of potentials, however, in this work we do not pursue this case.

As for non-interactive particles, we consider the system of indistinguishable and distinguishable particles. For indistinguishable particles the grand partition function is

[TABLE]

and for distinguishable particles it is

[TABLE]

where we use the index “a” and “b” to differentiate between the two cases. Both cases adapt periodic boundary conditions, $n_{L+1}=n_{1}$, which ensures that each site is equivalent.

The systems are analyzed using the transfer matrix method, the standard method for dealing with lattice models in 1d Lavis . The central object of the method is the transfer matrix, $T(n,n^{\prime})$, by means of which the partition function can be written as

[TABLE]

revealing chain structure of a partition function. Using matrix algebra, the partition function is shorthanded into

[TABLE]

Eigendecomposition of the transfer matrix, ${\bf T}={\bf Q}{\bf\Lambda}{\bf Q}^{T}$ (where ${\bf\Lambda}$ is the diagonal matrix with diagonal elements $\Lambda_{ii}=\lambda_{i}$ and ${\bf Q}$ is the square matrix whose $i$-th column is the eigenvector $\phi_{i}(n)$ of ${\bf T}$), further transforms the expression into

[TABLE]

where $\lambda_{i}^{L}$ are eigenvalues of the matrix ${\bf T}^{L}$, and ${\bf T}^{L}$ is the product matrix generated by multiplying ${\bf T}$ by itself $L$-times. If eigenvalues are ordered according to their modulus as $|\lambda_{1}|>|\lambda_{2}|>|\lambda_{3}|\dots$, and because in the thermodynamic limit, $L\to\infty$, $\Xi$ is dominated by the largest eigenvalue, the grand partition function simply becomes

[TABLE]

and the corresponding pressure per lattice site is given by

[TABLE]

We next use the transfer matrix framework to obtain the probability $p(n)$, that a given site $i$ is occupied by $n$ particles. The formal definition is

[TABLE]

and amounts to breaking the ring structure of Eq. (14) at a site $i=1$, giving rise to a linear chain. After the application of eigendecomposition (see Appendix A), the expression reduces to

[TABLE]

where $\phi_{1}(n)$ are elements of the dominant eigenvector corresponding to the largest eigenvalue $\lambda_{1}$, thus, $p(n)$ is properly normalized since the modulus of a vector $\phi_{1}(n)$ is $1$.

The transfer matrix for indistinguishable particles, corresponding to the partition function in Eq. (12), is

[TABLE]

and that for distinguishable particles corresponding to $Z$ in Eq. (13) is

[TABLE]

Eignevalues and the eigenvectors are calculated numerically using Mathematica. In principle, $T(n,n^{\prime})$ is an infinite matrix, but in practice the $20\times 20$ matrix is sufficient for most situations.

Fig. (1) shows a number of distributions $p(n)$ for indistinguishable particles, for $\beta K=0.1$ and $\alpha=1$, for different densities arranged in increasing order. The consecutive plots show gradual transformation of $p(n)$ into a bimodal structure, emerging at around $\rho=4$. The two peaks of the bimodal structure are at $n=0$ and $n\approx 2\rho$, suggesting an alternating structure of occupied versus empty sites, rather than a coexistence of vacuum cavities embedded in a fluid with density $2\rho$. The emergence of an alternating structure is analogous to the “stack” formations of penetrable-spheres discussed in the introduction.

A similar transformation into a bimodal structure occurs for distinguishable particles, see Fig. (2). The crossover, however, occurs at a higher density, $\rho\approx 10$. The explanation for this difference lies in different entropies of the two systems. For the case of distinguishable particles the adaptation of an ordered alternating structure entails larger loss of entropy.

To confirm the existence of an alternating structure, we consider the two-site probability, $p_{m}(n,n^{\prime})$, which is the probability that two sites separated by $m$ sites have the occupation number $n$ and $n^{\prime}$, and whose formal definition is

[TABLE]

Applying eigendecomposition (see Appendix B), we arrive at the expression

[TABLE]

where $p_{m}(n,n^{\prime})$ can be shown to be normalized,

[TABLE]

where the second term vanishes for any $k\neq 1$ as the consequence of orthonormality of the eigenvectors $\phi_{k}$.

We next define the quantity

[TABLE]

In absence of correlations between empty sites, $\Gamma_{m}=0$. On the other hand, if occupied and empty sites alternate, we expect

[TABLE]

In Fig. (3) we plot $|\lambda_{n}|$ for the density when the distribution starts to separate into bimodal structure and for the density where the distribution has a well developed bimodal structure. The eigenvalues alternate in sign as

[TABLE]

which is not captured by the figure which plots the data points for $|\lambda_{n}|$. The main result is that once the bimodal structure is established, the spectrum is dominated by the first two eigenvalues, $\lambda_{1}$ and $\lambda_{2}$, so that $\Gamma_{m}$ in Eq. (26) can be approximated as

[TABLE]

Furthermore, we find that $\phi_{2}(0)/\phi_{1}(0)\approx 1$, so that the correlations are determined solely by the ratio $\lambda_{2}/\lambda_{1}<0$. Since $\lambda_{1}$ and $\lambda_{2}$ have different sign, $\lambda_{2}/\lambda_{1}$ raised to odd power is negative, and raised to even power it is positive. This suggests

[TABLE]

Fig. (3) indicates that the first two eigenvalues converge, $|\lambda_{2}|\to|\lambda_{1}|$, but attain equality only in the limit $\rho\to\infty$. Since the condition $|\lambda_{2}|=|\lambda_{1}|$ implies long-range order, the phase transition in 1d does not take place at finite density Cuesta04 . Only for dimensionality $d>1$ a phase transition is possible.

$\Gamma_{m}$ plotted in Fig. (4) confirms the exponentially decaying correlations and accuracy of the ansatz in Eq. (30).

So far we have considered only the case $\alpha=1$, which can be regarded as representative of the Penetrable-Sphere model whose interaction strength remains the same as long as spheres are overlaped. In the subsequent section we look into other values of $\alpha$, especially, we examine the effect of reduced $\alpha$ on the structure of $p(n)$.

In Fig. (5) we plot several distributions $p(n)$, for $\beta K=0.1$ and fixed $\rho$, for decreasing values of $\alpha$. The results indicate gradual transformation of a bimodal into a mono-modal structure, implying the dissolution of an ordered alternating structre. For indistinguishable particles at $\alpha=0.55$ the bimodal structure is no longer there, however, it reappears if density is increased. The question is, what is the critical value of $\alpha$ below which the ordered alternating structure never arises for any density?

To answer this, we consider two idealized configurations. One configuration is uniform, with each site having the occupation number $n_{i}=\rho$. Another one have alternating occupations between $n_{i}=2\rho$ and $n_{i}=0$. From the Hamiltonian in Eq. (11), energy of each configuration is

[TABLE]

and

[TABLE]

and the energy gained by each site by replacing a homogenous structure with an alternating structure is

[TABLE]

The results indicate that for $\alpha\leq 0.5$ there is no longer an energy gain by adapting an ordered alternating structure. The condition $\alpha>0.5$, however, is sufficient but not necessary to have an alternating structure, as this condition was determined from energy considerations alone. In real systems the energy gain, by switching to more ordered structure, is accompanied by loss of entropy. In addition to the condition $\alpha>0.5$, we need the condition where energy is a dominant contribution of the free energy. This is attained for large interaction strength and large density.

To illustrate how the structural rearrangement of adapting an alternating structure affects thermodynamic quantities of a system, in Fig. (6) we plot pressure per lattice site as a function of $\alpha$. Initially the pressure increases linearly with $\alpha$. At the crossover, when the distribution $p(n)$ becomes bimodal, this trend abruptly changes and pressure begins to decrease, reflecting the structural rearrangement of a system that is concurrent with the reduction of internal tensions.

So far we have examined strictly weak interactions, and all the results for $\beta K=0.1$, where the transformation into ordered alternating structure requires large densities.

Additional feature of systems with strong interactions is that once an ordered alternative structure is attained, it is followed by additional regular transformations upon further squeezing of the system at around densities corresponding to integer values, that is, when the occupation of a lattice site changes from $n\to n+1$ ($n$ being an integer). This results in a steplike structure of the pressure isotherm, which disappears at weak interactions. Such a strong interaction case and the resulting regular transformations were carefully studied in Prestipino14 ; Prestipino15a . It was demonstrated by the researchers of the above references that the transformations at $n\to n+1$ correspond to sharp crossovers rather than representing a true phase transitions.

Limiting ourselves to weak interactions in the present work, we eliminate the complication of the steplike behavior and demonstrate the onset of an ordered alternating structure through the continuous-like shape of $p(n)$.

IV Two-component system

In this section we consider a two-component 1d lattice model, where particles of the same species repel, and those of the opposite species attract one another. The system Hamiltonian is

[TABLE]

where the two species are labeled as “$+$” and “$-$”, in analogy to charged systems. The first line is for the interactions between particles occupying the same site (the second term on that line subtracts self-interaction introduced in the first term), and the second line is for particle interactions occupying neighboring sites. We rewrite Eq. (34) as

[TABLE]

so that, apart for one term that in partition function is incorporated into the chemical potential, it is written in terms of ”charge” per lattice site, $s_{i}=n_{i}^{+}-n_{i}^{-}$.

As for the one-component system, we consider both indistinguishable and distinguishable particles. The partition function for the indistinguishable case is

[TABLE]

and that for the distinguishable one

[TABLE]

Both partition functions can be transformed into the summations over $s_{i}=n_{i}^{+}-n_{i}^{-}$ (see Appendix C),

[TABLE]

and

[TABLE]

where $\mu^{\prime}=\mu+K/2$ and $I_{n}(x)$ in the second equation is the modified Bessel function of the first kind,

[TABLE]

The corresponding transfer matrices are

[TABLE]

and

[TABLE]

Note that the matrices are similar to those in Eqs. (21,22) for the one-component system. What makes the two systems different is that the indices $s$ and $s^{\prime}$ in Eqs. (41,42) are for all integers, raising the possibility of the term $e^{-\beta\alpha Kss^{\prime}}$ to dominate the transfer matrix if $s$ and $s^{\prime}$ have opposite sign and to eventual divergence of the partition function. The divergence, however, can be switched off for sufficiently low $\alpha$.

For example, if we consider elements of the transfer matrix corresponding to $s^{\prime}=-s$,

[TABLE]

we find that the possibility of divergence in the limit $s\to\infty$ is prevented if $\alpha\leq 0.5$ (keeping in mind that $\mu^{\prime}<0$). For $\alpha>0.5$, no matter how negative $\mu^{\prime}$, the divergence can never be suppressed. If we consider in turn the case $s>0$ and $s^{\prime}=-1$,

[TABLE]

we discover that the divergence in the limit $s\to\infty$ never arises for any $\alpha$, as the expression is dominated by $e^{-\frac{\beta K}{4}s^{2}}$, which vanishes in the same limit. We conclude that divergent elements in the limit $s\to\infty$ are those that roughly satisfy $s^{\prime}\approx-s$.

The presence of the divergent terms implies that a charge, $\langle s\rangle$, at a single site becomes infinite, which, by the same token, implies that the occupation number at a single site diverges, pointing out to thermodynamic instability (also known as catastrophe).

To see how the presence of the divergence plays its role, in Fig. (7) we plot $|\lambda_{n}|$ for $\alpha=1$ as a function of $n$, for different values of $M$, where $M$ is the size of a transfer matrix $M\times M$. In the stable system, $|\lambda_{n}|$ converges in the limit $M\to\infty$. In the unstable system, $|\lambda_{n}|$ diverge in the same limit. In addition to blowing up of $|\lambda_{n}|$ as $M$ increases, we observe increased domination of the two initial eigenvalues, $|\lambda_{1}|\approx|\lambda_{2}|$, similar to what was seen for the one-component system in Fig. (3), and which indicates an ordered alternative structure. For the two-component system, the ordering implies an alternating occupation of each site by a species ”$+$” and ”$-$”.

For comparison, in Fig. (8) we plot $|\lambda_{n}|$ for a stable system at $\alpha=0.5$. The results indicate that the two leading eigenvalues are separated even in the limit of high $\rho$, therefore, never come to dominate. The stability in this case implies the absence of an ordered alternating structure, and the region of stability corresponds to $\alpha\leq 0.5$.

In search of possible structures in a stable system at $\alpha=0.5$, in Fig. (9) we plot the distributions $p(s)=\phi_{1}^{2}(s)$ for $\alpha=0$ and $\alpha=0.5$. For $\alpha=0.5$ the distributions appear broader but otherwise fail to develop a bimodal structure. Indistinguishable particles show less response to the variation with $\alpha$, which we attribute to a higher entropy cost in adapting structured configuration.

The most interesting structural feature of the two-component system at $\alpha=0.5$ is the formation of semi-stable pairs between particles of opposite species. In electrolytes these pairs are referred to as the Bjerrum pairs Yan93 . The presence of such pairs is evident in “charge” fluctuations, $\langle s^{2}\rangle$, plotted in Fig. (10) as a function of $\beta K$ for $\rho=10$. For large values of particle interactions the fluctuations are suppressed, indicating that particles interchange sites not as free particles but as as permanent pairs.

IV.1 significance of thermodynamic instability for real systems

As stated above, the two-component lattice model is thermodynamically unstable for $\alpha>0.5$. This criterion was derived from the grand canonical ensemble. To determine how this instability is manifested in real systems, we perform a sequence Monte Carlo simulations in a canonical ensemble for a two-component 1d lattice model for $\alpha>0.5$. The system is finite with periodic boundary conditions. The system size is $L=1000$, there are $N^{+}=N^{-}=10000$ particles, and the interaction strength is $\beta K=0.1$. All simulations start with randomly distributed particles.

Configuration snapshots in Fig. (11) reveal that a system collapses into a finite number of clusters, the so called thermodynamic catastrophe Ruelle66a ; Ruelle66b . As all the sites are equivalent, the clusters form by spontaneous nucleation. For $\alpha=1$ in Fig. (11) (a) particles are distributed over a several five-site clusters. As $\alpha$ decreases, larger clusters are preferred. A snapshot shown in Fig. (11) (b) for $\alpha=0.55$ indicates that an entire system exists as a single sixteen-site cluster. The cluster disintegrates for $\alpha\leq 0.54$, as there are not enough particles to reorganize into a larger cluster. This does not imply thermodynamic stability, however. Changing the system size would produce another collapse into larger clusters.

To understand the dependence of the cluster size on $\alpha$, we consider a number of clusters. The smallest possible cluster is comprised of two sites. The energy of this cluster is obtained from the Hamiltonian in Eq. (35). Assuming that two sites have the same occupation number $n/2$, where $n$ is the total number of particles in the cluster, the energy is

[TABLE]

where the first term, which is proportional to $n^{2}$, is positive if $\alpha<1$, therefore, the system does not collapse into a single cluster for $\alpha\leq 1$, as the cluster energy eventually becomes positive for large $n$.

Next we consider the three-site cluster. If the total number of particles in the cluster is $n$, and the occupation of the central site is $n_{0}$ and that of each flanking site is $n-n_{0}$ (we assume symmetricity of a cluster), then based on Eq. (35) the energy is

[TABLE]

This energy, optimized with respect to $n_{0}$, becomes

[TABLE]

At $\alpha=1/\sqrt{2}\approx 0.707$ the term proportional to $n^{2}$ becomes positive and the three-site cluster disintegrates as its energy becomes positive for large $n$.

We can repeat the same calculations for a four-site cluster, whose optimized energy is (optimized clusters are symmetric)

[TABLE]

and the cluster disintegrates at $\alpha=(\sqrt{5}-1)/2\approx 0.618$, as the first term becomes positive. Each consecutive cluster disintegrating at $\alpha$ closer to $0.5$. For example, a five-site cluster disintegrates at $\alpha=1/\sqrt{3}\approx 0.577$.

The observed thermodynamic catastrophe is not an artifact of a lattice model, and a similar behavior has been observed for a two-component penetrable-spheres Frydel17 . A snapshot for a 2d two-component Penetrable-Sphere model is shown in Fig. (12), which reveals a similar catastrophe, characterized by a system collapse into small number of large clusters.

On the other hand, thermodynamic catastrophe does not occur for a two-component Gaussian core model Frydel18 , indicating that this system does not fulfill the conditions of instability.

In the one-component system the criterion $\alpha>0.5$ tells us that a system under certain conditions of density and interaction strength can adopt an ordered alternating structure. The condition $\alpha>0.5$, therefore, is necessary but not sufficient for this situation. In the two-component system, on the other hand, the criterion $\alpha>0.5$ tells us that a system is thermodynamically unstable under any conditions, no matter what its density and the interaction strength, as long as a system is in thermodynamic limit. The criterion $\alpha>0.5$ for the two-component system, therefore, is necessary and sufficient.

IV.2 criterion for thermodynamic catastrophe

Thermodynamic catastrophe of a one-component system was first investigated by Ruelle and Fisher in Ruelle66a ; Ruelle66b . The condition for thermodynamic instability for these systems is $\int d{\bf r}\,u(r)<0$, or using the Fourier transformed pair potential $\tilde{u}(k)$, the same condition is stated as $\tilde{u}(0)<0$. This implies that the potential needs to have an attractive part and a non-divergent (soft) core. An example of such a potential is a double Gaussian potential investigated in Prestipino15b ; Prestipino16 .

Based on our results for a lattice model and its connection to real penetrable particles, we conclude that the condition for thermodynamic instability for a two-component system are provided by the LLWL criterion. That is, if $\tilde{u}(k)$ is negative for some value of $k$, then the two-component system is thermodynamically unstable. In turn, if $\tilde{u}(k)$ is positive everywhere, then the system is stable. This justifies why thermodynamic catastrophe is observed for the two-component Penetrable-Sphere but not the Gaussian core model.

V Conclusion

The present work investigates a 1d lattice model with multiple occupations as a simple representation of penetrable particles. Starting with non-interacting particles, we distinguish between different ways of counting configurations, by treating particles as either indistinguishable or distinguishable, leading to two different partition functions. The indistinguishable case is representative of growth models, and the distinguishable case is representative of liquids.

For a one-component case we discover two classes of behavior, depending on weather $\alpha>0.5$ or $\alpha\leq 0.5$. For $\alpha>0.5$, under the conditions of large $\rho$ and/or large $\beta K$, systems form an ordered and alternating structure of occupied versus empty sites. For $\alpha\leq 0.5$ such structural reorganization never take place. Because the condition $\alpha<0.5$ does not guarantee the presence of an alternating structure, we say that this condition is necessary but not sufficient.

For the two-component system, where particles of the same species repel and of opposite species attract each other, we find that systems with $\alpha>0.5$ are thermodynamically unstable, and such systems exhibit thermodynamic catastrophe Ruelle66a ; Ruelle66b . In this case, the criterion $\alpha>0.5$ is necessary and sufficient for thermodynamic instability, as long as the system is in thermodynamic limit. Even a dilute system with weak interactions eventually collapses.

Thermodynamic instability was observed in the Penetrable-Sphere model, but not in the Gaussian core model. Consequently, we conclude that the LLWL criterion devised for a one-component system apply to a two-component system for predicting the conditions of thermodynamic instability.

A shortcoming of a 1d lattice model is that it does not undergo a true phase transition. A more realistic representation of penetrable particles would require working in 2d or higher dimension, where such transition becomes feasible. Some aspects of a lattice model in higher dimension studied in Hansen04 for repulsive interactions. In future we plan to investigate a two-component lattice model in 2d or a higher dimension, in order to study a gas-liquid phase transition and the role of Bjerrum pairs in the transition mechanism Frydel18 .

Appendix A Derivation of $p(n)$ within the transfer matrix method

In Sec. (III) in Eq. (19) we provide the relation $p(n)=\phi_{1}^{2}(n)$ of the transfer matrix method, for the probability $p(n)$ that a single lattice site is occupied by $n$ particles. The relation can be rigorously derived from the formal definition

[TABLE]

although the derivational steps are omitted. The aim of this appendix is to provide the omitted derivation. Note that $p(n)$ is normalized by construction, $\sum_{n=0}^{\infty}p(n)=1$.

The first step is to rewrite the definition in Eq. (49) in shorthanded form, using the matrix algebra nomenclature that is more convenient for carrying out subsequent matrix operations,

[TABLE]

where $({\bf T}^{L})(n,n^{\prime})$ designates the $(n,n^{\prime})$ element of the matrix ${\bf T}^{L}$, where ${\bf T}^{L}$ is the product matrix generated by multiplying ${\bf T}$ by itself $L$-times. In the next step we apply the eigendecomposition

[TABLE]

where ${\bf Q}$ is the square matrix whose $k$-th column is the eigenvector $\phi_{k}(n)$ of the transfer matrix ${\bf T}$, ${\bf Q}^{-1}$ is the inverse of ${\bf Q}$ such that

[TABLE]

where ${\bf I}$ is the identity matrix. Eigendecomposition applied to ${\bf T}^{L}$ yields

[TABLE]

The probability $p(n)$ in Eq. (50) can now be written as

[TABLE]

Because the transfer matrix is real and symmetric, ${\bf Q}^{-1}={\bf Q}^{T}$, where ${\bf Q}^{T}$ is the transpose of ${\bf Q}$, and we write

[TABLE]

where we used $\Xi={\rm Tr}\,{\bf T}^{L}=\sum_{k=1}^{\infty}\lambda_{k}^{L}$. Since the columns of the matrix ${\bf Q}$ correspond to eigenvectors $\phi_{k}$, we get

[TABLE]

In the final step we take the thermodynamic limit, $L\to\infty$, in which the above expression reduces to

[TABLE]

which recovers the result of Eq. (19).

Appendix B Derivation of $p_{m}(n,n^{\prime})$ within the transfer matrix method

In this appendix we derive the expression for the probability $p_{m}(n,n^{\prime})$, that a number of particles at two lattice sites separated by $m$ sites is $n$ and $n^{\prime}$. The derived result appears in Eq. (24) but no derivation is provided. We start with the formal definition for $p_{m}(n,n^{\prime})$,

[TABLE]

Note that $p_{m}(n,n^{\prime})$ is normalized by construction, $\sum_{n=0}^{\infty}\sum_{n^{\prime}=0}^{\infty}p_{m}(n,n^{\prime})=1$. Using matrix algebra nomenclature, the above expression can be shorthanded into

[TABLE]

Then eigendecomposition yields

[TABLE]

and

[TABLE]

leading to

[TABLE]

Further simplifications follow from the fact that for a real and symmetric matrix ${\bf T}$, ${\bf Q}^{-1}={\bf Q}^{T}$, and the columns of the matrix ${\bf Q}$ correspond to eigenvectors $\phi_{k}$. This leads to

[TABLE]

where we used $\Xi=\sum_{k=1}^{\infty}\lambda_{k}^{L}$. The final reduction comes from taking the thermodynamic limit, $L\to\infty$, in which case only $k^{\prime}=1$ does not vanish, leading to the final result

[TABLE]

Using the relation $p(n)=\phi_{1}^{2}(n)$ of the previous section we get

[TABLE]

which agrees with Eq. (24).

Appendix C Reduction of the partition functions for the two-component system

We reduce the partition functions of the two-component systems, from the $2L$-summation to the $L$-summation, by considering the simple case $L=3$, however, the procedure is general and valid for any $L$. For distinguishable particles the partition function for three lattice sites, obtained from Eq. (36), is

[TABLE]

where for the sake of clarity we write down every term explicitly. We also recall that $\mu^{\prime}=\mu+K/2$. The above expression can be shortened by using $s_{i}=n_{i}^{+}-n_{i}^{-}$,

[TABLE]

which by itself does not yet transform the summation. In order to transform the above result into summation in terms of $s_{i}$ we have to carry out partial summations of some of the terms, leading to

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function. To complete the transformation we need to calculate

[TABLE]

where the solution is found to be

[TABLE]

and where to obtain it we summed up all the terms corresponding to a given $s_{i}$. For example, for $s_{i}=0,1,2$ using Eq. (68) we get

[TABLE]

[TABLE]

[TABLE]

Furthermore, it turns out that $f_{a}(1)=f_{a}(-1)$, $f_{a}(2)=f_{a}(-2)$, etc., confirming the validity of the expression in Eq. (69). Consequently, the transformed partition function of Eq. (LABEL:eq:D0) becomes

[TABLE]

The procedure applies to any number of lattice sites.

We apply a similar procedure for indistinguishable particles. For the case $L=3$, the partition function obtained from Eq. (37) is

[TABLE]

then transforming it into the summations over $s_{i}$ we get

[TABLE]

where it now remains to obtain

[TABLE]

The result turns out to be

[TABLE]

where

[TABLE]

is the modified Bessel function of the first kind. To validate this expression we proceed as before. The results for $s_{i}=0,1,2$ are

[TABLE]

[TABLE]

[TABLE]

then $f_{b}(1)=f_{b}(-1)$, $f_{b}(2)=f_{b}(-2)$, etc., confirming the result in Eq. (77). The transformed partition function becomes

[TABLE]

Again, the procedure is valid for any number of lattice sites.

Acknowledgements.

D.F. would like to acknowledge financial support of the Federico Santa Maria Technical University via the “Programa 3: Apoyo a la Instalación en Investigación”.