Critical behavior of hard squares in strong confinement

TL;DR

This paper analytically studies the critical behavior of a confined quasi-one-dimensional system of hard squares, revealing a phase transition with Ising universality class at a specific channel width and infinite pressure.

Contribution

It introduces a restricted orientation and position approximation to analytically describe the critical transition in confined hard squares, extending previous transfer operator methods.

Findings

Identifies a critical point at H ≈ 1.8284 and infinite pressure.

Shows the transition belongs to the 1D Ising universality class.

Predicts similar behavior for unrestricted orientations and positions.

Abstract

We examine the phase behavior of a quasi-one-dimensional system of hard squares with side-length $\sigma$, where the particles are confined between two parallel walls and only nearest neighbor interactions occur. As in our previous work (PRE, 94, 050603 (2016)), the transfer operator method is used, but here we impose a restricted orientation and position approximation to yield an analytic description of the physical properties. This allows us to study the parallel fluid-like to zigzag solid-like structural transition, where the compressibility and heat capacity peaks sharpen and get higher as $H \rightarrow H_c=2\sqrt{2}-1\approx 1.8284$ and $p \rightarrow p_c= \infty$. Here $H$ is the width of the channel measured in $\sigma$ units and $p$ is the pressure. We have found that this structural change becomes critical at the $(p_c,H_c)$ point. The obtained critical exponents belong to the universality class of the one-dimensional Ising model. We believe this behavior holds for the unrestricted orientational and positional case.

Equations

$$\int d{y}^{\prime }d{\phi }^{\prime }K(y,\phi ;y^{\prime },{\phi }^{\prime })\psi (y^{\prime },{\phi }^{\prime })=\lambda \psi (y,\phi ).$$

$$K(y,\phi ;y^{\prime },{\phi }^{\prime })=\frac{e^{-p{\sigma }_x(y,\phi ;y^{\prime },{\phi }^{\prime })}}{p},$$

$$|0⟩$$

$$|1⟩$$

$$|2⟩$$

$$|3⟩$$

$$K=\frac{1}{p}\left(\begin{array}{cccc}{e}^{-p}& e^{-p\frac{1+\sqrt{2}}{2}}& e^{-p}& e^{-p(\frac{3}{2}+\sqrt{2}-H)}\\ e^{-p\frac{1+\sqrt{2}}{2}}& e^{-p\sqrt{2}}& e^{-p(\frac{3}{2}+\sqrt{2}-H)}& e^{-p(2\sqrt{2}-H)}\\ e^{-p}& e^{-p(\frac{3}{2}+\sqrt{2}-H)}& e^{-p}& e^{-p\frac{1+\sqrt{2}}{2}}\\ e^{-p(\frac{3}{2}+\sqrt{2}-H)}& e^{-p(2\sqrt{2}-H)}& e^{-p\frac{1+\sqrt{2}}{2}}& e^{-p\sqrt{2}}\end{array}\right).$$

$$U_y=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right),$$

$$K_{\pm }=\left(\begin{array}{cc}{\Lambda }_1^{\pm }& V^{\pm }\\ V^{\pm }& {\Lambda }_2^{\pm }\end{array}\right):=$$

$$\frac{1}{p}\left(\begin{array}{cc}{e}^{-p}\pm e^{-p}& e^{-p\frac{1+\sqrt{2}}{2}}\pm e^{-p(\frac{3}{2}+\sqrt{2}-H)}\\ e^{-p\frac{1+\sqrt{2}}{2}}\pm e^{-p(\frac{3}{2}+\sqrt{2}-H)}& e^{-p\sqrt{2}}\pm e^{-p(2\sqrt{2}-H)}\end{array}\right).$$

$$U_o^{(\pm )}=\left(\begin{array}{cc}{a}_1^{\pm }& a_2^{\pm }\\ a_2^{\pm }& -a_1^{\pm }\end{array}\right),$$

$$a_1^{\pm }=\pm \sqrt{\frac{1-S^{\pm }}{2}}{a}_2^{\pm }=\sqrt{\frac{1+S^{\pm }}{2}}$$

$$S^{\pm }=\frac{{\Lambda }_2^{\pm }-{\Lambda }_1^{\pm }}{\sqrt{({\Lambda }_2^{\pm }-{\Lambda }_1^{\pm }{)}^2+(2V^{\pm }{)}^2}}.$$

$$U$$

$$=\sum _{s\in \{+,-\}}{U}_y|s{⟩}_y⟨s|\otimes U_o^{(s)},$$

$$|{\psi }_0⟩$$

$$|{\psi }_1⟩$$

$$|{\psi }_2⟩$$

$$|{\psi }_3⟩$$

$${\lambda }_0$$

$${\lambda }_1$$

$${\lambda }_2$$

$${\lambda }_3$$

$$g:=\frac{\beta G}{N}=-\log ({\lambda }_0).$$

$$⟨{\mathcal{A}}_n⟩=⟨{\psi }_0|A|{\psi }_0⟩,$$

$$G_A(n)$$

$$=\sum _{k\ge 1}{\left(\frac{{\lambda }_k}{{\lambda }_0}\right)}^n⟨{\psi }_0|A|{\psi }_k⟩⟨{\psi }_k|A|{\psi }_0⟩.$$

$$⟨i|O|j⟩=\cos (4{\phi }_i){\delta }_{i,j}=\left(\begin{array}{cccc}-1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right).$$

$$S_o:=⟨{\mathcal{O}}_n⟩=(a_2^{+}{)}^2-(a_1^{+}{)}^2=S^{+},$$

$$G_o(n)=(2a_1^{+}{a}_2^{+}{)}^2{\left(\frac{{\lambda }_1}{{\lambda }_0}\right)}^n,$$

$$⟨i|y|j⟩$$

$$G_y(n)$$

$${\left(a_1^{+}{a}_2^{-}\frac{H-1}{2}-a_2^{+}{a}_1^{-}\frac{H-\sqrt{2}}{2}\right)}^2{\left(\frac{{\lambda }_3}{{\lambda }_0}\right)}^n.$$

$${\eta }^{-1}=H\frac{\partial g}{\partial p}.$$

$${\kappa }_{{}_T}:=-\frac{1}{L}\frac{\partial L}{\partial p_x}⟹\frac{k_{B}T}{{\sigma }^2}{\kappa }_{{}_T}=-H^2\eta \frac{{\partial }^{2}g}{\partial p^2},$$

$$G_o(n)=A_o{e}^{-n/{\xi }_o}$$

$$A_o=(2a_1^{+}{a}_2^{+}{)}^2$$

Full text

I Introduction

The structural and dynamical properties of confined fluids can be very different from that of bulk fluids. Upon dimensional restriction first order phase transitions transform into continuous ones Koga et al. (2001); Gordillo et al. (2006); Duran-Olivencia and Gordillo (2009); Han et al. (2010); Mochizuki and Koga (2015), jamming and glass formation occur Ashwin and Bowles (2009); Yamchi et al. (2012); Ashwin et al. (2013); Godfrey and Moore (2014, 2015); Robinson et al. (2016) and significant changes in some transport properties arise Wei et al. (2000); Kollmann (2003); Lin et al. (2005). In addition, other phenomena such as wetting, surface ordering, and layering transitions may occur in the presence of confinement Dijkstra et al. (2001); de las Heras et al. (2005, 2006); Geigenfeind et al. (2015). A fundamental issue in confined systems is to understand how the dimensional reduction of spatial variables changes the nature of phase transitions. It is a well-known result of the van Hove’s theorem van Hove (1950); Cuesta and Sanchez (2004) that particles interacting with short range potentials do not exhibit a genuine phase transition in one dimension. Therefore, it is interesting to study how the first order phase transitions become continuous in slit-like pores and cylindrical tubes. Along this line several molecular fluids and colloidal systems have been examined by changing the pore-width in slit-like geometries Radhakrishnan et al. (2002); Qi et al. (2014) and the pore-diameter in nanotubes Mochizuki and Koga (2015); Takaiwa et al. (2008). Common results of these studies are that the first order phase transitions observed in bulk may become continuous if one of the dimensions of the pore is reduced to the order of the size of the particle. This happens with the solid-liquid transition of water in nanotubes Mochizuki and Koga (2015) and with the isotropic-nematic transition of rod-like colloids in slit-like pores Cosentino Lagomarsino et al. (2003). Confinement may yield new types of phases such as the hexatic Radhakrishnan et al. (2002); Engel et al. (2013) and the biaxial-nematic ones Dijkstra et al. (2001). Moreover, unusual phase transitions can occur such as liquid-liquid Bianco and Franzese (2014); Wang et al. (2015) and smectic-smectic de las Heras et al. (2005, 2006). The criticality of confined water has been also studied in slit-like pores, where the liquid-liquid phase transition of water monolayers terminates at a critical temperature. The observed criticality belongs to the universality class of the two-dimensional (2d) Ising-model Bernard and Krauth (2011).

Along this line, the motivation of our work is to study the structural and critical behavior of a 2d system of confined hard squares, where the squares are allowed to rotate and move freely in a narrow channel of width $H$. In our previous work Gurin et al. (2016), we obtained that a fluid phase with parallel alignment to the wall abruptly transforms into a solid-like zigzag structure upon increasing pressure. We have found that our system shows extremely similar features than others presenting first order phase transitions. These results were based on numerical solutions of transfer operator equations and simulations, being both in very good accordance with each other. It is important to note that with our previous methods the sharpness of the structural transition was impossible to study for $H\rightarrow H_{c}$, $H_{c}$ being the pore width where the close packing densities of the parallel fluid and the zigzag structures are the same.

In this paper we further study the above described system in a pure analytical way by discretizing the rotational and transversal positions. We are employing dimensionless lengths; they are understood in $\sigma$ units where $\sigma$ is the side length of the square, e.g. the positional coordinates of the centers of squares are $x,y=\mathrm{distance}/\sigma$, $H=\mathrm{(width\ of\ the\ channel)}/\sigma$ etc. We set $H<2$ to satisfy the first neighbor interaction condition. This means that we have a quasi-one-dimensional system of classical particles which can be handled by the transfer operator method. It is well known that the equilibrium statistical physics of our system can be given by solving the following eigenvalue equation

[TABLE]

The kernel $K$ of the above integral operator is given by

[TABLE]

where $\sigma_{x}(y,\varphi;y^{\prime},\varphi^{\prime})$ is the $x$ projection of the contact distance of the nearest neighbor particles with orientations $\varphi$, $\varphi^{\prime}$ and positions $y$, $y^{\prime}$; furthermore $\beta=(k_{B}T)^{-1}$ is the inverse temperature and $p$ is the dimensionless longitudinal pressure (or force), $p=\beta p_{x}H\sigma^{2}$, where $p_{x}$ is the real longitudinal two-dimensional pressure in energy per square length units.

Note that this operator is clearly compact since it fulfills $\int dyd\varphi\int dy^{\prime}d\varphi^{\prime}\,|K(y,\varphi;y^{\prime},\varphi^{\prime})|^{2}<\infty$, i.e. the kernel type is Hilbert–Schmidt. Moreover, $K$ (eq. (2)) is everywhere positive and therefore the operator is irreducible. Finally, the integral operator is a positive operator in the sense that the image of all non-negative functions is non negative. For this case, based on the Perron–Frobenius–Jentzsch theorem van Hove (1950); Cuesta and Sanchez (2004), it can be proved that the dominant eigenvalue ($\lambda_{0}$) of the operator is unique, there is a gap between $\lambda_{0}$ and the remaining part of the spectrum, and $\lambda_{0}$ is an analytic function of $p$ (and $H$) Cuesta and Sanchez (2004). Thus, phase transitions (both, first order and continuous), as traditionally defined in the framework of statistical physics, are definitely out of the question.

The aim of the present paper is twofold. On the one hand we want to shed some light on the reason why the numerical and simulation (or even experimental) studies seem to predict a discontinuous behavior in spite of the above cited analytic theory. On the other hand, there still exists the demand to know whether the system shows real divergences at least at infinite pressure. To this end, we need to define an analytically solvable model, i.e. a simpler one, so that the limiting cases can be studied in an exact manner. Therefore, here we study a model in which both the orientation and the $y$ position of the particles are restricted to discretized values and only the $x$ position is continuous. To capture the main features of the freely rotating case, we have found that the minimal model must have at least four states for all particles, see Fig. 1. Fortunately, this model shows qualitatively similar structures to the continuum model. Even though the discretization of the orientational and positional degrees of freedom rises some fundamental questions Gurin and Varga (2011), we believe that the main features of the continuous system are captured. We will go back to this point at the end of the paper.

In the next section we present our notation and the analytic solution of the discretized model based on the transfer operator method. In the third section we examine the behavior of the orientational order parameter, correlation lengths and response functions, such as the isobaric heat capacity and the isothermal compressibility. Based on these results, we show that the system has a 1d Ising-like critical point at infinite pressure and a special channel width $H=H_{c}$. In the light of this result we revisit the problem of the structural transition between parallel fluid and a zigzag solid-like structures and we conclude that in the vicinity of the critical point this transition is practically indistinguishable from a genuine first order phase transition. Unfortunately, this $H\rightarrow H_{c}$ study is numerically impossible for the freely rotating case Gurin et al. (2016). Finally, in the last section we summarize our results and discuss the possible relation with glassy and jammed behavior, mentioned frequently in the literature in connection with other quasi-one-dimensional hard particle models.

II Analytic solution of a simple model

In our model the orientation of a particle can have two different values: $i)$ the side of the square is parallel to the wall, this orientational state of a particle is denoted by $|{\square}\rangle_{o}$ or $ii)$ the diagonal of the square is parallel to the wall, which is denoted as $|{\Diamond}\rangle_{o}$. (The index $o$ refers to the “orientation” and the “ket” and ”bra” symbols follow the notation given in Ref. Yeomans (1992).) Furthermore, we assume that the particles are always in contact with one wall, see Fig. 1. Thus, if the channel width is $H$ then the $y$ coordinate of a particle in the orientational state $|{\square}\rangle_{o}$ can be $\pm(H-1)/2$, and the $y$ coordinate of a particle in the orientational state $|{\Diamond}\rangle_{o}$ can be $\pm(H-\sqrt{2})/2$. These $y$ positional states—regardless of the difference due to the orientation—are denoted by $|{+}\rangle_{y}$ and $|{-}\rangle_{y}$.

These four possible combination of orientations and $y$ positions form a basis in the four dimensional $\mathcal{H}_{y}\otimes\mathcal{H}_{o}$ space. These basis states can be written as

[TABLE]

and the transfer operator acts on the space spanned by the above orthonormal basis. The matrix of the transfer operator (the transfer matrix) in this basis can be written as 111Here we have supposed that $\sqrt{2}/2+1<H<2$. The upper bound is necessary, but the calculations can be extended without any problem in the $\sqrt{2}<H<\sqrt{2}/2+1$ domain when the only difference is that $K_{0,3}=K_{1,2}=K_{2,1}=K_{3,0}=e^{-p\frac{1+\sqrt{2}}{2}}/p$ which has no any interesting consequence. We will see that the $H<2\sqrt{2}-1$ case is unimportant.

[TABLE]

Certainly, the operator represented by this matrix is positive, irreducible and compact, therefore the theorem presented in Cuesta and Sanchez (2004) is valid for this case, i.e. a phase transition is impossible. The model under study has a simple property, namely that there is no “entanglement” between the orientational and the $y$ degrees of freedom in the sense which is clarified below. Therefore, matrix (4) can be diagonalized in an easy way. For this purpose we define a unitary operator $U_{y}\otimes 1_{o}$, where $U_{y}$ acts only on the $y$ degrees of freedom. Its matrix in the $|{+}\rangle_{y}$, $|{-}\rangle_{y}$ basis is

[TABLE]

and $1_{o}$ is the unit operator acting only on the orientational degrees of freedom. The unitary operator $U_{y}\otimes 1_{o}$ uncouples the $|{\psi_{+}}\rangle_{y}:=U_{y}|{+}\rangle_{y}=\frac{|{+}\rangle_{y}+|{-}\rangle_{y}}{\sqrt{2}}$ and $|{\psi_{-}}\rangle_{y}:=U_{y}|{-}\rangle_{y}=\frac{|{+}\rangle_{y}-|{-}\rangle_{y}}{\sqrt{2}}$ degrees of freedom. This means that $K$ has two pieces of two dimensional invariant subspaces: $\mathcal{H}_{y}^{+}\otimes\mathcal{H}_{o}$ and $\mathcal{H}_{y}^{-}\otimes\mathcal{H}_{o}$ (where $\mathcal{H}_{y}^{\pm}$ denotes the one dimensional subspace in $\mathcal{H}_{y}$ generated by the vector $|{\psi_{\pm}}\rangle_{y}$). In other words, the matrix of $K$ after the $U_{y}\otimes 1_{o}$ transformation is a block diagonal matrix, which has the following $2\times 2$ blocks in its diagonal:

[TABLE]

The above $K_{+}$ and $K_{-}$ symmetric matrices can be diagonalized by unitary operators $U_{o}^{(+)}$ and $U_{o}^{(-)}$, which act only on the orientational degrees of freedom and their matrices (in the $\{|{\square}\rangle,|{\Diamond}\rangle\}$ basis) can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

Now we can summarize our results detailed above. $K$ can be diagonalized by the unitary operator

[TABLE]

that is, the eigenvectors of $K$ are $|{\psi_{i}}\rangle=U|{i}\rangle$ and the corresponding eigenvalues are $\lambda_{i}=\langle{i}|U^{\dagger}KU|{i}\rangle$. In detail, the eigenvectors are

[TABLE]

and taking into account that $\Lambda_{1}^{-}=0$, the corresponding eigenvalues are

[TABLE]

It can be seen that for any value of $p$ and $H$ we have $\lambda_{0}>\lambda_{1}>\lambda_{2}>0>\lambda_{3}$, however $\lambda_{0}>|\lambda_{3}|>\lambda_{1}>\lambda_{2}$.

As mentioned above, the orientational and $y$ positional degrees of freedom are not “entangled”. We mean that all the eigenvectors have a form $|{\psi}\rangle_{y}\otimes|{\psi^{\prime}}\rangle_{o}$. That is the reason why we need to solve only two quadratic equations instead of a quartic one. This simple feature of the transfer operator does not hold for the freely rotating case. From the general transfer operator theory follows (see e.g. Yeomans (1992)) that the Gibbs free energy is given by

[TABLE]

Having a one-particle physical quantity $\mathcal{A}_{n}$, i.e. for a given microscopic state of the system $\mathcal{A}_{n}$ has four (in general) different values, $a(i)$, depending on the state only of the $n$-th particle which is labeled by $i$, we can define an operator of which matrix is diagonal in the basis given by Eqs. (3): $\langle{i}|A|{j}\rangle=a(i)\delta_{i,j}$. Now the expectation value of $\mathcal{A}_{n}$ (which is, in our case, certainly independent of the label of the particle, $n$, because we have no positionally dependent external fields) can be written as

[TABLE]

and the correlation function between $n$-th neighboring particles is given by

[TABLE]

The orientation and position of a particle give examples for one-particle properties. If we define $\varphi$ as the orientation of a particle in state $|{\square}\rangle$ such as $\varphi=\pi/4$ and in state $|{\Diamond}\rangle$ such as $\varphi=0$, then the operator $O$, related to the one particle quantity $\mathcal{O}_{n}=\cos(4\varphi_{n})$, can be represented by the matrix

[TABLE]

Using Eqs. (11a), (8) and (9) we find that the orientational order parameter, $S_{o}$, which is the expectation value of $\mathcal{O}_{n}$, can be expressed as

[TABLE]

and the correlation function as

[TABLE]

because from Eq. (16) it can be seen that $\langle{\psi_{2}}|O|{\psi_{0}}\rangle=\langle{\psi_{3}}|O|{\psi_{0}}\rangle=0$. In a similar way we can define the operator of the $y$ position, denoted also by $y$, which can be represented by the matrix

[TABLE]

From the above definition we get $\langle y\rangle=0$ and

[TABLE]

Note that the second term alternates sign with $n$ since $\lambda_{3}$ is the only negative eigenvalue.

III Analysis of the results

Figure 2 shows the pressure, $p$, as a function of the packing fraction, $\eta:=N/(HL)$. Here $L$ is the length of the channel (normalized by $\sigma$) along the $x$ axis. The equation of state can be calculated from the result (12a) using Eq. (13), yielding

[TABLE]

Although $p(\eta)$ is obtained from an analytic formula, it shows a plateau which increases its resemblance with a first order discontinuity as $H\rightarrow H_{c}=2\sqrt{2}-1\approx 1.8284$. Note that for the case of $H=1.83$ (see the inset of Fig. 2), it already turns almost impossible to plot the curve as a smooth function, and this feature worsens for $H\rightarrow H_{c}$. Therefore, although continuous, the system behavior cannot be practically distinguished from a genuine first order transition neither by simulations nor—in case this system could have an experimental realization—by real experiments.

The situation is similar to the case of a finite but large ($N\approx 10^{23}$) system. We know from the statistical physics that in a finite system all derivatives of the free energy are continuous. The singularities, in the mathematical sense, appear only when considering the thermodynamic limit. However, although real systems consist of a finite but very large number of particles, experiments clearly show all significant features of phase transitions. It is generally accepted that the freezing of one liter of water is a genuine phase transition even if the system is finite and so there is no mathematical singularity.

Before we discuss the reason for our system behavior we show how other thermodynamic properties also depict quasi-singularities. The isothermal (and longitudinal) compressibility, given by

[TABLE]

is shown in Fig. 3.

As can be seen, $\kappa_{{}_{T}}(p)$ peaks at the structural transition, while turning extremely sharp as $H\rightarrow H_{c}$.

Now we examine the order parameters and their correlation functions. The orientational order parameter (see Eq. 17) can be seen in Fig. 4. At very low pressure the system behaves as an ideal gas and both $\square$ and $\Diamond$ orientations turn equally probable. This is simply because we are considering the same number of $y$-positions for both $\square$ and $\Diamond$ orientations. Thus, this behavior differs from that of the freely $y$-positioning system, since in this last case parallel configurations are favored by entropy, i.e. there are more configurations for the $\square$-orientation than for the $\Diamond$-orientation for an isolated particle. With increasing pressure the orientation $\square$ is preferred because the $\square$-$\Diamond$ pair at contact has a large $x$-projection and the $\Diamond$-$\Diamond$ pair is favored only for large $\Diamond$ clusters Gurin et al. (2016). Then, at a given $\tilde{p}(H)$ pressure a structural change happens in the system and the orientation $\Diamond$ becomes more favored. Clearly, the $\square$ cannot hold as the packing fraction surpasses $\eta_{cp}^{\square}=1/H$, but the structural transition takes place slightly below this value. The orientational correlation function, according to Eq. (18), can be written as

[TABLE]

with

[TABLE]

and

[TABLE]

Note that Eq. (23) is not only asymptotically but exactly valid for any neighboring distance $n$. The amplitude of the orientational correlations, $A_{o}$, is depicted in Fig. 4 and the orientational correlation length, $\xi_{o}$, is shown in Fig. 4.

Te orientational correlation length increases monotonically with increasing pressure at small densities, but the amplitude of the correlation function decreases and the correlations almost vanish even for $p\lessapprox\tilde{p}$ if $H\gtrapprox H_{c}$. Then, the amplitude has a sharp peak when the correlation length is maximal, indicating the strong and large scale fluctuations in the system at this special point. Finally, both the amplitude and the correlation length go down as $p$ approaches infinity.

All the above described peculiar behavior happens at a $\tilde{p}(H)$ pressure which can be determined from the $\Lambda^{+}_{1}=\Lambda^{+}_{2}$ condition. When this equality is fulfilled at high pressure, the off-diagonal elements of $K_{+}$ can be neglected (because $V^{+}\ll\Lambda^{+}$ at high $p$). If we suppose, as an approximation, that $V^{+}=0$, then the two different eigenvalues of $K_{+}$, namely $\lambda_{0}$ and $\lambda_{1}$, cross each other as $p$ is increased, and the system yields a first order transition at $\tilde{p}$ from a $\square$ to a $\Diamond$ oriented structure. Under this approximation, the Perron–Frobenius arguments cannot longer be applied, because the transfer matrix has zero entries. Thus, the phase transition appears only by forcing the off-diagonal terms in $K_{+}$ to zero. By removing this imposition, the level crossing of the two eigenvalues is avoided, i. e. the two eigenvalues only approach each other and the system does not yield a first order transition. However, the closer the eigenvalues approach each other, the more the behavior of the system gets reminiscent to that of a phase transition. The strength of this reminiscence increases with decreasing $H$ while $H>H_{c}$. This happens due to the fact that when $H$ approaches $H_{c}$, $\tilde{p}(H)$ increases and the level crossing approximation (the negligence of $V^{+}$) becomes more valid.

The $\Lambda^{+}_{1}=\Lambda^{+}_{2}$ equation can be written, in more detail, as

[TABLE]

In case of large pressure ($e^{-p(H-\sqrt{2})}\ll 1$) equation (26) simplify as

[TABLE]

As we have mentioned, for a given $H$ the solution of Eq. 26 determines the $\tilde{p}(H)$ value of the pressure, where the system behaves similarly as presenting a first order transition. Hereafter we call $\tilde{p}(H)$ as the level crossing pressure. Alternatively, if we keep the pressure fixed and deal with $H$ as a control parameter, from Eq. (26) we get the level crossing value of the channel width, $\tilde{H}(p)$. From here on we use the notations $\tilde{\Lambda}(p)=\Lambda^{+}_{1}(\tilde{H}(p),p)=\Lambda^{+}_{2}(\tilde{H}(p),p)$ and $\tilde{V}(p)=V^{+}(\tilde{H}(p),p)$. It is easy to show that $\lim_{H\rightarrow H_{c}}\tilde{p}(H)=\infty$ and $\lim_{p\rightarrow\infty}(\tilde{V}(p)/\tilde{\Lambda}(p))=0$, which explain why the level crossing approximation becomes more valid as $H$ approaches $H_{c}$.

It is important to note that both the orientational correlation length and the compressibility truly diverge as $(p,H)\rightarrow(\infty,H_{c})$ in such a special way that the condition $\Lambda^{+}_{1}=\Lambda^{+}_{2}$ is always fulfilled. When we go on this level crossing line, for large $p$ (from which follows that $\tilde{V}\ll\tilde{\Lambda}$), the orientational correlation length goes with $p$ as

[TABLE]

from where we observe that the orientational correlation length diverges exponentially. Further in the text we argue that this is really a critical divergence.

III.1 The scaling property of the orientational correlation function and the Gibbs free energy

When we see more carefully the correlation function given by Eq. (23), we can observe that on the level crossing line the amplitude of the correlation function equals one ($\Lambda^{+}_{1}=\Lambda^{+}_{2}\Rightarrow S^{+}=0\Rightarrow 2a_{1}^{+}a_{2}^{+}=A_{o}^{1/2}=1$). That is, the orientational correlation function,

[TABLE]

depends only on $n/\xi_{o}$, as it is usual near a critical point, where the system shows a scaling behavior. In general, the scaling form of the correlation function is written as (see e.g. Ref Chaikin and Lubensky (1997) or Cardy (1996))

[TABLE]

By comparing Eqs. (29) and (30) we must conclude that in our case $d-2+\eta=0$.

We can now see that in the vicinity of $(p_{c},H_{c})$ not only the orientational correlation function but other physical properties show a scaling behavior when expressed in terms of the correlation length. This is a common practice in one dimensional systems, rather than expressing the properties in terms of the reduced temperature or pressure. In the limit of $p\rightarrow p_{c}$ (it follows that $\tilde{V}\ll\tilde{\Lambda}$) and $H\approx\tilde{H}\rightarrow H_{c}$ (it follows that $\Lambda^{+}_{1}\approx\Lambda^{+}_{2}\approx\tilde{\Lambda}$ and that $|\Lambda^{+}_{1}-\Lambda^{+}_{2}|/2\ll\tilde{V}$) one can get

[TABLE]

and

[TABLE]

Using Eq. (28) this last expression can be identified with the general scaling form for the singular part of the free energy Chaikin and Lubensky (1997); Cardy (1996)

[TABLE]

where $g_{0}$ is the regular part of $g$, $h$ is some kind of external field, and $X$ is a scaling function. By comparing Eqs. (32) and (33) we find that

[TABLE]

where $g_{0}=-\log\tilde{\Lambda}$ and the external field $h$ measures the distance from the level crossing line defined by Eq. (26),

[TABLE]

The scaling property of the correlation function and the free energy (Eqs. (29) and (34)), prove that in the vicinity of the $(p_{c},\ H_{c})$ point our model shows critical behavior. From Eqs. (33, 34) we must conclude that our system really behaves as a one dimensional system, $d=1$; moreover $\lambda=1$. The two independent exponents can be chosen as $\eta=1$ and $\lambda=1$, and the other usual exponents can be calculated from the scaling laws. Alternatively, we can compute directly the scaling of the isothermal compressibility and the specific heat. We observe from Eq. (34) that

[TABLE]

because $X(a)=1+2a^{2}$, and therefore $X^{\prime}(0)=0$. But $X^{\prime\prime}(0)=4$, and so the second derivative of $g$ is

[TABLE]

Now, taking into account Eq. (22) we can see that $\kappa_{{}_{T}}\sim\xi$. A comparison with the usual definition of the exponent $\bar{\gamma}=\gamma/\nu$, from which $\kappa_{{}_{T}}\sim\xi^{\bar{\gamma}}$, leads to the conclusion that $\bar{\gamma}=1$. Similarly, the isobar specific heat (isobar in the sense that the longitudinal pressure $p$ is constant, but also $H=const.$) can be written as

[TABLE]

and therefore, from Eqs. (36) and  (37) we obtain $c_{p}\sim\xi$, implying that the exponent $\bar{\alpha}=\alpha/\nu$, such that $c_{p}\sim\xi^{\bar{\alpha}}$, equals one.

By comparing Eqs. (26) and (35) one can conclude that $h=0$ means that the pressure is $\tilde{p}(H)$, i.e. the system is at the level crossing line where the order parameter is zero. Therefore, $\lim_{h\to 0}\langle S_{o}\rangle=0$. In other words, there is not a spontaneous formation of an ordered phase in our model (which is not surprising in quasi-1d). Consequently, the value of $\bar{\beta}=\beta/\nu$, the critical exponent related to the order parameter, cannot be computed directly. However, from the scaling laws we can determine its value. We have $\lambda=1$, $\bar{\alpha}=\bar{\gamma}=1$ and $\eta=1$. From the scaling law $\lambda=\bar{\beta}+\bar{\gamma}$ we conclude that $\bar{\beta}=0$. Alternatively, from Eq. (34) we have seen that the dimension is $d=1$, and from the scaling law $\bar{\beta}=\frac{1}{2}(d-2+\eta)$ we consistently obtain $\bar{\beta}=0$. The effect of $\bar{\beta}=0$ is observed in Fig. 4, that is, as $\tilde{p}\rightarrow p_{c}$ on the level crossing line, both $\langle S_{o}(H)\rangle|_{p=const.}$ and $\langle S_{o}(p)\rangle|_{H=const.}$ go to the step function. This is a very important difference compared with the usual critical points in three dimensions.

From the obtained critical exponents we can see that the $(p_{c},H_{c})$ critical point belongs to the universality class of the 1d Ising model. The important consequence of this result, as we discussed in the previous paragraph, is that the critical exponent $\bar{\beta}=0$ (contrasting with the usual $\beta=1/3$ value of 3d systems), implies that the order parameter is discontinuous at the $(p_{c},H_{c})$ point. Despite the fact that this discontinuity disappears at any finite pressure, it has a significant impact on the system behavior near the critical point. Namely, the order parameter and the density behave almost like the step function, and their derivatives, as the compressibility, have high peaks. Simulation results, as close as possible to the $(p_{c},H_{c})$ point and for the unrestricted $y$ and $\varphi$ system, show that the system behavior is indistinguishable from that of a first order transition Gurin et al. (2016).

The discontinuity can disappear at any finite $p$ because the thermodynamic quantities are singular at the $(p_{c},H_{c})$ fixed point, contrary to the case of a zero temperature discontinuity fixed point of the usual 3d systems, like the ferromagnets. In this last case the thermodynamic quantities are not singular at a discontinuity fixed point which separates the different phases, and perturbation theory should have a finite radius of convergence, therefore the reason of the discontinuity, the coexistence of the different phases must persist for some distance into the phase diagram. This coexistence line is terminated at a different, critical fixed point Cardy (1996). In 1d, these two fixed points, the critical and the discontinuity ones are merged in a unique fixed point. This is the origin of the peculiar behavior of our 1d system, which is reminiscent of a first order phase transition though every property can be expressed as an analytic function.

Let us now discuss the physical origin of the “external field” given by Eq. (35). By extending our model with an external field which favors the $\Diamond$ orientation against with the $\square$ one, including in the Hamiltonian a $-h^{\prime}\sum_{i}\cos(4\varphi_{i})$ term, this extra field $h^{\prime}$ simply appears as an additive term in Eq. (35). This points out that $h$ has exactly the same effect as an extra external field $h^{\prime}$ which is directly coupled to particle orientation. But in our system the special combination of two different external fields—the external pressure $p$ and the wall of the channel—results in an effective field which favors one orientation or the other. Moreover, the strength and the “direction” of this effective field depend on the values of $p$ and $H$, given that the longitudinal pressure favors the closely packed $\Diamond$ orientation (having vanishing $y$-fluctuations) and the wall the $\square$ orientation. In other words, the direction of this external field results from an entropic competition between the $x$ and $y$ fluctuations.

We would like to emphasize the following interesting feature of the $(p_{c},\ H_{c})$ point: the compressibility diverges as $p\rightarrow\infty$. This is quite special in a system consisting of only rigid particles. Usually $p\rightarrow\infty$ implies that the system approaches the close packing structure while the compressibility goes to zero. But in our system, at $(p_{c},\ H_{c})$ the close packing structure is degenerated, because the $\square$ and the $\Diamond$ orientations have the same close packing density. This is the key feature of this point: the system cannot decide between these two competitive structures. In this sense, it is very similar to a spontaneous symmetry breaking. Note that other systems may show this peculiar point. For instance, hard anisotropic particles get spatially and orientationaly ordered at close packing, whereas spheres (or disks) get only spatially ordered. In the limit of small anisotropy, both, the orientationaly ordered and disordered structures also have the same packing fraction, and the structure gets again degenerated  Bautista-Carbajal et al. (2013); Bautista-Carbajal and Odriozola (2014).

III.2 The $y$ positional correlations

Up to this point we have focused only on the orientational correlations in our model, which can be totally described by the $2\times 2$ transfer matrix $K_{+}$ given by Eq. (6). Thus, the analogy with the 1d Ising model comes as no surprise because this last model can be described by a $2\times 2$ transfer matrix, too. However, in addition to the continuous longitudinal translational degrees of freedom, every particle has 4 discrete possible states, instead of the 2 states of the Ising model.

To evaluate the $y$ positional correlation function we have to take into account the 2nd and 3rd eigenvalues and eigenvectors of the transfer operator, and so the above mentioned $2\times 2$ matrix, $K_{+}$, is not enough anymore. The expectation value $\langle y\rangle$ is zero, as we have mentioned at the end of Sec. II. Its correlation function, according to Eq. (20), can be written as

[TABLE]

and the corresponding results are depicted in Fig. 5. This last equation is again exactly valid for any neighboring distance $n$. The amplitude of the $y$ correlations, $A_{y}=(a_{1}^{+}a_{2}^{-}\frac{H-1}{2}-a_{2}^{+}a_{1}^{-}\frac{H-\sqrt{2}}{2})^{2}$ and $A_{y}^{\prime}=(a_{1}^{+}a_{1}^{-}\frac{H-1}{2}+a_{2}^{+}a_{2}^{-}\frac{H-\sqrt{2}}{2})^{2}$ are shown in Fig. 5. The correlation length, $\xi_{y}=(-\log(|\lambda_{3}|/\lambda_{0}))^{-1}$ and also $\xi_{y}^{\prime}=(-\log(\lambda_{2}/\lambda_{0}))^{-1}$ are shown in Fig. 5.

From Fig. 5 it can be seen that the positional correlations show mainly a non-alternating behavior when $p<\tilde{p}(H)$, but suddenly become alternating for $p>\tilde{p}(H)$. The picture is as follows: when the orientation of the particles is $\square$, then the correlation is weak, non-alternating, and short range (nevertheless, the alternating part has relatively long correlation length and small amplitude, which is due to the presence of few $\Diamond$ particles.) When the orientation of the particles is mainly $\Diamond$, then the correlation is strong, alternating, and long range.

It is interesting that the correlation length $\xi_{y}$ suddenly increases at $\tilde{p}$ but has no peak; it is a monotonic function and diverges exponentially with $\tilde{p}$. This kind of divergence has been also observed in other systems of hard body particles such as rectangular Kantor and Kardar (2009) or V-shaped Khandkar et al. (2017) particles confined to a line. The divergence of $\xi_{y}$ can be explained easily as follows. If—due to a fluctuation—a particle changes its position in the $\Diamond$ oriented zig-zag phase, then the neighboring particles are forced to follow that change to avoid forming domain walls in the zig-zag structure. This means that the correlation among the $y$ positions is strong. This effect becomes stronger with increasing pressure because the cost of domain walls increases, therefore the correlation length increases, too.

It is interesting to see that at the level of the transfer operator the reason of this sudden change of $\xi_{y}$ is due to the fact that $\lambda_{0}>|\lambda_{3}|>\lambda_{1}$, where $\lambda_{0}$ and $\lambda_{1}$ produce an avoided level crossing, thus $|\lambda_{3}|$ gets stuck between $\lambda_{0}$ and $\lambda_{1}$. The orientational and $y$ positional degrees of freedom are not “entangled”, but the orientational level crossing has some impact on the positional behavior.

Finally, in our opinion, this divergence of $\xi_{y}$ is not so interesting as the divergence of $\xi_{o}$, because the latter causes an interesting behavior at finite pressure (the peculiar behavior of the equations of state, see Fig. 2) and also at infinite pressure (diverging compressibility in spite of dealing with a hard system), but the former has not such consequences. The reason for this is that the orientational order is coupled to both physical external fields, the wall and the longitudinal pressure—see Eq. 35—, while the alternating $y$ positional order is not. The zig-zag positional order would be coupled to an alternating external field, i.e. the Hamiltonian should contain a $h^{\prime\prime}\sum_{i}(-1)^{i}y_{i}$ term, which is physically unusual, $h^{\prime\prime}$ cannot be tuned by $p$ or $H$, therefore the related quantities (susceptibility etc.) are uninteresting.

IV Summary and further discussions

We have shown that near the $H_{c}=2\sqrt{2}-1$, $p_{c}=\infty$ point the Gibbs free energy and the orientational correlation function show a scaling behavior, so this is a critical point. We have calculated all the critical exponents and we have found that our model is in the same universality class as the 1d Ising model. One can argue that there is no real critical point at any finite pressure, however, there is a real (experimentally observable) critical behavior in the vicinity of the critical point, which is indeed located in the physically meaningless (experimentally unreachable) parameter regime ($p_{c}=\infty$). At this fixed point the critical exponent $\bar{\beta}=0$, implying that the order parameter has a discontinuity. This means that this fixed point unifies the feature of a usual critical fixed point and a discontinuity fixed point of three dimensional systems. Therefore, in the vicinity of this point the system behaves very similar as showing a first order transition, however, at the same time, the peaks of the compressibility and the specific heat are typically like that of a critical system. The free energy is singular at $(p_{c},H_{c})$, therefore the discontinuity disappears at any finite $p$, turning all thermodynamic functions analytic.

In spite of the similarities between our system and the 1d Ising model, we want to emphasize some differences in the underlying physics. We have mentioned that in our model every particle has continuum translational $x$ degrees of freedom and four possible discrete states while the Ising spins are localized and have only two different discrete states, but more important, in the Ising model the thermal fluctuations can change the direction of a single spin alone at any temperature. Certainly, when the temperature is small, the probability of spin flipping is very small too, but possible, irrespectively to the states of its neighbors. In our model the ”flipping” of a square from the rotational state parallel to the wall ($\square$) into the other state when its diagonal is parallel to the wall ($\Diamond$) is impossible at high densities without the disturbance of its neighbors. At high densities fluctuations can ”flip” a square only together with many other neighbors, what is more, the positions of the squares have to be changed at the same time. Only collective motions can change the orientational state of the particles, which makes an important difference.

We would like to add that this kind of model, namely hard particles confined into a narrow, quasi-one-dimensional channel, is often regarded as a simple model to study the glassy or jamming phenomenon, see for example Refs. Yamchi et al. (2012); Ashwin et al. (2013); Godfrey and Moore (2015); Robinson et al. (2016). Most of these works focus on disks confined by somewhat wider channels than ours, but the observed glassy behavior is similar since only collective rearrangements of particles can increase the density, and with increasing pressure these collective rearrangements become less probable; the system is sticked into the so called locally jammed states. These works report interesting features such as an isobaric heat capacity maximum and a corresponding diverging pressure for some density below the maximal packing fraction. Then, they conclude that these features are related to some kind of fragile-strong fluid crossover or avoided phase transition, which is phenomenologically reminiscent of the bulk glass transition. Here we would like to provide an alternative point of view. We think that the almost singular behavior of the thermodynamic quantities below the close packing density can be a marker of a critical behavior related to a critical point at the close packing density. The system studied in our work gives an example for this possibility. The existence of a critical point at $(p_{c},H_{c})$ has consequences at finite $\tilde{p}$ pressures for $H>H_{c}$, namely: $c_{p}$ and $\kappa_{T}$ have large peaks, and the pressure goes up suddenly at a given (below the close packing) density. These features are very similar to those reported in connection with the fragile-strong fluid crossover or avoided phase transition. This behavior can be explained coherently by the existence of a fixed point at infinite pressure, which unifies the properties of a usual critical fixed point and a discontinuity fixed point.

It remains an open question what is the relationship between the present criticality and the possible glass/jam behavior. Here we would like to point out the differences. First of all, it is a long standing question whether the jamming transition has some sign in the equilibrium properties of the system or not. We emphasize that we studied only equilibrium properties. On the other hand, to understand the critical behavior in our model it is enough to take into account only two competitive structures. This contrasts with a glass, where the system have not only two, but many (usually very much) almost stable but actually metastable states. Moreover, in our system, the reason of the existence of metastable states is the presence of the confining walls, which strongly decrease the room for rearranging configurations. Finally, the existence of a merged critical and discontinuity fixed point is typical in one dimension. Therefore, the extrapolation of the conclusions deduced from such quasi-1d systems to bulk 3d systems is, from our point of view, strongly questionable.

Another issue is the analogy between the present and the continuous models. How relevant are the presented results for the freely rotating and moving case, when all the degrees of freedom are continuous? We have no exact answer, but we strongly believe that for the case where the $y$ degrees of freedom are continuous with discretized orientation, the long range orientational properties of our model can be effectively described by a $2\times 2$ matrix and, as a consequence, the 1d Ising like critical point is preserved. But even in this case the system cannot be completely resolved so easily. As we have mentioned, the orientational and $y$ positional degrees of freedom turn “entangled” and it becomes not so trivial to construct the effective $2\times 2$ matrix. On the other hand, treating the continuous rotational degrees of freedom is mathematically more subtle. Nevertheless, physical considerations suggest that the underlying reason for the “almost” singular behavior at pressure $\tilde{p}(H)$ is not the discrete nature of the ($y$ and orientational) degrees of freedom. What it is more, the extra degrees of freedom enhance the singular-like behavior even for not so strong confinement. In the discrete system studied in this paper there is nothing interesting when $H>1.9$, and the curves just start to show singular-like behavior below $H=1.85$, see e.g. Fig. 2 or 3. We know that the situation is more exciting in wider channels when $y$ is continuous Gurin et al. (2016). The reason is that the negative constant part of the external field $h$ in Eq. (35) comes from the $y$ positional fluctuations. When $y$ is continuous, fluctuations have more room as the phase space is larger. Therefore $h$ is more negative and favors more strongly the $\square$ orientation. The emerging picture is as follows: $y$ fluctuations favor the $\square$ orientation. Therefore, by increasing the number of $y$-degrees of freedom the “transition” pressure, $\tilde{p}$, goes up, approaching the critical point. Our numerical results show Gurin et al. (2016) that the freely rotating case shifts its “transition” pressure further more to high values, although in this case it is not trivial the reason why. In this unrestricted system we have no prove of the existence of a critical point, but we strongly believe on its presence.