Article pubs.acs.org/JPCB

Re-entrant Phase Behavior in Confined Two-Patch Colloidal Particles S. Sokołowski*,† and Y. V. Kalyuzhnyi‡ †

Department for the Modelling of Physico-Chemical Processes, Faculty of Chemistry, MCS University, 20031 Lublin, Poland Institute for Condensed Matter Physics, Svientsitskoho 1, 79011 Lviv, Ukraine

‡

ABSTRACT: We present a version of density functional approach for the system of patchy colloidal particles confined in slitlike pores with hard walls. Each particle possesses two off-center sites of the types A and B, and in addition to single A−A and B−B bonds, formation of the double A−B−A and B−A−B bonds is allowed. The proposed approach is based on the fundamental measure theory and the second order perturbation theory of Wertheim. For the model in question, a re-entrant phase behavior in a bulk system has been found [Kalyuzhnyi Y. V.; Cummings, P. T., J. Chem. Phys. 2013, 139, 104905] . Our calculations revealed that the re-entrant phase diagrams are also observed in confined systems. The upper critical temperature decreases with the pore width, while the lower critical temperature increases very slightly.

■

INTRODUCTION During the last three decades, substantial effort has been devoted to the investigation of the structure, thermodynamic properties, and phase behavior of uniform, as well as nonuniform, associating fluids (see, for example, refs 1−15 and the references quoted therein). A great majority of the investigations were based on the models that combine hardsphere repulsion (or a weak spherically symmetric repulsive− attractive interaction) with strong, short-ranged orientationally dependent attraction. Models of this type belong to the class of the so-called primitive models of associating fluids and were originally introduced some time ago to study molecular systems.16−20 More recently, the interest in the application of the primitive models of association was shifted toward description of systems involving patchy colloidal particles. Patchy colloidal particles can be modeled, assuming that a set of off-center bondable sites are located on the surface of each sphere.21−27 Similarly, as in the case of associating fluids, theoretical description of structure and thermodynamic properties of systems involving patchy colloids can be obtained from thermodynamic perturbation theory of Wertheim.28 In a great majority of works, the first-order thermodynamic perturbation theory has been employed. This theory allows only for the formation of single bonds between relevant sites. However, recently a generalized version of the Wertheim’s theory that allows for the multiple bonding between sites has also been developed.26,29,31−33 Quite recently a patchy colloidal model with unusual reentrant phase behavior was developed and studied via computer simulations and the first-order thermodynamic perturbation theory of Wertheim.34−37 In accordance with this model, a patchy colloidal particle was a hard sphere with two types of off-center attractive sites. There were two sites of one type and n (e.g., n = 9) sites of the second type. Since the theory was based on the first-order perturbation theory, only the existence © 2014 American Chemical Society

of single bonds between appropriate sites has been allowed. The re-entrant transition is the result of a competition between formation of different aggregates (associates). However, to observe the re-entrant transition, it was necessary to assume the existence of a big number of sites, characterized by appropriate energies of association. An alternative model exhibiting a re-entrant transition was proposed by Kalyuzhnyi and Cummings.33 In their approach, a colloidal particle was a hard sphere with only two off-center sites. Each of these sites, however, could accommodate up to two bonds per site. Therefore, in order to describe the properties of the system, the second-order perturbational theory of Wertheim was used. The sites were of a different type and similarly as in the previous model,34−37 the re-entrant phase behavior appears due to the competition between chain formation and the formation of a network. The Wertheim’s theory of association has also found a wide application in the field of inhomogeneous associating fluids. Although the first attempts to describe nonuniform associating fluids were based on solutions of the integral equations,38−41 a great majority of the results were obtained using density functional approaches.11−13,42−48 In density functional approaches, the grand canonical thermodynamic potential is minimized to obtain an equation for the positional-dependent local density. Perhaps the most elegant version of the density functional approach was proposed by Yu and Wu,48 since their approach for associating hard spheres is entirely based on the fundamental measure theory.49,50 Recently, the approach of Yu and Wu was also used by Gnan et al. to study properties of patchy colloidal particles close to a surface.51 Received: April 18, 2014 Revised: June 12, 2014 Published: July 11, 2014 9076

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

excess term, we employ a perturbation approach, according to which the excess intrinsic Helmholtz energy for a system of inhomogeneous associating hard spheres consists of two terms, the first one describes the hard-sphere reference system and the second one accounts for perturbation:

In this work, we generalize the theory proposed by Kalyuzhnyi and Cummings33 to the case of nonuniform systems. Similarly as in the original work, we assume the existence of two sites and allow for the formation of double bonds between appropriate sites. Since the system exhibits reentrant phase behavior in bulk system, our main interest is to check how the confinement in slitlike pores influences the phase diagrams. The extension of the theory follows the ideas of Yu and Wu,48 which have been proven to be very successful in developing a version of the density functional approach based on the first-order thermodynamic perturbation theory.

Fex /kBT =

MODEL AND THEORY We assume the following form of the interaction potential between a pair5,33 of particles (1)

nα(r) =

where r12 is the magnitude of the vector r12 connecting the centers of molecules 1 and 2, and ωi, where i = 1,2, are the orientations of molecules 1 and 2 relative to the vector r12. In this work, the nonassociative part of the pair energy, unon(r12), is just a hard-sphere potential and the hard-sphere diameter is σ. The association energy, uas(r12,ω1,ω2) is the sum of the energies of interactions between the sites K on the molecule 1 and the sites L on the molecule 2, uas(r12,ω1,ω2) = ∑K,LuKL(r12,ω1,ω2), where uKL(r12 , ω1 , ω2 ) ≡ uKL(r12 , θ1 , θ2) ⎧ (c) (c) (c) ⎪− ε KL , σ < r12 < rKL , θ1 < θKL , θ2 < θLK , =⎨ ⎪ otherwise ⎩ 0,

=

(6)

n 3 ln(1 − n3) n1n2 + 2 1 − n3 36πn32

n23

n ·n − V 1 V 2 + n2 1 − n3 36πn3(1 − n3)2 nV 2 ·nV 2 nV 2 ·nV 2 ln(1 − n3) − n2 2 12πn3 12πn3(1 − n3)2

+ (2)

(7)

where the dependence of all the functions on r has been omitted for the sake of brevity. For a bulk system of density ρb, the excess free energy due to association can be derived from the second-order perturbation theory of Wertheim.33 The resulting general expressions can be found in ref 33 in this work; however, we restrict our attention (c) (c) (c) (c) (c) to a special case when θ(c) AA = θBB, θAB = θBA, θAA = θBB and εAA = εBB. This imposed symmetry of bonds leads to a great simplification of the expression for the bulk excess free-energy contribution due to association, f b,as/kBT, ⎡ ⎤ 1 fb,as /kBT = ρb ⎢2 ln X − X + ρb2 X3E3 + 1⎥ ⎣ ⎦ 2

(8)

Here X is the fraction of particles that are bonded at most at one site. It is evaluated by solving the cubic equation 3 2 ρ E3X3 + ρb E2X2 + X − 1 = 0 2 b

(3) 52

(9)

In the above, the quantities E2 and E3 are respectively connected with the volumes available for single (A − A or B − B) or double (A − B − A or B − A − B) bondings. They are given by

According to density functional methodology, the theory provides an approximation for the grand thermodynamic potential,

∫ d rρ(r)[v(z) − μ]

∫ d r′ρ(r′)wα′(r − r′)

fhs /kBT = −n0 ln(1 − n3) +

In the above, is the binding distance, −εKL is the depth of the potential, the lower indices K and L denote conical squarewell sites and take the values A and B, and θ1 (θ2) is the angle between the line connecting the centers of the particles 1 and 2 and the line connecting the center of the particle 1 (2) and its (c) site K (L). In general, θ(c) KL ≠ θLK, but, of course, εKL = εLK and (c) (c) rKL = rLK. The number of the sites L on particle 2 that can be bonded to the site K on particle 1 depends on the values of the parameters of the potential (2), cf. refs 26,29, and 33. Following Kalyuzhnyi and Cummings,33 we consider the model with θ(c) KK (c) ≤ 30° for K = A,B and 30° ≤ θ(c) AB and θBA ≤ 35.3° In this case, each site of the type A (B) can be bonded either to two sites of the type B (A) or to only one site of the same type A (B). The particles are confined in a slitlike pore of width H, thus the external potential is

Ω=F+

∫ d r′ρ(r′)wα(r − r′), nα′(r)

where the definitions of the weight functions, wα(r − r′) and wα′(r − r′), can be found elsewhere.48−50 The hard-sphere Helmholtz energy density consists of contributions from scalar weighted densities and vector weighted densities,50

r(c) KL

⎧ 0, σ /2 < z < H − σ /2, v (z ) = ⎨ ⎩∞ , otherwise

(5)

where f hs(r) and fas(r) are the free-energy densities due to hardsphere and associative interactions, respectively. The hard-sphere reference system free-energy contribution is evaluated using the modified fundamental measure theory.50 As in the original fundamental measure theory of Rosenfeld,49 the expression for f hs(r) is written in terms of four scalar, nα(r) (α = 0, 1, 2, and 3), and two vector, nα′(r) (α′ = V1 and V2), weighted densities. They are defined as

■

u(r1, r2) = unon(r12) + uas(r12 , ω1 , ω2 )

∫ d r[fhs (r) + fas (r)]

E2 = E2(AA) + E2(AB)

(4)

as a functional of the number density ρ(r) at a given value of the configurational chemical potential μ and the temperature T. The Helmholtz free energy, F, is a sum of an ideal and excess terms, F = Fid + Fex. The exact expression for the ideal contribution is Fid = kBTρ(r){ln[ρ(r)] − 1}. To evaluate the

(10)

with (c) E2(KL) = πσ 2[exp(εKL /kBT ) − 1][rKL − σ] (c) 2 [1 − cos θKL ] yhs (σ )

9077

(11)

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

and E3 = σ 6[exp(εAB/kBT ) − 1]2 yhs2 (σ )ych Γ

(12)

In the above, yhs(σ) is the value of the two-particle cavity correlation function at the contact, r = σ, yhs (σ ) =

2 nb,2σ nb,2 σ2 1 + + 1 − nb,3 4(1 − nb ,3)2 72(1 − nb,3)3

1 + achnb,3 +

(1 − nb,3)

⎧ (c) ⎪1, , σ < r < rAB λ (r ) = ⎨ ⎪ ⎩ 0, otherwise

(19)

⎧ (c) (c) ⎪1, θ1 < θAB and θ2 < θAB , U(ω1 , ω2 ) = ⎨ ⎪ ⎩ 0, otherwise

The expression for the variable ych and constant Γ have been derived in refs 33 and 53; they depend on the limiting value of the angle θ(c) AB. In general ych =

(18)

and (13)

2 bchnb,3 3

⎧ 0, r ≤ σ , eR(r ) = ⎨ ⎩1, r > σ

(20)

and A is the surface area. Similarly as in ref 30, evaluation of the integral 17 for different values of z has been carried out using Monte Carlo integration. Figure 1 shows the results obtained

(14) 53

where ach and bch are constants. The quantities nb,2 and nb,3 are the weighted densities n2 and n3, evaluated for a uniform bulk density ρb. In fact, nb,3 = (π/6)ρbσ3 is the bulk packing fraction, and n2,b = πρbσ2. The method of generalization of eqs 8−14 to the case of nonuniform systems used in this work follows the development of Yu and Wu. 48 Formally, to extend eqs 8−14 to inhomogeneous systems, one might replace nb,3, nb,2, and ρb by n3(r), n2(r) and n0(r), respectively, but this approach would include only the scalar-weighted densities in calculating the associative Helmholtz energy. In order to take into account the vector-weighted densities, one has to introduce the proportional factor ζ(r) = [1 − nV2(r)·nV2(r)]/n22(r) to “correct” the weighted densities n0(r) and n2(r), as well as the squared density n22(r). In other words, we replace nb,2, ρb, and n2b,2 everywhere in eqs 8−14 with n2(r)ζ(r), n0(r)ζ(r), and n22(r)ζ(r), respectively, as well as nb,3 with n3(r). Thus, eqs 8 and 13 become

Figure 1. Examples of the integral Γ(z) (eq 17) for H/σ = 3 (solid lines) and for a single wall (dashed line). The results are for θ(c) AB = 35° (black curves) and 45° (red curve) and for r(c) AB = 1.1σ.

{

fas (r)/kBT = ζ(r)n0(r) 2 ln X(r) − X(r) 1 + [ζ(r)n0(r)]2 X3(r)E3(r) + 1 2

}

(c) (c) for r(c) AB = 1.1σ and for θAB = 35° and θAB = 45°. At the distances (c) z > r(c) AB + σ/2 and for the pores H > 2rAB + σ, the function Γ(z)

(15)

becomes equal to its bulk counterpart. The obtained results of

and

Γ(z) close to the wall can be interpolated (e.g., by splines), and

ζ(r)n2(r)σ 1 yhs (r; σ ) = + 1 − n3(r) 4[1 − n3(r)]2 +

the interpolating formula is next used in the calculations of E3, according to eq 12. We show here the function Γ(z) only for z

ζ(r)n22(r)σ 2 72[1 − n3(r)]3

≤ H/2, since Γ(H/2 − z) = Γ(z). The approximations for fas are ad hoc and one can consider a

(16)

Equation 14 is modified in the same manner. The evaluation of the free energy (15) requires the knowledge of X(r). Similarly, as for uniform systems, this is done by solving eq 12 at each point r. However, one should also put emphasis on the quantity Γ that enters eq 12. This is a geometric quantity that characterizes the bond volume. In the case of nonuniform systems, we have to account for the decrease in bond volume near the wall. Formally (cf. 30),

more sophisticated theory, following the developments of Marshall and Chapman.25,26 However, also the latter approach involves several simplifying assumptions the accuracy of which is, in fact, difficult to be assessed. Therefore, we have decided to propose the computationally simplest approach. Note that a similar approximation to that used here was previously

1 Γ(z) = δ(z 2 − z)e(z1)e(z 2)e(z 3) (4π )3 Aσ 2 λ(r12)λ(r23)eR (r13) U (ω1, ω2 )U (ω2 , ω3)d(1)d(2)d(3) × 2 2 r12 r23 (17)

∫

proposed by Malijevsky et al.54 Once an expression for the intrinsic Helmholtz energy has been constructed, minimization of the grand potential, eq 4 with respect to the local density, δΩ/δρ(r) = 0, leads to the

where e(z) = exp[−v(z)/kT],

Euler−Lagrange equation 9078

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B ⎧ μ − v(r) 1 − ρ(r) = exp⎨ k T k ⎩ B BT ⎪

⎪

[∑ α

∑ α′

∂(fhs (r) + fas (r)) ∂nα(r) ∂(fhs (r) + fas (r)) ∂nα ′(r)

Article

we show an example of the dependence of the excess grand canonical potential, ΔΩ = Ω − Ωb (per unit surface area A), on the chemical potential. Here Ωb is the bulk grand canonical potential, calculated from the bulk counterpart of eq 4. We have also displayed the dependence of Ωb on μ within the bulk phase transition. Part b of Figure 2, however, shows the dependence of the average density inside the pore on the chemical potential. The average density is defined as follows

∫ dr′ wα(r − r′) +

⎫ wα ′(r − r′)]⎬ ⎭ ⎪

⎪

(21)

⟨ρ*⟩ =

The expression for the chemical potential (μ) of an uniform system is given in ref 33.

∫0

H*

dz*ρ*(z*)/H *

(22)

(i.e., we average the local density only over the entire pore). In the above z* = z/σ and ρ(z*) = ρ(z*)σ3. The transition in the confined system takes place at the chemical potential μs higher than the chemical potential at the bulk transition point, μc. The adsorption isotherm exhibits the hysteresis loop and the value of μs is closer to the desorption than to the adsorption branch of the isotherm. The observed phenomenon is called the capillary evaporation and has been already observed in the case of Lennard-Jones fluids confined in pores with hard walls. For μ < μs and, in particular for μ < μb, the adsorption ⟨ρ*⟩ is extremely low. This is the consequence that for μ < μs, the local densities in the pore are lower than the bulk densities of an uniform system in equilibrium with the confined one, cf. Figure 3a, where we have displayed the functions g(z) = ρ(z)/ρb for the pore of H* = 5 and at different values of μ/kBT (i.e, at different bulk system densities). Except for very high bulk densities (ρ*b ), all the density profiles are depleted at the pore walls (cf. Figure 3, panels a and b). This is due to strong associative interactions between the particles. Only at very high values of ρ*b , when the repulsive interactions play a more significant role, the development of high local density maxima at the walls is observed. This effect is quite analogous to that reported previously.48 For μ < μs, the functions g(z) are lower than 1 in the entire pore and attain their maximum values at the pore center, z* = H/2. With μ approaching μs the functions g(z) decrease and become very low at the state conditions close to the transition point. Similarly as in the case of bulk uniform systems, the phase behavior of the confined system results from competition between formation of the clusters due to bonding between different sites (the AB-bonding) and due to bonding between the same sites (the AA and BB bonding). The AB bonding causes the formation of a network with up to four bonds per particle, while the bonding between the same sites destroys this network and leads to the formation of chains with up to two bonds per particle. Following previous works,33 to analyze the structure of the clusters formed, we introduce the fractions of the particles in different bonding states. The symbol Xij(z) denotes the fraction of the particles with i and j times bonded A and B sites, respectively. Expressions for the fractions are

■

RESULTS AND DISCUSSION The solution of the local density equation, eq 21, was carried out by using an iterational procedure. During iterations, the evaluation of the integrals involving the weight functions wα and wα was based on the application of the fast fourier transform technique and the subroutines from the FFTW library.55 In order to solve the cubic eq 9, the NewtonRalphson method was employed. Our calculations were performed on a grid of the size of 0.005σ, and the iterations were continued until the maximum deviation between two consecutive iterations was not lower than 10−9 percent. Following the work of Kalyuzhnyi and Cummings,33 all the (c) (c) (c) results were obtained assuming that θ(c) AA = θBB = 5°, θAB = θBA= (c) (c) (c) 35°, rAB = 1.1σ, and rAA = rBB = 1.01σ. We use the parameters σ and ε = εAB as units of length and energy, respectively, and introduce the ratio of energies k = εAA/ε = εBB/ε. The reduced temperature is defined as usual, T* = kBT/ε, and the reduced bulk density is ρb* = ρbσ3. We already know33 that the model defined above with k = 1.7 exhibits a re-entrant transition in the bulk phase, and the upper and lower critical temperatures are T*c,u ≈ 0.0858 and T*c,l ≈ 0.04277, respectively. It is of interest to study how the confinement of the particles in slitlike pores of the width H* = H/σ modifies the phase behavior of the system. In Figure 2a,

X 00(z) = X2(z),

Figure 2. (a) Example of the excess grand canonical potential, ΔΩ/ HAkBT (black line) and of the bulk grand canonical potential, Ωb/ HAkBT, vs the chemical potential (red line). Points b and s indicate the bulk and surface phase transitions, respectively. (b) Average density of the confined fluid, ⟨ρ*⟩ versus the chemical potential. Vertical green line shows the equilibrium transition. Thin vertical red line indicates the value of the chemical potential at the bulk transition, μc/kBT. The part of the isotherm at the chemical potentials μ/kBT ≤ μc/kBT (marked as blue line) has been multiplied by 1000. All the calculations are for k = 1.7, H* = 5, and at T* = 0.07008.

X10(z) = X 01(z) = ζ(z)n0(z)X3(z)E2(z) + ζ(z)n02(z)X 4(z)E3(z), X 20(z) = X 02(z) =

1 ζ(z)n02(z)X 4(z)E3(z) 2 (23)

and 9079

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

the another one [X12(z)] is also present. Therefore, before the transition, the system contains mainly chains and unbonded particles. This situation is quite analogous as to the bulk system.33 After the transition, the amount of free particles decreases and is almost zero in the entire system. Among all Xij’s, the ratio X11(z) is the highest and at the pore center it assumes the value ≈ 0.63. The last value is close to that in the bulk system.33 Also, after the transition, the ratios X21(z) and X22(z) significantly increase: 2X12(H/2) ≈ 0.32 and X22(H/2) ≈ 0.04. We can conclude that after the transition the particles form chains with some amount of branched chains and some amount of network-forming clusters. Figure 4 shows the results similar to those in Figures 2 and 3, but at lower temperature, T* = 0.05. The bulk transition, as

Figure 3. (a) Left panel: the profiles g(z) at bulk densities before the bulk phase transition, ρb* = 5 × 10−5 (μ/kBT = −0.1114979 × 102, solid line), 1 × 10−4 (μ/kBT = −0.1093745 × 102, dashed line), 1.5 × 10−4 (μ/kBT = −0.1084389 × 102, dotted line), and 1.8 × 10−4 (μ/kBT = −0.1080816 × 102 dash-dotted line). Right panel: the profiles g(z) at bulk densities just before the capillary evaporation, ρb* = 0.298 [μ/kBT = −0.1070723 × 102, dashed line, to make this curve visible, the values of g(z) were multiplied by 500], at the bulk density just after the capillary evaporation, ρ*b = 0.299 (μ/kBT = −0.1070540 × 102, solid line) and at ρb* = 0.8 (μ/kBT = 6.73743, dashed dotted line). (b) The fractions Xij(z). Left panel: the fractions X00(z) (black lines), X10(z) (red lines), and X11(z) (blue lines). Right panel: X20(z) (black lines), X12(z) (red lines), and X22(z) (blue lines). The calculations are for bulk densities just before, ρb* = 0.298 (dashed lines) and just after, ρb* = 0.299 (solid lines) capillary evaporation. For some selected curves, the values of Xij(z) were multiplied by 10 or 100, as indicated in the figure. The calculations are for k = 1.7, H* = 5, and at T* = 0.07008.

Xij(z) =

X i 0 ( z ) X 0j ( z ) X 00(z)

,

i ≠ j, i, j ≠ 0

Figure 4. (a) The same as in Figure 1b, but at T* = 0.05. (b) The functions g(z) for bulk densities before (dashed line) and after (solid line) capillary evaporation. (c and d) The functions Xij(z) for (c) ij = 00, 10, and 11 and (d) ij = 20, 12, and 22. The code of the lines is the same as in Figure 2 (panels a and b, respectively). Dashed lines are for ρ*b = 0.094 (μ/kBT = −0.193674 × 102); solid lines are for ρ*b = 0.296 (μ/kBT = −0.1936476 × 102). The calculations are for k = 1.7, H* = 5, and at T* = 0.05.

well as the capillary evaporation, occur at much lower values of the chemical potential than at T* = 0.07008 and the depletion of the local density profiles, ρ(z), is now more pronounced. Before capillary evaporation, as well as after capillary evaporation, the ratio X00(z) is low and in two of those states only the ratio X11(z) is pronounced. Thus, at T* = 0.05, a tendency to form multiply connected network, as well as branched chains, is much lower than at T* = 0.08, and the system contains almost exclusively chains. Summary of our calculations for the system with k = 1.7 is shown in Figure 5. We display here the phase diagrams for the pores of the width H* = 3, 5, and 10, together with the bulk phase diagram. Panel a shows the dependencies of the difference Δμ = μs − μc on temperature for H* = 3, 5, 10, and additionally, for H* = 8. Panel b illustrates the dependence of Δμ on H* at T* = 0.075. Finally, part c shows the phase diagram in the density−temperature plane. The difference between the values of the chemical potential at the capillary evaporation point and at the bulk transition, Δμ, increases with a decrease of the pore width. For the transition inside wider pores, H* > 10, the chemical potential μs is not very distinct from μc, (i.e., the capillary evaporation occurs close to the bulk transition).

(24)

Note that with the definitions (23) and (24), the limiting ground-state (i.e., at T → 0) behavior of the coefficients Xij(z) remains the same as in the uniform system, cf. equation (25) of ref 33. In Figure 3 (panels c and d), we show the functions Xij(z) calculated along the density profiles from Figure 3b. The results displayed here are at T* = 0.07008. Before the capillary evaporation, the fraction of free particles, X00(z), is relatively high, and at the pore center it is close to 0.1, similarly as in a bulk system (cf. Figure 3 of ref 33). Before the transition, the system contains mainly particles singly bonded at one of the sites [their ratio gives the coefficient 2X10(z)] and particles singly bonded at each of the sites [X11(z)], but a small amount of the particles doubly bonded at one site and singly bonded at 9080

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

attractive potential energy results from the formation of bonds. At higher temperatures, the formation of clusters due to the existence of doubly bondable sites occurs, while at lower temperatures, the formation of chains prevails. The possibility to form double bonds describes the geometric parameter Γ(z) (cf. Figure 1). At the pore walls, this possibility is delimited and, therefore, the effect of a decrease of the upper critical temperature for narrower pores has the same origin as in the case of simple confined fluids. Each particle possesses only two associative sites, and the chains formed by the association can assume a nearly planar configuration. Therefore, the lower critical temperatures are only marginally sensitive to the pore width. The second series of calculations has been carried out for the systems characterized by k = 0. In this case, only the formation of AB bonds is possible. We already know that this system exhibits the gas−liquid-like transition in the bulk, and the bulk critical temperature is ≈ 0.08728.33 In Figure 7a, we have plotted examples of adsorption isotherms, evaluated at T* = 0.08 in two pores of the width H*

Figure 5. (a) Phase diagrams in the Δμ/kBT − T* plane. Pore widths are given in the figure. (b) Dependence of Δμ/kBT on the pore width at T* = 0.075. (c) Phase diagrams in the density−temperature plane for H* = 3, 5, 10 and for the bulk system. All the calculations are for k = 1.7.

For all cases, the capillary evaporation is characterized by lower, (T*l ) and upper (T*u ) critical temperatures, similarly as the bulk phase diagram. The upper critical temperature decreases with a decrease of the pore width. This is particularly well seen for the narrowest pore of H* = 3. In the latter case, the upper critical temperature is close to 0.0825, while the bulk upper critical temperature is T*u ≈ 0.0858. The lower critical temperature, however, marginally increases with decreasing pore width. For H* = 3, it is very close to 0.0440, while in the bulk phase it is equal to T*l ≈ 0.04377. In Figure 6, we show

Figure 7. (a) Adsorption isotherms, (b) local densities, and the fractions Xij(z) [(c) ij = 00, 10, and 11 and (d) ij = 20, 12, and 22] at chemical potential before (ρb* = 0.27, μ/kBT = −0.837774 × 101, dashed lines) and after ρ*b = 0.27, μ/kBT = −0.807775 × 101, solid lines). Thick lines are for H* = 3, and thin lines are for H* = 5. The color codes in (c and d) are the same as in Figure 3 (panels c and d, respectively). The curves that are invisible on the figure scale are omitted. The calculations are for k = 0 and at T* = 0.08.

= 5 and H* = 3. The chemical potential at the equilibrium, μs, is higher in the narrower pore, similarly as in the case of the system investigated previously. Figure 7b displays the density profiles in both pores at chemical potentials before and after the capillary evaporation. The profiles for ratified phase are low, and the profiles for dense adsorbed phase exhibit strong depletion at the walls. Parts c and d of Figure 7 show how the fractions Xij(z) change during capillary evaporation. Before capillary evaporation, the fraction of nonbonded particles [X00(z)] is high, especially at the pore walls. It is higher for the narrower pore. After capillary evaporation, almost no unbonded particles are present inside the pores. Before the transition the fractions, 2X10(z) and X11(z), that indicate the formation of dimers and chains are relatively high, but all remaining fractions Xij are very small. After capillary evaporation, the fractions X00(z) and X10(z) become low. Although after the transition the content of chains increases, the most striking difference in

Figure 6. Adsorption isotherms at temperatures in the vicinity of the upper critical temperatures. All the calculations are for k = 1.7 and H* = 5.

the adsorption isotherms at the temperatures very close to the upper critical temperature in the pore of H* = 5. We see that the temperature T* = 0.0845 is slightly below the upper critical temperature. In fact, at T* = 0.0846, the hysteresis loop vanishes, so the upper critical temperature in that pore is between 0.0845 and 0.0846. In the case of molecules interacting via spherically symmetric attractive−repulsive potentials, the decrease of the critical temperature results from a decrease of the possible number of nearest neighbors of particles and thus the potential energy, due to geometrical confinement. For the model in question, the 9081

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

Within the presented approach the values of upper critical temperature decrease with decreasing pore width, similarly as in the case of simple confined fluids. The influence of the pore width on the lower critical temperature is small, and we observe its very small increase with decreasing values of H*. We attribute this behavior to the fact at lower temperatures the particles form chains that can be located almost in a plane. Therefore, geometrical constraints imposed by the walls have nearly negligible influence on the lower critical temperatures. The modified fundamental measure density functional theory provides accurate free-energy functional for three-dimensional hard spheres with the correct properties of dimensional crossover [i.e., the functional is accurate for hard spheres between narrow plates (two-dimensional systems), and inside narrow cylindrical pores (one-dimensional systems)].56 However, in order to check if the associative free energy correctly, and even qualitatively, describes this crossover, it would be necessary to develop a counterpart of the theory of ref 33 for the strictly two-dimensional system. To our knowledge, up to now, there has not been found experimental single component systems presenting a re-entrant phase behavior, though such a behavior is common in binary molecular mixtures.57 However, soft matter and biology offer several examples of particles arranging themselves into aggregates of different geometry.58 In particular, one can expect that Janus-type colloids,59,60 particles forming chains, branched, and ringlike structures, can exhibit a quite complex phase behavior.34,35 Of course, in order to verify the correctness of theoretical predictions either for bulk and for confined systems, it would profitable to perform computer simulation studies of the model in question. Such a study will be carried out in the future. The approach outlined in the present work opens new perspectives of research. In particular, it would be of interest to check for a possibility of the occurrence of surface phase transitions like wetting and layering transitions in the systems involving two-site hard-sphere colloids in contact with a single, attractive wall. Next, as we know, the presence of additional van der Waals type of attraction between colloidal particles modifies the phase behavior in bulk systems. Attractive van der Waals interactions will also modify the phase diagrams of confined fluids and their dependence on the width of pores. All these problems are currently under study.

comparison with the cases considered previously (Figures 2b and 3, panels c and d) is the formation of a network that is reflected by relatively high values of X12(z) and X22(z). The phase diagram for the systems with k = 0 is shown in Figure 8. Similarly as in the case of k = 1.7, the critical

Figure 8. Phase diagrams in the density−temperature plane for the bulk system and for pores of the widths H* = 5 and 3. Right panel magnifies the gaslike branches of the diagram. The calculations are for k = 0. Examples of the integral Γ(z) (eq 17) for H/σ = 3 (solid lines) and for a single wall (dashed line). The results are for θ(c) AB = 35° (black curves) and 45° (red curve) and for θ(c) AB = 1.1σ.

temperatures for confined systems are the same as for the bulk system. The average densities of liquidlike confined phases (⟨ρ*⟩) are considerably lower than the density in the bulk phase, while the densities of gaslike confined phases are slightly higher than their bulk counterparts. In order to visualize the latter statement in panel b, we have plotted parts of the gaslike branches of the phase diagrams using the logarithmic density scale. The explanation for the decrease of the critical temperature on for narrower poresis the same as that given above for the system with k = 1.7.

■

CONCLUSIONS In this paper, we extended the theory of associating hard-sphere patchy colloidal particles, proposed recently33 to the case of nonuniform systems. Each colloidal particles possesses two doubly bondable sites of different types: A-type and B-type sites. Double bonding between different types of sites leads to the formation of a network, while single bonding between the sites promotes chain formation. The proposed theory can be regarded as a combination of the Fundamental Measure Theory and Wertheim’s TPT1 and TPT2 theory to describe association. The extension of the Wertheim approach to the case of nonuniform systems is carried out by using the method proposed by Yu and Wu.48 We study the phase behavior of the system confined in slitlike pores of different width with hard walls and at different values of the association energy for bonding between particular sites. Similarly as for bulk system, there exist a range of k = εAA/ εAB values, where the system exhibits a re-entrant type phase behavior, with lower and upper critical temperatures. The phase transitions in confined systems take place at chemical potentials higher than the bulk chemical potential at the transition point. Thus, the observed phenomenon can be identified as the capillary evaporation.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS S.S. acknowledges support of the European Union under the PIRSES STCSCMBS Grant 268498.

■

REFERENCES

(1) Nezbeda, I.; Kolafa, J.; Kalyuzhnyi, Y. V. Primitive Model of Water. 2. Theoretical Results for the Structure and Thermodynamic Properties. Mol. Phys. 1989, 68, 143−160. (2) Chapman, W. G. Prediction of the Thermodynamic Properties of Associating Lennard-Jones Fluids: Theory and Simulation. J. Chem. Phys. 1990, 93, 4299−4304.

9082

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

(3) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709−1721. (4) Kalyuzhnyi, Y. V.; Holovko, M. F.; Protsykevytch, I. A. Solution of the Associative Percus-Yevick Approximation for the N-Component Mixture of Dimerizing Hard-Spheres. Chem. Phys. Lett. 1993, 215, 1− 4. (5) Kalyuzhnyi, Y. V.; Cummings, P. T. On the Relation Between the Wertheim’s Two-Density Integral Equation Theory for Associating Fluids and Chandler-Silbey-Ladanyi Integral Equation Theory for SiteSite Molecular Fluids. J. Chem. Phys. 1996, 104, 3325−3328. (6) Lin, C. T.; Kalyuzhnyi, Y. V.; Stell, G. Primitive Models of Chemical Association. III. Totally Flexible Sticky Two-Point Model for Multicomponent Heteronuclear Fixed-Chain-Length Polymerization. J. Chem. Phys. 1998, 108, 6513−6524. (7) Müller, E. A.; Gubbins, K. E. Molecular-Based Equations of State for Associating Fluids: A Review of SAFT and Related Approaches. Ind. Eng. Chem. Res. 2001, 40, 2193−2211. (8) Economou, I. G. Statistical Associating Fluid Theory: A Successful Model for the Calculation of Thermodynamic and Phase Equilibrium Properties of Complex Fluid Mixtures. Ind. Eng. Chem. Res. 2002, 41, 953−962. (9) Vlc̆ek, L.; Slovak, J.; Nezbeda, I. Thermodynamic Perturbation Theory of the Second Order: Implementation for Models with Double-Bonded Sites. Mol. Phys. 2003, 101, 2921−2927. (10) Vlc̆ek, L.; Nezbeda, I. Thermodynamics of Simple Models of Associating Fluids: Primitive Models of Ammonia, Methanol, Ethanol and Water. Mol. Phys. 2004, 102, 771−781. (11) Wu, J. Density Functional Theory for Liquid Structure and Thermodynamics. In Molecular Thermodynamics of Complex Systems, Series: Structure and Bonding; Lu, X., Hu, Y., Eds.; Springer: New York, 2009; Vol. 131, pp 1−74. (12) Bryk, P.; Pizio, O.; Sokołowski, S. Density Functional Theory for Inhomogeneous Associating Chain Fluids. J. Chem. Phys. 2006, 125, 024909-1−024909-10. (13) Malo, B. M.; Pizio, O.; Patrykiejew, O.; Sokołowski, S. Adsorption and Phase Transitions in a Two-Site Associating Lennard-Jones Fluid Confined to Energetically Heterogeneous SlitLike Pores; Application of the Density Functional Method. J. Phys.: Condens. Matter 2001, 13, 1361−1379. (14) Kalyuzhnyi, Y. V.; Iacovella, C. R.; Docherty, H.; Holovko, M.; Cummings, P. T. Network Forming Fluids: Yukawa Square-Well mPoint Model. J. Stat. Phys. 2011, 145, 481−506. (15) Khan, S.; Bhandary, D.; Singh, J. K. Surface Phase Transitions of Multiple-Site Associating Fluids. Mol. Phys. 2012, 110, 1241−1248. (16) Bol, W. Monte-Carlo Simulations of Fluid Systems of Waterlike Molecules. Mol. Phys. 1982, 45, 605−616. (17) Smith, W. R.; Nezbeda, I. A Simple Model for Associated Fluids. J. Chem. Phys. 1984, 81, 3694−3699. (18) Chapman, W. G.; Gubbins, K. E. Theory and Simulation of Associating Liquid-Mixtures. Fluid Phase Equilib. 1986, 29, 337−346. (19) Johnson, J. K.; Gubbins, K. E. Phase-Equilibria for Associating Lennard-Jones Fluids from Theory and Simulation. Mol. Phys. 1992, 77, 1033−1053. (20) Ghonasgi, D.; Chapman, W. G. Theory and Simulation for Associating Fluids with 4 Bonding Sites. Mol. Phys. 1993, 79, 291−311. (21) Kern, N.; Frenkel, D. Fluid-Fluid Coexistence in Colloidal Systems with Short-Ranged Strongly Directional Attraction. J. Chem. Phys. 2003, 118, 9882−9889. (22) Tavares, J. M.; Almarza, N. G.; Telo da Gama, M. M. ThreeDimensional Patchy Lattice Model: Ring Formation and Phase Separation. J. Chem. Phys. 2014, 140, 044905-1−044905-9. (23) Marshall, B. D.; Chapman, W. G.; Telo da Gama, M. M. Classical Density Functional Theory for Associating Fluids in Orienting External Fields. Phys. Rev. E 2013, 88, 060301(R)-1− 060301(R)-4. (24) Dias, C. S.; Araujo, N. A. M.; Telo da Gama, M. M. Mixtures of Functionalized Colloids on Substrates. J. Chem. Phys. 2013, 139, 154903-1−154903-6.

(25) Marshall, B. D.; Chapman, W. G. Thermodynamic Perturbation Theory for Associating Fluids with Small Bond Angles: Effects of Steric Hindrance, Ring Formation, and Double Bonding. Phys. Rev. E 2013, 87, 052307-1−052307-12. (26) Marshall, B. D.; Chapman, W. G. A Density Functional Theory for Patchy Colloids Based on Wertheim’s Association Theory: Beyond the Single Bonding Condition. J. Chem. Phys. 2013, 138, 044901-1− 044901-8. (27) Sciortino, F.; Bianchi, E.; Douglas, J. F.; Tartaglia, P. SelfAssembly of Patchy Particles into Polymer Chains: A Parameter-Free Comparison Between Wertheim Theory and Monte Carlo Simulation. J. Chem. Phys. 2007, 126, 194903-1−194903-10. (28) Wertheim, M. Thermodynamic Perturbation-Theory of Polymerization. J. Chem. Phys. 1987, 87, 7323−7331. (29) Marshall, B. D.; Ballal, D.; Chapman, W. G. Publisher’s Note: Wertheim’s Association Theory Applied to One Site Patchy Colloids: Beyond the Single Bonding Condition (vol 137, 104909, 2012). J. Chem. Phys. 2012, 137, 129902-1. (30) Marshall, B. D.; Ballal, D.; Chapman, W. G. Wertheim’s Association Theory Applied to one Site Patchy Colloids: Beyond the Single Bonding Condition. J. Chem. Phys. 2012, 137, 104909-1− 104909-7. (31) Kalyuzhnyi, Y. V.; Protsykevitch, I. A.; Cummings, P. T. LiquidGas Phase Behavior of Stockmayer Fluid with High Dipolar Moment. Condens. Matter Phys. 2007, 10, 553−562. (32) Y. Kalyuzhnyi, Y. V.; Docherty, H.; Cummings, P. T. Resummed Thermodynamic Perturbation Theory for Central Force Associating Potential. Multi-patch Models. J. Chem. Phys. 2010, 133, 044502; 2011, 135, 014501-1−014501-12. (33) Kalyuzhnyi, Y. V.; Cummings, P. T. Two-Patch Colloidal Model with Re-entrant Phase Behaviour. J. Chem. Phys. 2013, 139, 104905-1− 104905-6. (34) Russo, J.; Tavares, J. M.; Teixeira, P. I. C.; Telo da Gama, M. M.; Sciortino, F. Reentrant Phase Diagram of Network Fluids. Phys. Rev. Lett. 2011, 106, 085703-1−085703-4. (35) Russo, J.; Tavares, J. M.; Teixeira, P. I. C.; Telo da Gama, M. M.; Sciortino, F. Re-entrant Phase Behaviour of Network Fluids: A Patchy Particle Model with Temperature-Dependent Valence. J. Chem. Phys. 2011, 135, 034501-1−034501-13. (36) Tavares, J. M.; Teixeira, P. I. C. Phase Diagrams of Particles with Dissimilar Patches: X-Junctions and Y-Junctions. J. Phys.: Condens. Matter 2012, 24, 284108-1−284108-7. (37) Rovigatti, L.; Tavares, J. M.; Sciortino, F. Self-Assembly in Chains, Rings, and Branches: A Single Component System with Two Critical Points. Phys. Rev. Lett. 2013, 111, 168302-1−168302-5. (38) Holovko, M. F.; Vakarin, E. V.; Duda, Y. Y. The Structure of a Sticky Interface. Chem. Phys. Lett. 1995, 233, 420−423. (39) Trokhymchuk, A.; Pizio, O.; Henderson, D.; Sokołowski, S. The Reference Inhomogeneous Percus-Yevick Approximation for a Dimerizing Fluid Near a Hard Wall. Chem. Phys. Lett. 1996, 262, 33−40. (40) Vakarin, E.; Duda, Y.; Holovko, M. Cooperative Adsorption of Network Forming Fluids onto Crystalline Surfaces: Structure and Connectivity of the Interface. J. Chem. Phys. 1997, 107, 5569−5581. (41) Borówko, M.; Sokołowski, S.; Pizio, O. Nonuniform Associating Fluids. In Computational Methods in Surface and Colloid and Colloid Science; Borówko, M., Ed.; Marcel Deker, Inc.: New York, 2000; pp 167−244. (42) Huerta, A.; Pizio, O.; Bryk, P.; Sokołowski, S. Application of the Density Functional Method to Study Phase Transitions in an Associating Lennard-Jones Fluid Adsorbed in Energetically Heterogeneous Slit-Like Pores. Mol. Phys. 2000, 98, 1859−1869. (43) Huerta, A.; Pizio, O.; Sokołowski, S. Phase Transitions in an Associating, Network-Forming, Lennard-Jones Fluid in Slit-Like Pores. II. Extension of the Density Functional Method. J. Chem. Phys. 2000, 112, 4286−4295. (44) Malo, B. M.; Huerta, A.; Pizio, O.; Sokołowski, S. Phase Behavior of Associating Two- and Four-Bonding Sites Lennard-Jones 9083

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084

The Journal of Physical Chemistry B

Article

Fluid in Contact with Solid Surfaces. J. Phys. Chem. B 2000, 104, 7756−7763. (45) Patrykiejew, A.; Sokołowski, S. Adsorption of Associating Fluids on Solid Surfaces: Wetting Transition from Density Functional Theory. J. Phys. Chem. B 1999, 103, 4466−4473. (46) Segura, C. J.; Vakarin, E. V.; Chapman, W. G.; Holovko, M. A Comparison of Density Functional and Integral Equation Theories vs Monte Carlo Simulations for Hard Sphere Associating Fluids near a Hard Wall. J. Chem. Phys. 1998, 108, 4837−4848. (47) Marshall, B. D.; García-Cuéllar, A. J.; Chapman, W. G. A Perturbation Density Functional Theory for the Competition Between Inter and Intramolecular Association. J. Chem. Phys. 2012, 136, 154103-1−154103-11. (48) Yu, Y. X.; Wu, J. A Fundamental-Measure Theory for Inhomogeneous Associating Fluids. J. Chem. Phys. 2002, 116, 7094− 7103. (49) Rosenfeld, Y. Free-Energy Model for the Inhomogeneous HardSphere Fluid Mixture and Density-Functional Theory of Freezing. Phys. Rev. Lett. 1998, 63, 980−983. (50) Yu, Y. X.; Wu, J. Structures of Hard-Sphere Fluids from a Modified Fundamental-Measure Theory. J. Chem. Phys. 2002, 117, 10156−10164. (51) Gnan, N.; de las Heras, D.; Tavares, J. M.; Telo da Gama, M. M.; Sciortino, F. Properties of Patchy Colloidal Particles Close to a Surface: A Monte Carlo and Density Functional Study. J. Chem. Phys. 2012, 137, 084704-1−084704-10. (52) Lutsko, J. F. Recent Developments in Classical Density Functional Theory. Adv. Chem. Phys. 2010, 144, 1−92. (53) Kalyuzhnyi, Y. V.; Marshall, B. D.; Chapman, W. G.; Cummings, P. T. Second-Order Resummed Thermodynamic Perturbation Theory for Central-Force Associating Potential: Multi-Patch Colloidal Models. J. Chem. Phys. 2013, 139, 044909-1−044909-9. (54) Malijevsky, A.; Bryk, P.; Sokołowski, S. Density Functional Approach for Inhomogeneous Star Polymer Fluids. Phys. Rev. E 2005, 72, 032801-1−032801-4. (55) Documentation. fftw v3.3.4; http://www.fftw.org/ #documentation. (56) Rosenfeld, Y.; Schmidt, M.; Lö wen, H.; Tarazona, P. Fundamental-Measure Free-Energy Density Functional for Hard Spheres: Dimensional Crossover and Freezing. Phys. Rev. E 1997, 55, 4245−4263. (57) Avlund, A. S. Extension of Association Models to Complex Chemicals. Ph.D. Thesis, Center for Energy Resources Engineering Department of Chemical and Biochemical Engineering Technical University of Denmark Kongens Lyngby, 2011 (http://orbit.dtu.dk/ en/projects/extension-of-association-models-to-complexchemicals%2803ead894-8bd2-4cfc-aae8-b62183f9dcf8%29.html). (58) Jones, R. Soft Condensed Matter; Oxford University Press: New York, 2002. (59) Rosenthal, G. Theory and Computer Simulations of Amphiphilic Janus Particles. Ph.D. Thesis, Technical University Berlin, 2012 (http://opus.kobv.de/tuberlin/volltexte/2012/3587/). (60) Sciortino, F.; Giacometti, A.; Pastore, G. A Numerical Study of One-Patch Colloidal Particles: From Square-Well to Janus. Phys. Chem. Chem. Phys. 2010, 12, 11869−11877.

9084

dx.doi.org/10.1021/jp503826p | J. Phys. Chem. B 2014, 118, 9076−9084