Finite-size effects and thermodynamic limit in one-dimensional Janus fluids

New materials chemical technology allows for the synthesis of colloidal-size particles with patches exhibiting an interaction pattern different from that of the rest of the surface [1–3]. When the patch occupies a hemisphere, we are in the presence of so-called Janus particles [3–8].

One-dimensional fluids play an important role in statistical mechanics because they often offer integrable systems [9–34]. In a recent paper [35], two of us derived the exact equilibrium thermodynamic and structural properties of one-dimensional Janus fluids in the thermodynamic limit (TL). The system consisted in a binary mixture of two-face N_i = x_i N particles of species i = 1, 2, where x_i is the mole fraction of species i and N is the total number of particles. See figure 1 for a sketch of the system. In this type of systems (henceforth referred to as quenched), the number of particles (N₁ and N₂) with each face orientation is kept fixed but of course one needs to average over all possible microscopic configurations to obtain macroscopic quantities. Interestingly, the theoretical predictions for quenched systems agreed excellently well with Monte Carlo (MC) simulations for annealed systems (where at each MC attempt a particle is assigned the face orientation 1 or 2 with probabilities q₁ and q₂ = 1 − q₁, respectively) with N = 500.

Zoom In Zoom Out Reset image size

The investigation of [35] stimulates a few questions: (i) can the exact derivation of the Gibbs free energy in the TL (N → ∞) be extended to quenched and/or annealed finite-N systems? (ii) Does the quenched ↔ annealed equivalence break down at finite N? (iii) Can those theoretical predictions be validated by MC simulations? (iv) Is the dependence of the average mole fraction ⟨x₁⟩ on the probability q₁ robust with respect to N in annealed MC simulations for biased situations (${q}_{1}\ne \frac{1}{2}$)? The main aim of this paper is to address those questions. As will be seen, the answers are affirmative in all the cases.

The remainder of this paper is organized as follows. Section 2 presents the derivation of the configuration integral, and hence of the Gibbs free energy G, for a finite-size quenched binary mixture in the isothermal–isobaric ensemble. Those results are then used in section 3 to derive G for an annealed fluid. Since the exact results in sections 2 and 3 apply to any choice of the two nearest-neighbor interaction potentials ϕ₁₁ = ϕ₂₂ = ϕ₂₁ and ϕ₁₂ (see figure 1), the expressions are particularized in section 4 to the Kern–Frenkel model [36], where ϕ₁₁ and ϕ₁₂ are the hard-rod and square-well potentials, respectively. The theoretical results are validated and confirmed by MC simulations in section 5, where also the case of biased annealed systems is addressed. Finally, the main results of the work are summarized in section 6. The most technical parts of the paper are relegated to five appendices.

2.1. The system

Let us consider a one-dimensional binary fluid mixture made of N₁ particles of species 1 (right 'spin') and N₂ = N − N₁ particles of species 2 (left 'spin') on a line of length L (see figure 1). Henceforth, we will use Latin and Greek indices for species and particles, respectively. A particular spatial configuration will be denoted as x ≡ {x_α ; α = 1, 2, ..., N}. Analogously, a particular spin (or species) configuration will be denoted as s ≡ {s_α ; α = 1, 2, ..., N}, where s_α = 1, 2 represents the spin of particle α. Since we are considering a quenched mixture, the number of possible spin configurations are restricted by the constraint

$$\begin{equation}\sum\limits _{\alpha =1}^{N}\enspace {\delta }_{{s}_{\alpha },1}={N}_{1}.\end{equation} \tag{ 2.1 }$$

The total number of allowed spin configurations is $\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)$.

We assume that the rods are impenetrable and that their interaction is restricted to nearest neighbors. Given s and x , the total potential energy can be written as

$$\begin{equation}{{\Phi}}_{N}(\boldsymbol{s},\boldsymbol{x})=\sum\limits _{\alpha =1}^{N-1}\enspace {\phi }_{{s}_{\alpha },{s}_{\alpha +1}}({x}_{\alpha +1}-{x}_{\alpha })+\omega {\phi }_{{s}_{N},{s}_{1}}({x}_{1}+L-{x}_{N}),\end{equation} \tag{ 2.2 }$$

where, without loss of generality, we assume that particles 1, 2, ..., N are ordered from left to right. In equation (2.2), ω = 1 if periodic boundary conditions are applied and ω = 0 otherwise (open systems).

2.2. Isothermal–isobaric partition function

In the isothermal–isobaric ensemble, the partition function is [32, 37]

$$\begin{equation}{\mathcal{Z}}_{{N}_{1},{N}_{2}}(\beta ,\gamma )={\mathcal{Z}}_{{N}_{1},{N}_{2}}^{\text{id}}(\beta ,\gamma ){\mathcal{Q}}_{{N}_{1},{N}_{2}}(\beta ,\gamma ),\end{equation} \tag{ 2.3 }$$

where

$$\begin{equation}{\mathcal{Z}}_{{N}_{1},{N}_{2}}^{\text{id}}(\beta ,\gamma )=\frac{{C}_{{N}_{1},{N}_{2}}}{{L}_{\text{ref}}{\left[{{\Lambda}}_{1}(\beta )\right]}^{{N}_{1}}{\left[{{\Lambda}}_{2}(\beta )\right]}^{{N}_{2}}},\quad {C}_{{N}_{1},{N}_{2}}(\gamma )\equiv \left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right){\gamma }^{-(N+1)},\end{equation} \tag{ 2.4 }$$

is the ideal-gas partition function and

$$\begin{equation}{\mathcal{Q}}_{{N}_{1},{N}_{2}}(\beta ,\gamma )=\frac{1}{{C}_{{N}_{1},{N}_{2}}(\gamma )}\sum\limits _{\boldsymbol{s}}^{\prime }{\int }_{0}^{\infty }\mathrm{d}L\enspace {\mathrm{e}}^{-\gamma L}\underset{0< {x}_{1}< \dots < {x}_{N}< L}{\int {\mathrm{d}}^{N}x}\enspace {\mathrm{e}}^{-\beta {{\Phi}}_{N}(\boldsymbol{s},\boldsymbol{x})}\end{equation} \tag{ 2.5 }$$

is the configuration integral. Here, β ≡ 1/k_B T (k_B and T being the Boltzmann constant and the absolute temperature, respectively) and γ ≡ βp (p being the pressure). In equation (2.4), L_ref is a reference length (introduced to make ${\mathcal{Z}}_{N}^{\text{id}}$ dimensionless) and ${{\Lambda}}_{i}(\beta )\equiv h\sqrt{\beta /2\pi {m}_{i}}$ is the thermal de Broglie wavelength (h being the Planck constant and m_i being the mass of a particle of species i). In equation (2.5), the prime in the summation denotes the constraint (2.1). Note that, by construction, ${\mathcal{Q}}_{{N}_{1},{N}_{2}}=1$ if Φ_N = 0.

Let us make ${\mathcal{Q}}_{{N}_{1},{N}_{2}}$ more explicit. First,

$$\begin{align}{\mathcal{Q}}_{{N}_{1},{N}_{2}}& =\frac{1}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{\boldsymbol{s}}^{\prime }{\int }_{0}^{\infty }\mathrm{d}L\enspace {\mathrm{e}}^{-\gamma L}{\int }_{0}^{L}\mathrm{d}{x}_{1}{\int }_{{x}_{1}}^{L}\mathrm{d}{x}_{2}\dots {\int }_{{x}_{N-1}}^{L}\mathrm{d}{x}_{N}\enspace {\mathrm{e}}^{-\beta {{\Phi}}_{N}(\boldsymbol{s},\boldsymbol{x})}\\ & =\frac{1}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{\boldsymbol{s}}^{\prime }{\int }_{0}^{\infty }\mathrm{d}{x}_{1}{\int }_{{x}_{1}}^{\infty }\mathrm{d}{x}_{2}\dots {\int }_{{x}_{N-1}}^{\infty }\mathrm{d}{x}_{N}{\int }_{{x}_{N}}^{\infty }\mathrm{d}L\enspace {\mathrm{e}}^{-\gamma L-\beta {{\Phi}}_{N}(\boldsymbol{s},\boldsymbol{x})},\end{align} \tag{ 2.6 }$$

where in the second step we have changed the order of integration. Next, we perform the change of variables {x₁, x₂, ..., x_N , L} → {x₁, r₂, ..., r_N , r_N+1}, where (see figure 2)

$$\begin{equation}{r}_{i}\equiv {x}_{i}-{x}_{i-1}\quad (i=2,\dots ,N),\quad {r}_{N+1}\equiv {x}_{1}+L-{x}_{N}.\end{equation} \tag{ 2.7 }$$

Note that $L={\sum }_{\alpha =2}^{N+1}\enspace {r}_{\alpha }$. With this change of variables, equation (2.6) becomes

$$\begin{align}{\mathcal{Q}}_{{N}_{1},{N}_{2}}& =\frac{1}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{\boldsymbol{s}}^{\prime }\left[\prod\limits _{\alpha =2}^{N}{\int }_{0}^{\infty }\mathrm{d}{r}_{\alpha }\enspace {\mathrm{e}}^{-\gamma {r}_{\alpha }-\beta {\phi }_{{s}_{\alpha -1},{s}_{\alpha }}({r}_{\alpha })}\right]\\ & \quad \times {\int }_{0}^{\infty }\mathrm{d}{x}_{1}{\int }_{{x}_{1}}^{\infty }\mathrm{d}{r}_{N+1}{\mathrm{e}}^{-\gamma {r}_{N+1}-\beta \omega {\phi }_{{s}_{N},{s}_{1}}({r}_{N+1})}\\ & =\frac{1}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{\boldsymbol{s}}^{\prime }\left[\prod\limits _{\alpha =2}^{N}{{\Omega}}_{{s}_{\alpha -1},{s}_{\alpha }}(\beta ,\gamma )\right]\left[-\frac{\partial {{\Omega}}_{{s}_{N},{s}_{1}}(\beta \omega ,\gamma )}{\partial \gamma }\right],\end{align} \tag{ 2.8 }$$

where

$$\begin{equation}{{\Omega}}_{ij}(\beta ,\gamma )\equiv {\int }_{0}^{\infty }\mathrm{d}r\enspace {\mathrm{e}}^{-\gamma r-\beta {\phi }_{ij}(r)}.\end{equation} \tag{ 2.9 }$$

Zoom In Zoom Out Reset image size

Henceforth, we particularize to open systems (ω = 0), so that

$$\begin{equation}{\mathcal{Q}}_{{N}_{1},{N}_{2}}=\frac{{\gamma }^{-2}}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{\boldsymbol{s}}^{\prime }\prod\limits _{\alpha =2}^{N}{{\Omega}}_{{s}_{\alpha -1},{s}_{\alpha }}.\end{equation} \tag{ 2.10 }$$

Given a spin configuration s , let us call n_ij ( s ) the number of pairs ij. Thus,

$$\begin{equation}\prod\limits _{\alpha =2}^{N}{{\Omega}}_{{s}_{\alpha -1},{s}_{\alpha }}={{\Omega}}_{11}^{{n}_{11}(\boldsymbol{s})}{{\Omega}}_{22}^{{n}_{22}(\boldsymbol{s})}{{\Omega}}_{12}^{{n}_{12}(\boldsymbol{s})}{{\Omega}}_{21}^{{n}_{21}(\boldsymbol{s})}.\end{equation} \tag{ 2.11 }$$

Obviously, n₁₁ + n₂₂ + n₁₂ + n₂₁ = N − 1. If we call w(n₁₁, n₂₂, n₁₂, n₂₁) the number of spin configurations with n_ij pairs ij, equation (2.10) can be rewritten as

$$\begin{equation}{\mathcal{Q}}_{{N}_{1},{N}_{2}}=\frac{{\gamma }^{-2}}{{C}_{{N}_{1},{N}_{2}}}\sum\limits _{{n}_{11},{n}_{22},{n}_{12},{n}_{21}}w({n}_{11},{n}_{22},{n}_{12},{n}_{21})\enspace {{\Omega}}_{11}^{{n}_{11}}{{\Omega}}_{22}^{{n}_{22}}{{\Omega}}_{12}^{{n}_{12}}{{\Omega}}_{21}^{{n}_{21}}.\end{equation} \tag{ 2.12 }$$

Table 1 shows the possible values of n_ij and w for the simple example of N₁ = 4 and N₂ = 2.

Table 1. Spin configurations s for N₁ = 4 and N₂ = 2, organized according to the number (n_ij ) of pairs ij. The number of spin configurations sharing the same values of n_ij is given by w({n_ij }); analogously, w₁₂(n₁₂) is the number of spin configurations sharing the same n₁₂, regardless of the values of n₁₁, n₂₂, and n₂₁.

In general, the evaluation of the number of combinations w({n_ij }) is quite hard. On the other hand, since in the end we will apply the results to the Kern–Frenkel Janus model [36], we can particularize to the case where ϕ₁₁(r) = ϕ₂₂(r) = ϕ₂₁(r), what implies Ω₁₁ = Ω₂₂ = Ω₂₁, so that equation (2.12) reduces to

$$\begin{equation}{\mathcal{Q}}_{{N}_{1},{N}_{2}}=\frac{{\gamma }^{N-1}}{\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)}\sum\limits _{{n}_{12}=0}^{\mathrm{min}\left\{{N}_{1},{N}_{2}\right\}}\enspace {w}_{12}({n}_{12}){{\Omega}}_{11}^{N-1-{n}_{12}}{{\Omega}}_{12}^{{n}_{12}},\end{equation} \tag{ 2.13 }$$

where w₁₂(n₁₂) stands for the number of spin configurations with n₁₂ pairs 12.

To determine w₁₂(n₁₂), imagine that we enumerate particles of each species i = 1 and 2 from left to right as α_i = 1, ..., N_i . Then, each pair of type 12 can be identified with a label (α₁, α₂). Thus, given a number n₁₂, each compatible spin configuration s is characterized by n₁₂ pairs of the form (α₁, α₂). For example, if N₁ = 4 and N₂ = 2 (table 1), the spin configuration s = {112 121} has n₁₂ = 2 pairs: (α₁, α₂) = (2, 1) and (3, 2), while the spin configuration s = {211 121} has a single n₁₂ pair: (α₁, α₂) = (3, 2). There is a one-to-one correspondence between the n₁₂ pairs of the form (α₁, α₂) and the associated spin configuration s . As a consequence, the number of spin configurations w₁₂(n₁₂) with n₁₂ pairs of type 12 is given by the number of ways of choosing the n₁₂ labels α₁ out of N₁ possible values and the n₁₂ labels α₂ out of N₂ possible values. Therefore,

$$\begin{equation}{w}_{12}({n}_{12})=\left(\genfrac{}{}{0pt}{}{{N}_{1}}{{n}_{12}}\right)\left(\genfrac{}{}{0pt}{}{{N}_{2}}{{n}_{12}}\right).\end{equation} \tag{ 2.14 }$$

As a test of consistency, note that the total number of spin configurations is recovered as ${\sum }_{{n}_{12}=0}^{\mathrm{min}\left\{{N}_{1},{N}_{2}\right\}}\enspace {w}_{12}({n}_{12})=\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)$. Finally, the configuration integral is

$$\begin{equation}{\mathcal{Q}}_{{N}_{1},{N}_{2}}=\frac{{\left(\gamma {{\Omega}}_{11}\right)}^{N-1}}{\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)}{{\Xi}}_{{N}_{1},{N}_{2}},\hspace{25.0pt}{{\Xi}}_{{N}_{1},{N}_{2}}\equiv \sum\limits _{n=0}^{\mathrm{min}\left\{{N}_{1},{N}_{2}\right\}}\enspace {\xi }_{{N}_{1},{N}_{2}}(n),\end{equation} \tag{ 2.15 }$$

where

$$\begin{equation}{\xi }_{{N}_{1},{N}_{2}}(n)\equiv \left(\genfrac{}{}{0pt}{}{{N}_{1}}{n}\right)\left(\genfrac{}{}{0pt}{}{{N}_{2}}{n}\right){\left(1-R\right)}^{-n},\quad R\equiv 1-\frac{{{\Omega}}_{11}}{{{\Omega}}_{12}}.\end{equation} \tag{ 2.16 }$$

Interestingly, ${{\Xi}}_{{N}_{1},{N}_{2}}$ can be formally rewritten in terms of the hypergeometric function:

$$\begin{equation}{{\Xi}}_{{N}_{1},{N}_{2}}{=}_{2}{F}_{1}\left(-{N}_{1},-{N}_{2};1,\frac{1}{1-R}\right).\end{equation} \tag{ 2.17 }$$

2.3. Gibbs free energy, internal energy, and equation of state

The finite-size Gibbs free energy G_N (T, p, x₁) is related to the partition function ${\mathcal{Z}}_{{N}_{1},{N}_{2}}(\beta ,\gamma )$ as ${G}_{N}=-{k}_{\text{B}}T\enspace \mathrm{ln}\enspace {\mathcal{Z}}_{{N}_{1},{N}_{2}}$ [32, 37]. According to equations (2.3), (2.4) and (2.15), the finite-size Gibbs free energy per particle g_N = G_N /N can be decomposed as ${g}_{N}={g}_{N}^{\text{id}}+{g}_{N}^{\text{ex}}$, with

$$\begin{equation}\beta {g}_{N}^{\text{id}}={x}_{1}\enspace \mathrm{ln}\left(\gamma {{\Lambda}}_{1}\right)+{x}_{2}\enspace \mathrm{ln}\left(\gamma {{\Lambda}}_{2}\right)-{N}^{-1}\enspace \mathrm{ln}\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)+{N}^{-1}\enspace \mathrm{ln}\left(\gamma {L}_{\text{ref}}\right),\end{equation} \tag{ 2.18a }$$

$$\begin{equation}\beta {g}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right)\mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)-{N}^{-1}\enspace \mathrm{ln}\frac{{{\Xi}}_{{N}_{1},{N}_{2}}}{\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)}.\end{equation} \tag{ 2.18b }$$

By viewing g_N as a function of β and γ (instead of as a function of T and p), it is easy to obtain the average volume (length) per particle (v_N ) and the excess energy per particle (u_N ) at finite N as

$$\begin{equation}{v}_{N}={\left(\frac{\partial \beta {g}_{N}}{\partial \gamma }\right)}_{\beta }={v}_{N}^{\text{id}}+{v}_{N}^{\text{ex}},\hspace{25.0pt}{u}_{N}={\left(\frac{\partial \beta {g}_{N}}{\partial \beta }\right)}_{\gamma }={u}^{\text{id}}+{u}_{N}^{\text{ex}}.\end{equation} \tag{ 2.19 }$$

From equations (2.18a) and (2.18b), one has

$$\begin{equation}{v}_{N}^{\text{id}}=\frac{1+{N}^{-1}}{\gamma },\hspace{25.0pt}{u}^{\text{id}}=\frac{1}{2\beta },\end{equation} \tag{ 2.20a }$$

$$\begin{equation}{v}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }-{N}^{-1}\frac{\partial \enspace \mathrm{ln}\enspace {{\Xi}}_{{N}_{1},{N}_{2}}}{\partial R}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 2.20b }$$

$$\begin{equation}{u}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }-{N}^{-1}\frac{\partial \enspace \mathrm{ln}\enspace {{\Xi}}_{{N}_{1},{N}_{2}}}{\partial R}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma },\end{equation} \tag{ 2.20c }$$

where, in view of equation (2.17),

$$\begin{equation}\frac{\partial {{\Xi}}_{{N}_{1},{N}_{2}}}{\partial R}={\frac{{N}_{1}{N}_{2}}{{(1-R)}^{2}}\enspace }_{2}{F}_{1}\left(-{N}_{1}+1,-{N}_{2}+1;2,\frac{1}{1-R}\right).\end{equation} \tag{ 2.21 }$$

2.4. Limit N → ∞

Equations (2.18b), (2.20b) and (2.20c) provide the excess quantities for any finite N. It is important to take the limit N → ∞ to obtain the TL expressions and their first finite-N corrections.

In appendix A, it is proved that, for large N at fixed mole fractions,

$$\begin{equation}{{\Xi}}_{{N}_{1},{N}_{2}}\approx \frac{{\mathrm{e}}^{N{\bar{\psi }}_{0}}}{\sqrt{2\pi N{y}_{0}(2-{y}_{0}/{x}_{1}{x}_{2})}},\end{equation} \tag{ 2.22 }$$

where

$$\begin{equation}{\bar{\psi }}_{0}=-{x}_{1}\enspace \mathrm{ln}\left(1-\frac{{y}_{0}}{{x}_{1}}\right)-{x}_{2}\enspace \mathrm{ln}\left(1-\frac{{y}_{0}}{{x}_{2}}\right),\quad {y}_{0}=\frac{1-\sqrt{1-4{x}_{1}{x}_{2}R}}{2R}.\end{equation} \tag{ 2.23 }$$

As a consistency test, note that in the case of equal interactions (R → 0), one has y₀ → x₁ x₂ and ${\bar{\psi }}_{0}\to -{x}_{1}\enspace \mathrm{ln}\enspace {x}_{1}-{x}_{2}\enspace \mathrm{ln}\enspace {x}_{2}$, so that ${{\Xi}}_{{N}_{1},{N}_{2}}\to {({x}_{1}^{{N}_{1}}{x}_{2}^{{N}_{2}}\sqrt{2\pi N{x}_{1}{x}_{2}})}^{-1}$. The latter expression is not but the Stirling approximation of $\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)$, as it should be.

Thus, from equation (2.18b) we obtain

$$\begin{equation}\beta {g}_{N}^{\text{ex}}\approx \beta {g}_{\text{TL}}^{\text{ex}}+{N}^{-1}\enspace \mathrm{ln}\left[\gamma {{\Omega}}_{11}\sqrt{(2-{y}_{0}/{x}_{1}{x}_{2}){y}_{0}/{x}_{1}{x}_{2}}\right],\end{equation} \tag{ 2.24 }$$

where

$$\begin{equation}\beta {g}_{\text{TL}}^{\text{ex}}=-\mathrm{ln}(\gamma {{\Omega}}_{11})-{\bar{\psi }}_{0}-{x}_{1}\enspace \mathrm{ln}\enspace {x}_{1}-{x}_{2}\enspace \mathrm{ln}\enspace {x}_{2}\end{equation} \tag{ 2.25 }$$

and we have taken into account that ${N}^{-1}\enspace \mathrm{ln}\left(\genfrac{}{}{0pt}{}{N}{{N}_{1}}\right)\approx -{x}_{1}\enspace \mathrm{ln}\enspace {x}_{1}-{x}_{2}\enspace \mathrm{ln}\enspace {x}_{2}-{N}^{-1}\enspace \mathrm{ln}\sqrt{2\pi N{x}_{1}{x}_{2}}$. Obviously, ${g}_{\text{TL}}^{\text{ex}}$ is the excess Gibbs free energy per particle in the TL. That quantity was evaluated by a completely independent route in [35] with the result

$$\begin{equation}\beta {g}_{\text{TL}}^{\text{ex}}=-\mathrm{ln}(\gamma {{\Omega}}_{11})-\mathrm{ln}\frac{1+\sqrt{1-4{x}_{1}{x}_{2}R}}{2\sqrt{1-R}}+\vert {x}_{1}-{x}_{2}\vert \enspace \mathrm{ln}\frac{\vert {x}_{1}-{x}_{2}\vert +\sqrt{1-4{x}_{1}{x}_{2}R}}{(\vert {x}_{1}-{x}_{2}\vert +1)\sqrt{1-R}}.\end{equation} \tag{ 2.26 }$$

Taking into account the identity (see appendix B for a proof)

$$\begin{align}{\bar{\psi }}_{0}& =-{x}_{1}\mathrm{ln}\enspace {x}_{1}-{x}_{2}\enspace \mathrm{ln}\enspace {x}_{2}+\mathrm{ln}\frac{1+\sqrt{1-4{x}_{1}{x}_{2}R}}{2\sqrt{1-R}}\\ & \quad -\vert {x}_{1}-{x}_{2}\vert \enspace \mathrm{ln}\frac{\vert {x}_{1}-{x}_{2}\vert +\sqrt{1-4{x}_{1}{x}_{2}R}}{(\vert {x}_{1}-{x}_{2}\vert +1)\sqrt{1-R}},\end{align} \tag{ 2.27 }$$

it is obvious that equations (2.25) and (2.26) are equivalent. Note, however, that equation (2.25) is more compact than equation (2.26).

As for the average volume and internal energy per particle, application of equation (2.19) yields

$$\begin{equation}{v}_{\text{TL}}^{\text{ex}}=-{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }-\frac{{y}_{0}^{3}/{x}_{1}{x}_{2}}{(1-{y}_{0}/{x}_{1})(1-{y}_{0}/{x}_{2})}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 2.28a }$$

$$\begin{equation}{u}_{\text{TL}}^{\text{ex}}=-{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }-\frac{{y}_{0}^{3}/{x}_{1}{x}_{2}}{(1-{y}_{0}/{x}_{1})(1-{y}_{0}/{x}_{2})}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma },\end{equation} \tag{ 2.28b }$$

$$\begin{equation}{v}_{N}^{\text{ex}}-{v}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }+\frac{{N}^{-1}}{2}\frac{(1-{y}_{0}/{x}_{1}{x}_{2}){y}_{0}^{2}/2{x}_{1}{x}_{2}}{{(1-{y}_{0}/2{x}_{1}{x}_{2})}^{2}}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 2.28c }$$

$$\begin{equation}{u}_{N}^{\text{ex}}-{u}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }+\frac{{N}^{-1}}{2}\frac{(1-{y}_{0}/{x}_{1}{x}_{2}){y}_{0}^{2}/2{x}_{1}{x}_{2}}{{(1-{y}_{0}/2{x}_{1}{x}_{2})}^{2}}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma }.\end{equation} \tag{ 2.28d }$$

Note that, while u^id has no finite-N contribution, this is not so for ${v}_{N}^{\text{id}}$. According to equation (2.20a), ${v}_{N}^{\text{id}}={v}_{\text{TL}}^{\text{id}}+{\left(\gamma N\right)}^{-1}$, with ${v}_{\text{TL}}^{\text{id}}={\gamma }^{-1}$.

2.5. Equimolar mixture

In the special case of an equimolar binary mixture (${x}_{1}={x}_{2}=\frac{1}{2}$), equations (2.25), (2.28a) and (2.28b) become

$$\begin{equation}\beta {g}_{\text{TL}}^{\text{ex}}=-\mathrm{ln}\left[\frac{\gamma {{\Omega}}_{11}}{2}\left(1+\frac{1}{\sqrt{1-R}}\right)\right],\end{equation} \tag{ 2.29a }$$

$$\begin{equation}{v}_{\text{TL}}^{\text{ex}}=-{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }-\frac{1-\sqrt{1-R}}{2R(1-R)}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 2.29b }$$

$$\begin{equation}{u}_{\text{TL}}^{\text{ex}}=-{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }-\frac{1-\sqrt{1-R}}{2R(1-R)}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma }.\end{equation} \tag{ 2.29c }$$

Analogously, equations (2.24), (2.28c) and (2.28d) simplify to

$$\begin{equation}{g}_{N}^{\text{ex}}-{g}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}\enspace \mathrm{ln}\left[2\gamma {{\Omega}}_{11}\frac{\left(1-\sqrt{1-R}\right){\left(1-R\right)}^{1/4}}{R}\right],\end{equation} \tag{ 2.30a }$$

$$\begin{equation}{v}_{N}^{\text{ex}}-{v}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }-{N}^{-1}\frac{{\left(1-\sqrt{1-R}\right)}^{2}}{4R(1-R)}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 2.30b }$$

$$\begin{equation}{u}_{N}^{\text{ex}}-{u}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }-{N}^{-1}\frac{{\left(1-\sqrt{1-R}\right)}^{2}}{4R(1-R)}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma }.\end{equation} \tag{ 2.30c }$$

In the case of (unbiased) annealed systems, the total number of particles (N) is fixed but the number of particles (N₁ or N₂) with either spin orientation species is allowed to take any value between 0 and N. Thus, the associated configuration integral is

$$\begin{equation}{\mathcal{Q}}_{N}(\beta ,\gamma )=\frac{1}{{C}_{N}(\gamma )}\sum\limits _{{N}_{1}=0}^{N}\enspace \sum\limits _{\boldsymbol{s}}^{\prime }{\int }_{0}^{\infty }\mathrm{d}L\enspace {\mathrm{e}}^{-\gamma L}\underset{0< {x}_{1}< \dots < {x}_{N}< L}{\int {\mathrm{d}}^{N}x}\enspace {\mathrm{e}}^{-\beta {{\Phi}}_{N}(\boldsymbol{s},\boldsymbol{x})},\end{equation} \tag{ 3.1 }$$

where now ${C}_{N}(\gamma )={\sum }_{{N}_{1}=0}^{N}\enspace {C}_{{N}_{1},{N}_{2}}={2}^{N}{\gamma }^{-(N+1)}$ to guarantee that ${\mathcal{Q}}_{N}=1$ if Φ_N = 0.

By following the same steps as those followed to arrive to equation (2.15), we now get

$$\begin{equation}{\mathcal{Q}}_{N}=\frac{{\left(\gamma {{\Omega}}_{11}\right)}^{N-1}}{{2}^{N}}{{\Xi}}_{N},\quad {{\Xi}}_{N}\equiv \sum\limits _{{N}_{1}=0}^{N}\enspace {{\Xi}}_{{N}_{1},{N}_{2}}.\end{equation} \tag{ 3.2 }$$

Consequently,

$$\begin{equation}\beta {g}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right)\mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)+\mathrm{ln}\enspace 2-{N}^{-1}\enspace \mathrm{ln}\enspace {{\Xi}}_{N},\end{equation} \tag{ 3.3a }$$

$$\begin{equation}{v}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }-{N}^{-1}\frac{\partial \enspace \mathrm{ln}\enspace {{\Xi}}_{N}}{\partial R}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 3.3b }$$

$$\begin{equation}{u}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }-{N}^{-1}\frac{\partial \enspace \mathrm{ln}\enspace {{\Xi}}_{N}}{\partial R}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma },\end{equation} \tag{ 3.3c }$$

where we recall that the quantity R is defined by the second equality in equation (2.16).

In the limit of large N it is proved in appendix C that

$$\begin{equation}{{\Xi}}_{N}\approx {\left(1+\frac{1}{\sqrt{1-R}}\right)}^{N}\frac{1+\sqrt{1-R}}{2}.\end{equation} \tag{ 3.4 }$$

Therefore,

$$\begin{equation}\beta {g}_{N}^{\text{ex}}-\beta {g}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}\enspace \mathrm{ln}\frac{2\gamma {{\Omega}}_{11}}{1+\sqrt{1-R}},\end{equation} \tag{ 3.5a }$$

$$\begin{equation}{v}_{N}^{\text{ex}}-{v}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }+{N}^{-1}\frac{1-\sqrt{1-R}}{2R\sqrt{1-R}}{\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta },\end{equation} \tag{ 3.5b }$$

$$\begin{equation}{u}_{N}^{\text{ex}}-{u}_{\text{TL}}^{\text{ex}}\approx {N}^{-1}{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }+{N}^{-1}\frac{1-\sqrt{1-R}}{2R\sqrt{1-R}}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma },\end{equation} \tag{ 3.5c }$$

where the TL quantities are given by equations (2.29a)–(2.29c).

Comparison between equations (2.30a)–(2.30c) and equations (3.5a)–(3.5c) shows that, although the quenched and annealed systems are equivalent in the TL, they differ in their respective finite-size corrections.

Thus far, except for the constraint to nearest neighbors, the interaction potentials ϕ₁₁(r) and ϕ₁₂(r) are arbitrary. In the special case of isotropic interactions, one has ϕ₁₁(r) = ϕ₁₂(r), so that R = 0. In that case,

$$\begin{equation}{{\Xi}}_{{N}_{1},{N}_{2}}=\left(\genfrac{}{}{0pt}{}{{N}_{1}}{{N}_{2}}\right),\quad {{\Xi}}_{N}={2}^{N},\hspace{25.0pt}{\mathcal{Q}}_{{N}_{1},{N}_{2}}={\mathcal{Q}}_{N}={\left(\gamma {{\Omega}}_{11}\right)}^{N-1},\end{equation} \tag{ 4.1a }$$

$$\begin{equation}\beta {g}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right)\mathrm{ln}\left(\gamma {{\Omega}}_{11}\right),\end{equation} \tag{ 4.1b }$$

$$\begin{equation}{v}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta },\quad {u}_{N}^{\text{ex}}=-\left(1-{N}^{-1}\right){\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }.\end{equation} \tag{ 4.1c }$$

Thus, the finite-size effects become almost trivial if the interactions are isotropic and, of course, no distinction between quenched and annealed systems remains.

The situation becomes much more interesting in the genuine Janus case ϕ₁₁(r) ≠ ϕ₁₂(r). We take now the well-known Kern–Frenkel model [7, 36, 38–41], in which case ϕ₁₁(r) and ϕ₁₂(r) correspond to the hard-rod and square-well potentials, respectively, i.e.

$$\begin{equation}{\phi }_{11}(r)=\begin{cases}\infty ,\quad & r< \sigma ,\\ 0,\quad & r > \sigma ,\end{cases}\hspace{25.0pt}{\phi }_{12}(r)=\begin{cases}\infty ,\quad & r< \sigma ,\\ -{\epsilon},\quad & \sigma < r< \lambda \sigma ,\\ 0,\quad & r > \lambda \sigma ,\end{cases}\end{equation} \tag{ 4.2 }$$

where λ ⩽ 2. Henceforth, we take σ = 1, = 1, and /k_B = 1 as units of length, energy, and temperature, respectively. Therefore,

$$\begin{equation}{{\Omega}}_{11}=\frac{{\mathrm{e}}^{-\gamma }}{\gamma },\hspace{25.0pt}{{\Omega}}_{12}={\mathrm{e}}^{\beta }\frac{{\mathrm{e}}^{-\gamma }}{\gamma }-\left({\mathrm{e}}^{\beta }-1\right)\frac{{\mathrm{e}}^{-\lambda \gamma }}{\gamma },\hspace{25.0pt}R={\left\{1+\frac{1}{\left({\mathrm{e}}^{\beta }-1\right)\left[1-{\mathrm{e}}^{-(\lambda -1)\gamma }\right]}\right\}}^{-1},\end{equation} \tag{ 4.3a }$$

$$\begin{equation}{\left(\frac{\partial \enspace \mathrm{ln}\left(\gamma {{\Omega}}_{11}\right)}{\partial \gamma }\right)}_{\beta }=-1,\hspace{25.0pt}{\left(\frac{\partial \enspace \mathrm{ln}\enspace {{\Omega}}_{11}}{\partial \beta }\right)}_{\gamma }=0,\end{equation} \tag{ 4.3b }$$

$$\begin{align}& {\left(\frac{\partial R}{\partial \gamma }\right)}_{\beta }={\left(1-R\right)}^{2}\left({\mathrm{e}}^{\beta }-1\right)(\lambda -1){\mathrm{e}}^{-(\lambda -1)\gamma },\hspace{25.0pt}{\left(\frac{\partial R}{\partial \beta }\right)}_{\gamma }={\left(1-R\right)}^{2}{\mathrm{e}}^{\beta }\left[1-{\mathrm{e}}^{-(\lambda -1)\gamma }\right].\end{align} \tag{ 4.3c }$$

5.1. Equimolar quenched and unbiased annealed systems

In order to confirm the theoretical results provided by equations (2.20b) and (2.20c) for quenched systems and by equations (3.3b) and (3.3c) for (unbiased) annealed systems, we have performed isothermal–isobaric MC simulations. To make contact between the annealed and quenched results in the TL, we have considered equimolar mixtures (${x}_{1}=\frac{1}{2}$) in the latter case. Moreover, the Kern–Frenkel model (4.2) with λ = 1.2 is chosen. Some technical details about the simulation method are given in appendix D.

Tables 2 and 3 give the MC results of v_N and $-{u}_{N}^{\text{ex}}$, respectively, for p = 0.6, T = 1 and 0.2, and N = 4, 10, 20, and 100. Tables 2 and 3 also include the exact theoretical values given by equations (2.20b) and (3.3b) for v_N and by equations (2.20c) and (3.3c) for $-{u}_{N}^{\text{ex}}$. The deviations from the TL values are displayed in figures 3 and 4, which also include the asymptotic behaviors obtained from equations (2.30b) and (2.30c) for (equimolar) quenched systems and from equations (3.5b) and (3.5c) for (unbiased) annealed systems.

Table 2. Values of the average volume (length) per particle, v_N , in equimolar quenched mixtures and in annealed systems for N = 4, 10, 20, and 100. In all cases, λ = 1.2 and p = 0.6. The TL values are v_TL = 2.6000 and 1.2265 at T = 1 and 0.2, respectively.

	T = 1				T = 0.2
	Quenched		Annealed		Quenched		Annealed
N	Exact	MC	Exact	MC	Exact	MC	Exact	MC
4	2.7658	2.77(2)	2.7819	2.80(2)	1.0502	1.050(3)	1.0547	1.063(4)
10	2.6664	2.68(1)	2.6728	2.69(1)	1.1540	1.150(4)	1.1591	1.152(4)
20	2.6332	2.646(5)	2.6364	2.647(5)	1.1903	1.189(3)	1.1936	1.193(3)
100	2.6067	2.612(8)	2.6073	2.623(8)	1.2194	1.218(2)	1.2200	1.219(1)

Table 3. Absolute values of the excess energy per particle, $-{u}_{N}^{\text{ex}}$, in equimolar quenched mixtures and in annealed systems for N = 4, 10, 20, and 100. In all the cases, λ = 1.2 and p = 0.6. The TL values are $-{u}_{\text{TL}}^{\text{ex}}=0.067\enspace 20$ and 0.4421 at T = 1 and 0.2, respectively.

	T = 1				T = 0.2
	Quenched		Annealed		Quenched		Annealed
N	Exact	MC	Exact	MC	Exact	MC	Exact	MC
4	0.068 15	0.0690(8)	0.051 83	0.0510(6)	0.4820	0.481(2)	0.4635	0.461(2)
10	0.067 52	0.0677(4)	0.061 05	0.0610(4)	0.4664	0.468(3)	0.4453	0.447(2)
20	0.067 35	0.0676(3)	0.064 12	0.0645(3)	0.4539	0.453(2)	0.4402	0.442(2)
100	0.067 23	0.0674(3)	0.066 58	0.0668(3)	0.4441	0.444(2)	0.4416	0.439(2)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can observe from tables 2 and 3 and figures 3 and 4 that the simulations nicely confirm our theoretical results. The differences between quenched and annealed finite-size corrections are much more important for the energy than for the volume. In the latter case, there is a change of the sign of v_N − v_TL when decreasing temperature from T = 1 to T = 0.2. Interestingly, except for the energy at low temperature (T = 0.2), the asymptotic behaviors given by equations (2.30b), (2.30c), (3.5b) and (3.5c) apply very well for any N, including N = 4.

5.2. Biased annealed systems

The MC simulations for annealed systems presented above are unbiased in the sense that, even though the identities of the particles are not fixed and thus the mole fraction x₁ is a fluctuating quantity, no preference to either spin orientation is imposed, so that $\langle {x}_{1}\rangle =\frac{1}{2}$. As a consequence, the unbiased annealed results become equivalent to the equimolar quenched ones in the TL.

On the other hand, it is possible to carry out biased annealed simulations by introducing a parameter $q\ne \frac{1}{2}$ which favors one of the two possible spin orientations (see appendix D). As observed in [35], the average value ⟨x₁⟩ ≡ ⟨x⟩ does not coincide with q, but a natural question arises as to whether or not the inequality ⟨x⟩ ≠ q is a finite-size artifact.

To address that question, we have performed MC simulations for biased annealed systems with q = 0.55, 0.65, 0.75, 0.85, and 0.95. As before, we have fixed λ = 1.2, p = 0.6, and temperatures T = 1 and 0.2. As for the number of particles, the values N = 50 and 200 have been chosen. The results are displayed in figure 5, which shows that the data with N = 50 and 200 practically coincide. Therefore, the property ⟨x⟩ ≠ q (actually, $\frac{1}{2}\leqslant \langle x\rangle \leqslant q$ or $q\leqslant \langle x\rangle \leqslant \frac{1}{2}$) and the dependence ⟨x⟩(q) are robust with respect to N and must hold in the TL. While the derivation of the exact function ⟨x⟩(q) seems to be rather involved and lies outside of the scope of this work, we have constructed a simple heuristic approximation in appendix E. Figure 5 shows that equation (E.7) with a = 10 displays an excellent agreement with the MC data.

Zoom In Zoom Out Reset image size

In the MC simulations for biased annealed systems we have also evaluated the specific volume (v) and the excess internal energy per particle (u^ex). Once the robustness of the relationship ⟨x⟩(q) has been checked, one can take q as a parameter and plot v and u^ex as functions of the mole fraction ⟨x⟩. This is done in figure 6. While in the case T = 1 the mapped range is 0.55 ≲ ⟨x⟩ ≲ 0.94, the range shrinks to 0.51 ≲ ⟨x⟩ ≲ 0.63 if T = 0.2. Again, a very weak influence of N is observed. As a matter of fact, comparison with the exact theoretical results for non-equimolar mixtures in the TL (see equations (2.28a) and (2.28b)) presents a very good agreement. It is worth mentioning that v exhibits a rather weak dependence on the mole fraction, with a local minimum at $\langle x\rangle =\frac{1}{2}$. On the other hand, the excess energy u^ex is much more sensitive to ⟨x⟩, vanishing at ⟨x⟩ = 0 and ⟨x⟩ = 1, as expected.

Zoom In Zoom Out Reset image size

This paper has focused on the study of finite-size effects on the thermodynamic quantities of Janus fluids confined to one-dimensional configurations. Two classes of systems (quenched and annealed) have been considered. In the quenched case, the fraction x_i of particles with a particular face (or spin) orientation is kept fixed. On the other hand, particles can flip their orientations in annealed systems, so that the mole fraction x_i fluctuates around a value $\langle {x}_{i}\rangle =\frac{1}{2}$ (unbiased case, ${q}_{i}=\frac{1}{2}$) or $\langle {x}_{i}\rangle \ne \frac{1}{2}$ (biased case, ${q}_{i}\ne \frac{1}{2}$).

Our study allows us to answer affirmatively the four questions initially posed in section 1:

(i) Can the exact derivation of the Gibbs free energy in the TL (N → ∞) be extended to quenched and/or annealed finite-N systems?

By working on the isothermal–isobaric ensemble with open boundary conditions, we have been able to derive exactly the configuration integral (and hence the Gibbs free energy, the specific volume, and the internal energy) for quenched systems with arbitrary values of number of particles N, mole fraction x₁, temperature T, pressure p, and nearest-neighbor interactions ϕ₁₁ and ϕ₁₂. The results are summarized by equations (2.15)–(2.20c).

The exact results for quenched systems are next exploited to get the finite-size quantities for unbiased annealed systems, as given by equations (3.2)–(3.3c).

(ii) Does the quenched ↔ annealed equivalence break down at finite N?

The exact results referred to in the previous point apply to any finite N. An interesting problem consists in taking the limit N → ∞ in order to obtain well-defined expressions for the thermodynamic quantities in the TL, as well as the first N⁻¹-correction. This is done in appendices A and C, the correction results being given by equations (2.24), (2.28c) and (2.28d) for the quenched case and by equations (3.5a)–(3.5c) for the unbiased annealed case.

The quenched quantities in the TL are provided by equations (2.25), (2.28a) and (2.28b). As proved in appendix B, equation (2.25) is equivalent to (but more compact than) the Gibbs free energy derived in [35] from a completely different method. While in [35] the thermodynamic results were derived directly in the TL from the structural correlation functions, here they have been derived by carefully taking the limit N → ∞ from the configuration integral. The equivalence between both routes reinforces the exact character of the results.

The results for equimolar quenched systems and those for unbiased annealed systems agree in the TL (equations (2.29a)–(2.29c)), but they differ in the first N⁻¹-correction (compare equations (2.30a)–(2.30c) with equations (3.5a)–(3.5c)). Therefore, the quenched ↔ annealed equivalence does break down at finite N.

(iii) Can those theoretical predictions be validated by MC simulations?

The conclusions summarized by the two preceding points apply to any choice of the interaction potentials ϕ₁₁ and ϕ₁₂. In order to validate them by simulations, we have specialized to the Kern–Frenkel model [36], as defined by equation (4.2). MC results have been measured for a well range λ = 1.2, a common pressure p = 0.6, two temperatures (T = 1 and 0.2), and four values of the number of particles (N = 4, 10, 20, and 100). As shown by figures 3 and 4, the agreement is very good. Interestingly, except for the case of the internal energy at T = 0.2, the deviations from the TL values closely follow the N⁻¹ rule even for system sizes as small as N = 4.

(iv) Is the dependence of the average mole fraction ⟨x⟩ on the probability q robust with respect to N in annealed MC simulations for biased situations ($q\ne \frac{1}{2}$)?

The finite-size corrections mentioned above for annealed systems apply to unbiased situations. In particular, in each MC step an attempt to assign the orientation identity i = 1 to a given particle is carried out with a probability $q=\frac{1}{2}$, what results in an average mole fraction $\langle x\rangle =\frac{1}{2}$. The procedure can be extended in a straightforward way to a biased choice $q\ne \frac{1}{2}$, which gives rise to $\langle x\rangle \ne \frac{1}{2}$. The naive expectation would be ⟨x⟩ = q, but preliminary results in [35] showed that either $\frac{1}{2}< \langle x\rangle < q$ or $\frac{1}{2} > \langle x\rangle > q$, depending on whether $q > \frac{1}{2}$ or $q< \frac{1}{2}$, respectively. One might reasonable wonder whether the property ⟨x⟩ ≠ q is a finite-size effect that would disappear in the TL.

However, our MC results provide strong evidence about the robustness of the inequality ⟨x⟩ ≠ q and the dependence of ⟨x⟩ on q (see figure 5). This can be qualitatively explained as follows. In the quenched case, the configuration integral presents a peaked local maximum at N₁ = N/2, i.e. $x=\frac{1}{2}$, as can be seen from equations (2.15), (2.22) and (E.1). For annealed systems, this competes against a weight function w_N (x) exhibiting a peaked local maximum at x = q. The annealed probability density P_N (x) is proportional to the product of both functions and then it has a peaked maximum at an intermediate value x = ⟨x⟩. Based on these arguments, a heuristic approach has been put forward in appendix E. Its theoretical predictions (with a single fitting parameter a = 10 independent of T and q) agree excellently well with MC simulations, as figure 5 shows.

As a bonus of the biased annealed simulations, and given the weak influence of N observed in figure 5, we have compared the measured MC values of volume and energy with the theoretical exact results in the TL as functions of the mole fraction. The results displayed by figure 6 show again an excellent agreement.

To put our findings in a proper context, some of their limitations should be remarked. First, the theoretical results have been obtained for open boundary conditions (ω = 0 in equation (2.2)). As shown by equation (2.8), application of periodic boundary conditions (ω = 1) significantly hampers the quest for an exact treatment at finite N. While the choice of the boundary conditions (open or periodic) becomes irrelevant in the TL, finite-size effects are affected by such a choice.

A second limitation arises from the use of the isothermal–isobaric ensemble rather than the standard canonical ensemble. Of course, the partition function and its associated configuration integral can be formally written in the canonical ensemble [consider equation (2.5) with the integration over L removed], but then it is much more difficult to reduce the problem to a purely combinatorial one at finite N, as happens, however, with equations (2.10)–(2.13). One might believe that it would be possible to get the finite-size Helmholtz free energy from the finite-size Gibbs free energy derived here by means of the conventional Legendre transformation. However, this transformation is justified in the TL only and washes out finite-size effects, as we have checked by comparison with canonical MC simulations (not shown).

Third, we have not addressed in the present paper the problem of deriving the exact relationship between ⟨x⟩ and q in biased annealed systems, even in the TL. The theoretical approach in appendix E is heuristic and depends upon a parameter a whose value must be obtained by a fitting procedure. It would be very interesting to analyze in detail the random walk represented by the annealed MC simulations and derive the dependence ⟨x⟩(q), at least in the TL. However, this goal is outside of the scope of the present work.

The last limitation refers to the choice of the one-dimensional geometry itself. Of course, two- and three-dimensional systems are much more realistic, but the one-dimensional setting, apart from being applicable to single-file confinement situations, has the enormous advantage of allowing for the derivation of nontrivial exact results. For instance, we have explicitly shown in a clean way that the first corrections to the TL values are of order N⁻¹, as usually assumed in the literature to get rid of finite-size effects and extrapolate the simulation data to the TL. Moreover, exact results are utterly important to test simulation methods and/or theoretical approaches that can then be extended to scenarios where exact solutions are absent.

A S acknowledges financial support from Grant PID2020-112936GB-I00 funded by MCIN/AEI/10.13039/501100011033, and from Grants IB20079 and GR18079 funded by Junta de Extremadura (Spain) and by ERDF: A way of making Europe.

In this appendix, we prove that the function ${{\Xi}}_{{N}_{1},{N}_{2}}$ defined in equations (2.15) and (2.16) reduces to equation (2.22) in the limit N → ∞.

First, application of the Stirling approximation $x!\approx \sqrt{2\pi x}{(x/\mathrm{e})}^{x}$ yields

$$\begin{equation}{\xi }_{{N}_{1},{N}_{2}}(n=Ny)\approx \mathrm{exp}\left[N\psi (y)\right],\hspace{25.0pt}\psi (y)={\psi }_{0}(y)+{N}^{-1}{\psi }_{1}(y),\end{equation} \tag{ A.1 }$$

where

$$\begin{equation}{\psi }_{0}(y)=-{x}_{1}\enspace \mathrm{ln}\left(1-\frac{y}{{x}_{1}}\right)-{x}_{2}\enspace \mathrm{ln}\left(1-\frac{y}{{x}_{2}}\right)+y\enspace \mathrm{ln}\frac{({x}_{1}-y)({x}_{2}-y)}{{y}^{2}(1-R)},\end{equation} \tag{ A.2 }$$

$$\begin{equation}{\psi }_{1}(y)=-\mathrm{ln}\left[2\pi Ny\sqrt{\left(1-\frac{y}{{x}_{1}}\right)\left(1-\frac{y}{{x}_{1}}\right)}\right].\end{equation} \tag{ A.3 }$$

Equating to zero the first derivative of ψ(y) with respect to y, one can find that the maximum value of ψ(y) corresponds to

$$\begin{equation}{y}_{\mathrm{max}}\enspace \approx {y}_{0}+{N}^{-1}{y}_{1},\end{equation} \tag{ A.4 }$$

where

$$\begin{equation}{y}_{0}=\frac{1-\sqrt{1-4{x}_{1}{x}_{2}R}}{2R},\hspace{25.0pt}{y}_{1}=-\frac{1+(4{y}_{0}-3){y}_{0}/2{x}_{1}{x}_{2}}{2-{y}_{0}/{x}_{1}{x}_{2}}.\end{equation} \tag{ A.5 }$$

Note that ψ₀'(y₀) = 0 and ${y}_{1}=-{\psi }_{1}^{\prime }({y}_{0})/{\psi }_{0}^{{\prime\prime}}({y}_{0})$, where the second derivative of the ψ₀(y) is

$$\begin{equation}{\psi }_{0}^{{\prime\prime}}(y)=-\frac{2-y/{x}_{1}{x}_{2}}{y(1-y/{x}_{1})(1-y/{x}_{2})}.\end{equation} \tag{ A.6 }$$

Note also that the last term on the right-hand side of equation (A.2) vanishes at y = y₀, so that ${\bar{\psi }}_{0}\equiv {\psi }_{0}({y}_{0})$ is given by equation (2.23).

As a second step, let us expand ψ(y) around y = y_max to get

$$\begin{equation}\psi (y)\approx \psi ({y}_{\mathrm{max}})+\frac{{\psi }^{{\prime\prime}}({y}_{\mathrm{max}})}{2}{\left(y-{y}_{\mathrm{max}}\right)}^{2}.\end{equation} \tag{ A.7 }$$

Next, we replace the sum in ${{\Xi}}_{{N}_{1},{N}_{2}}$ by an integral:

$$\begin{align}{{\Xi}}_{{N}_{1},{N}_{2}}& \approx N{\int }_{-\infty }^{\infty }\mathrm{d}y\enspace {\xi }_{{N}_{1},{N}_{2}}(Ny)\\ & \approx N{\mathrm{e}}^{N\psi ({y}_{\mathrm{max}})}{\int }_{-\infty }^{\infty }\mathrm{d}y\enspace {\mathrm{e}}^{\frac{N{\psi }^{{\prime\prime}}({y}_{\mathrm{max}})}{2}{\left(y-{y}_{\mathrm{max}}\right)}^{2}}\\ & =N{\mathrm{e}}^{N\psi ({y}_{\mathrm{max}})}\sqrt{\frac{2\pi }{-N{\psi }^{{\prime\prime}}({y}_{\mathrm{max}})}},\end{align} \tag{ A.8 }$$

where in the second step use has been made of equation (A.7). Finally, taking into account that ψ (y_max) ≈ ψ₀(y₀) + N⁻¹ ψ₁(y₀) and ${\psi }^{{\prime\prime}}({y}_{\mathrm{max}})\approx {\psi }_{0}^{{\prime\prime}}({y}_{0})$, equation (A.8) becomes

$$\begin{equation}{{\Xi}}_{{N}_{1},{N}_{2}}\approx N{\mathrm{e}}^{N{\psi }_{0}({y}_{0})+{\psi }_{1}({y}_{0})}\sqrt{\frac{2\pi }{-N{\psi }_{0}^{{\prime\prime}}({y}_{0})}}.\end{equation} \tag{ A.9 }$$

Insertion of equations (A.3) and (A.6) into equation (A.9) yields equation (2.22).

While ${\bar{\psi }}_{0}$ is expressed in terms of y₀ (see equation (2.23)), the right-hand side of equation (2.27) is expressed in terms of R. The latter quantity is related to y₀ by the identities

$$\begin{equation}R=\frac{{y}_{0}-{x}_{1}{x}_{2}}{{y}_{0}^{2}},\quad \sqrt{1-4{x}_{1}{x}_{2}R}=\frac{2{x}_{1}{x}_{2}}{{y}_{0}}-1,\quad \sqrt{1-R}=\frac{\sqrt{({x}_{1}-{y}_{0})({x}_{2}-{y}_{0})}}{{y}_{0}},\end{equation} \tag{ B.1 }$$

$$\begin{equation}\frac{1+\sqrt{1-4{x}_{1}{x}_{2}R}}{2\sqrt{1-R}}=\frac{{x}_{1}{x}_{2}}{\sqrt{({x}_{1}-{y}_{0})({x}_{2}-{y}_{0})}},\end{equation} \tag{ B.2 }$$

$$\begin{equation}\frac{\vert {x}_{1}-{x}_{2}\vert +\sqrt{1-4{x}_{1}{x}_{2}R}}{(\vert {x}_{1}-{x}_{2}\vert +1)\sqrt{1-R}}=\frac{{x}_{2}}{{x}_{1}}\sqrt{\frac{{x}_{1}-{y}_{0}}{{x}_{2}-{y}_{0}}},\end{equation} \tag{ B.3 }$$

where, without loss of generality, we have assumed x₁ ⩾ x₂ in equation (B.3).

The right-hand side of equation (2.27) can be rewritten as

$$\begin{align}\mathrm{r}.\mathrm{h}.\mathrm{s}.& =-{x}_{1}\enspace \mathrm{ln}\left[{x}_{1}\frac{2\sqrt{1-R}}{1+\sqrt{1-4{x}_{1}{x}_{2}R}}\frac{{x}_{1}-{x}_{2}+\sqrt{1-4{x}_{1}{x}_{2}R}}{({x}_{1}-{x}_{2}+1)\sqrt{1-R}}\right]\\ & \quad -{x}_{2}\enspace \mathrm{ln}\left[{x}_{2}\frac{2\sqrt{1-R}}{1+\sqrt{1-4{x}_{1}{x}_{2}R}}\frac{({x}_{1}-{x}_{2}+1)\sqrt{1-R}}{{x}_{1}-{x}_{2}+\sqrt{1-4{x}_{1}{x}_{2}R}}\right]\\ & =-{x}_{1}\enspace \mathrm{ln}\left(1-\frac{{y}_{0}}{{x}_{1}}\right)-{x}_{2}\enspace \mathrm{ln}\left(1-\frac{{y}_{0}}{{x}_{2}}\right),\end{align} \tag{ B.4 }$$

where we have made use of equations (B.2) and (B.3). Comparison with equation (2.23) closes the proof of equation (2.27).

The method is analogous to the one followed in appendix A. The quantities ${\bar{\psi }}_{0}$ and y₀ defined in equation (2.23) are functions of the mole fraction x₁. It can be checked that ${\bar{\psi }}_{0}$ presents a maximum at ${x}_{1}=\frac{1}{2}$. Expanding in powers of ${x}_{1}-\frac{1}{2}$,

$$\begin{equation}{\bar{\psi }}_{0}\approx \mathrm{ln}\left(1+\frac{1}{\sqrt{1-R}}\right)-\frac{2}{\sqrt{1-R}}{\left({x}_{1}-\frac{1}{2}\right)}^{2}.\end{equation} \tag{ C.1 }$$

Combination of equations (2.22) and (C.1) yields

$$\begin{equation}{{\Xi}}_{{N}_{1},{N}_{2}}\approx {\left(1+\frac{1}{\sqrt{1-R}}\right)}^{N}\frac{1+\sqrt{1-R}}{\sqrt{2\pi N\sqrt{1-R}}}{\mathrm{e}}^{-2N{\left({x}_{1}-\frac{1}{2}\right)}^{2}/\sqrt{1-R}}.\end{equation} \tag{ C.2 }$$

As a second step, for large N the summation of ${{\Xi}}_{{N}_{1},{N}_{2}}$ over N₁ can be approximated by an integral over x₁:

$$\begin{equation}\sum\limits _{{N}_{1}=0}^{N}\enspace {{\Xi}}_{{N}_{1},{N}_{2}}\approx {\left(1+\frac{1}{\sqrt{1-R}}\right)}^{N}\frac{1+\sqrt{1-R}}{\sqrt{2\pi N\sqrt{1-R}}}N{\int }_{-\infty }^{\infty }\mathrm{d}{x}_{1}\enspace {\mathrm{e}}^{-2N{\left({x}_{1}-\frac{1}{2}\right)}^{2}/\sqrt{1-R}}.\end{equation} \tag{ C.3 }$$

This finally gives equation (3.4).

Since our exact finite-size results are found in the isothermal–isobaric ensemble and the Legendre transform 'washes out' the finite-size effects, we found it necessary to perform our numerical experiments also in the isothermal–isobaric ensemble [42]. Moreover, in order to find agreement with our theoretical exact results, open boundary conditions were used. Of course, only in the TL open and periodic boundary conditions become equivalent.

We performed two kinds of MC experiments, which we label as MCa and MCq for annealed and quenched systems, respectively.

The MCa transition rule consists of single particle MC moves (one MC step), which are the combination of a particle position displacement x_α → x_α + (2η − 1)δ, where η is a pseudo-random number in [0, 1] and δ < σ is the maximum displacement (to be kept fixed during the whole simulation to preserve detailed balance) and a particle assignment to species i = 1, 2 with probability q_i (where q₁ = q and q₂ = 1 − q). Open boundary conditions were enforced by generating a new position until it falls inside the segment x_α ∈ [−L/2, L/2]. According to the Metropolis algorithm [43, 44] the move is accepted with probability ${\mathrm{e}}^{-\beta {\Delta}{{\Phi}}_{N}}$, ΔΦ_N being the change in potential energy due to the combined move. This would be enough in the canonical ensemble, while in the isothermal–isobaric ensemble we also need to perform a volume move. The latter is computationally the most expensive one, since it requires a full energy calculation at each attempt and therefore should be used with a low frequency during the run. We chose 30% for the frequency of the volume move in all our simulations. For the transition and acceptance probability for this volume move, see for example [42].

In contrast to the MCa case, in the MCq simulations the particles are assigned an identity i = 1, 2 with probability x_i = q_i from the start and the species assignment is never changed afterwards. The MCq transition rule consists of single particle MC moves that amount to a particle position displacement with δ > σ (note that this condition may be relieved in dimensions higher than one), which is accepted with probability ${\mathrm{e}}^{-\beta {\Delta}{{\Phi}}_{N}}$, ΔΦ_N being the change in potential energy due to the displacement. Again, in the isothermal–isobaric ensemble we also have the volume move [42].

Notice that we can obtain the same result for quenched systems by using a third simulation strategy that we will call MCaq. The MCaq transition rule consists of single particle MC moves that are the combination of a particle position displacement (with δ > σ), which is accepted with probability ${\mathrm{e}}^{-\beta {\Delta}{{\Phi}}_{N}}$ (where ΔΦ_N is the change in potential energy due to the displacement only), followed by a particle assignment to species i = 1, 2 with probability q_i , which is always accepted and therefore completely disentangled from the displacement move. As before, we also have the volume move [42] in the isothermal–isobaric ensemble.

In all cases we chose δ so to have acceptance ratios as close as possible to $\frac{1}{2}$. The equilibration time for MCa was much longer than for MCq.

Given an observable $\mathcal{O}$, its statistical-mechanical average $\langle \mathcal{O}\rangle$ was evaluated by averaging $\mathcal{O}$ over a sufficiently large number of MC configurations after a sufficiently long equilibration time. The measured observables were the mole fraction $x={N}^{-1}{\sum }_{\alpha =1}^{N}\enspace {\delta }_{{s}_{\alpha },1}$, the specific volume (or reciprocal density) v = L/N, and the excess internal energy per particle u^ex = Φ_N /N.

The statistical error on $\langle \mathcal{O}\rangle$ is as usual given by ${\sigma }_{\langle \mathcal{O}\rangle }=\sqrt{{\sigma }_{\mathcal{O}}^{2}{\tau }_{\mathcal{O}}/M}$, where M is the number of MC steps, ${\sigma }_{\mathcal{O}}^{2}$ is the intrinsic variance of $\mathcal{O}$, and ${\tau }_{\mathcal{O}}$ is the correlation time for the observable $\mathcal{O}$ [44]. The latter quantity depends crucially on the transition rule and has a minimum value equal to 1 if one can move so far in configuration space that successive values become uncorrelated. In general, the number of independent steps which contribute to reducing the error bar is not M but $M/{\tau }_{\mathcal{O}}$. Hence, to determine the true statistical error in the random walk, one needs to estimate the correlation time. To do this, it is very important that the total length of the random walk be much greater than ${\tau }_{\mathcal{O}}$. Otherwise, the result and its error bar will not be reliable. In general, there is no mathematically rigorous procedure to determine ${\tau }_{\mathcal{O}}$, so that usually one must determine it from the random walk itself. It is a good practice occasionally to carry out very long runs to test that the results are well converged. In order to equilibrate the random walk, we generally found it necessary to use 10⁶ MC steps at high temperature (T = 1) and 2 × 10⁷ MC steps at low temperature (T = 0.2), and collect averages over M = 10⁵ MC steps.

From equations (C.2) and (3.4), we have that, for large N, the probability that the mole fraction x₁ lies between x and x + dx in the unbiased annealed system is

$$\begin{equation}{P}_{N}(x)\mathrm{d}x=\frac{1}{{{\Xi}}_{N}}\sum\limits _{{N}_{1}=Nx}^{N(x+\mathrm{d}x)}{{\Xi}}_{{N}_{1},{N}_{2}}\approx \frac{N{{\Xi}}_{Nx,N(1-x)}}{{{\Xi}}_{N}}\mathrm{d}x\approx \frac{{\mathrm{e}}^{-2N{\left(x-\frac{1}{2}\right)}^{2}/\sqrt{1-R}}}{\sqrt{\pi \sqrt{1-R}/2N}}\mathrm{d}x.\end{equation} \tag{ E.1 }$$

Obviously, $\langle x\rangle =\frac{1}{2}$.

Imagine now a biased annealed system where each value of x = N₁/N is weighed with a certain function w_N (x) centered around a value $x=q\ne \frac{1}{2}$. In that case,

$$\begin{equation}{P}_{N}(x)\propto {w}_{N}(x){{\Xi}}_{Nx,N(1-x)},\end{equation} \tag{ E.2 }$$

which, for large N, would be extremely peaked around a value (comprised between $\frac{1}{2}$ and q) that coincides with the average $\langle x\rangle ={\int }_{0}^{1}\mathrm{d}x\enspace x{P}_{N}(x)$. Thus, the value ⟨x⟩ can be determined as the solution to the equation

$$\begin{align}0& =\frac{\partial }{\partial x}\underset{N\to \infty }{\mathrm{lim}}\enspace {N}^{-1}\left[\mathrm{ln}\enspace {w}_{N}(x)+\mathrm{ln}\enspace {{\Xi}}_{Nx,N(1-x)}\right]\\ & =\frac{\partial }{\partial x}\underset{N\to \infty }{\mathrm{lim}}\enspace {N}^{-1}\enspace \mathrm{ln}\enspace {w}_{N}(x)+\frac{\partial {\bar{\psi }}_{0}}{\partial x},\end{align} \tag{ E.3 }$$

where in the second step we have made use of equation (2.22). Note that here, in contrast to equation (E.1), we need to take into account the full dependence of ${\bar{\psi }}_{0}$ on x because the solution to equation (E.3) is not, in general, close to $\frac{1}{2}$. According to equation (2.23),

$$\begin{equation}\frac{\partial {\bar{\psi }}_{0}}{\partial x}=-\mathrm{ln}\left[1-\frac{1-\sqrt{1-4x(1-x)R}}{2xR}\right]+\mathrm{ln}\left[1-\frac{1-\sqrt{1-4x(1-x)R}}{2(1-x)R}\right].\end{equation} \tag{ E.4 }$$

The simplest choice for the weight function w_N (x) is the binomial distribution ${w}_{N}(x)=\left(\genfrac{}{}{0pt}{}{{N}_{\text{eff}}}{{N}_{\text{eff}}x}\right){q}^{{N}_{\text{eff}}x}{(1-q)}^{{N}_{\text{eff}}(1-x)}$, where N_eff ≡ Nb, b being an effective factor accounting for the expected dependence of w_N (x) on the thermodynamic state (T and p). In that case,

$$\begin{equation}\underset{N\to \infty }{\mathrm{lim}}{N}^{-1}\enspace \mathrm{ln}\enspace {w}_{N}(x)=b\left[x\enspace \mathrm{ln}\frac{q}{x}+(1-x)\mathrm{ln}\frac{1-q}{1-x}\right],\end{equation} \tag{ E.5 }$$

$$\begin{equation}\frac{\partial }{\partial x}\underset{N\to \infty }{\mathrm{lim}}{N}^{-1}\enspace \mathrm{ln}\enspace {w}_{N}(x)=b\enspace \mathrm{ln}\frac{q(1-x)}{x(1-q)}.\end{equation} \tag{ E.6 }$$

Therefore, equation (E.3) becomes

$$\begin{align}0& =-\mathrm{ln}\left[1-\frac{1-\sqrt{1-4x(1-x)R}}{2xR}\right]+\mathrm{ln}\left[1-\frac{1-\sqrt{1-4x(1-x)R}}{2(1-x)R}\right]\\ & \quad +a\sqrt{1-R}\enspace \mathrm{ln}\frac{q(1-x)}{x(1-q)},\end{align} \tag{ E.7 }$$

where we have taken $b=a\sqrt{1-R}$, a being a constant to be empirically determined. A simple and yet optimal value is a = 10.