Article Cite This: Langmuir 2017, 33, 11366-11376
pubs.acs.org/Langmuir
Phase Behavior of Magnetic Nanocolloids of Different Sizes Suspended in an Apolar Solvent Stefanie M. Wandrei,*,† Dannielle G. McCarthy,‡ and Martin Schoen†,§ †
Downloaded via NEW MEXICO STATE UNIV on July 3, 2018 at 21:48:05 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
Stranski-Laboratorium für Physikalische und Theoretische Chemie, Fakultät für Mathematik und Naturwissenschaften, Technische Universität Berlin, Straße des 17. Juni 115, 10623 Berlin, Germany ‡ Department of Chemistry, Stanford University, 333 Campus Dr., Stanford, California 94305, United States § Department of Chemical and Biomolecular Engineering, Engineering Building I, Box 7905, North Carolina State University, 911 Partners Way, Raleigh, North Carolina 27695, United States ABSTRACT: We employ classical density functional theory (DFT) to investigate the phase behavior and composition of binary mixtures; each compound consists of hard spheres of different sizes with superimposed dispersion attraction. In addition to the dispersion attraction, molecules of one component carry an additional three-dimensional magnetic “spin” where the orientation-dependent spin−spin interaction is accounted for by the Heisenberg model. We are treating the excess free energy using a modified mean-field approximation (second virial coefficient) for the orientation-dependent pair correlation function. Depending on the concentration of the magnetic particles, the strength of the spin−spin coupling, and the size ratio of the particles, the model predicts the formation of ordered (polar) phases in addition to the more conventional gas and isotropic liquid phases. Key features of our model are a particle-size dependent shift of the gas− liquid critical point (critical temperature and density) and a change in the width of the phase diagram. In the near-critical region, the latter can be analyzed quantitativly in terms of an ef fective critical exponent βeff that may differ from the classical critical 1 exponent β = 2 ; the classical value is attained in the immediate vicinity of the critical point as it must. The deviation between βeff and β can be linked to nontrivial composition effects along the phase boundaries.
■
INTRODUCTION In the context of (fluidic) soft matter, systems of considerable interest are often composed of molecules that carry an internal vectorial degree of freedom. This is because, generally speaking, these fluids can form ordered phases under certain thermodynamic conditions. Examples of vectorial internal degrees of freedom are electromagnetic (point) dipoles or classical “spins” in various dimensions.1−3 Among these spin fluids, the Heisenberg fluid ranks prominently.4−12 The Heisenberg fluid can be considered as a minimalistic but sufficiently realistic model to describe magnetic ordering in fluids.4−11 In addition, Heisenberg-like interactions arise in the context of amphiphilic colloids composed of antithetic materials (Janus particles).11,13,14 Compared with spin fluids, anisotropic media composed of dipolar (hard or soft) spheres have been more extensively studied in the past.15−18 These systems serve as models for more complex (e.g., capped) magnetic particles.19,20 As far as the phase behavior of these systems is concerned, numerical methods of modern statistical physics have frequently been employed. These comprise Monte Carlo and molecular dynamics simulations5,6,14,15,18 as well as classical density functional theory (DFT) based upon (approximate) freeenergy functional.9,16,17 In the overwhelming number of previous studies on anisotropic fluids, though not exclusively, single-component © 2017 American Chemical Society
systems have been considered. Despite this restriction, spin and dipolar fluids generally exhibit a rather rich topology of phase diagrams. For example, three generic phase-diagram topologies have been established.4,9,10,14,16,21 If the anisotropic interactions are relatively weak, one observes coexistence between a gas (G) and a polar liquid (P) phase at sufficiently low temperatures T; at higher T, a G phase coexists with an isotropic liquid (I). A line of second-order phase transitions (critical line) separates the I from the P phase ending in a critical end point. For an intermediate strength of the anisotropic interactions, the critical end point is transformed into a tricritical point; in addition to GI and GP phase coexistence, one observes such a coexistence between I and P phases. In this case, a triple point exists between all three fluidic phases. For these weak and intermediate strengths of the anisotropic interactions, the phase diagram terminates in a GI critical point. This is not so in the limit of strong anisotropic interactions. In this case, the GI critical point becomes metastable and is submerged in the twophase region. Special Issue: Tribute to Keith Gubbins, Pioneer in the Theory of Liquids Received: June 9, 2017 Revised: July 28, 2017 Published: August 1, 2017 11366
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Langmuir
Article
■
From this brief description of the three generic phasediagram topologies, one intuitively expects the situation to become even more complex if, instead of one-component fluids, binary mixtures are studied. There is a long-standing interest in binary mixtures. Already, back in 1916, Bramley22 noted that “the study of the physical properties of binary mixtures has engaged the attention for a great many years of a large number of chemists and physicists.” In fact, the investigation of binary mixtures continues to be an active field of research until today.23 Examples of recent research on binary mixtures include studies on colloidal particles dispersed in liquid-crystalline carrier fluids,24 block copolymers in a nonsolvent, 25 ferrofluids consisting of magnetic Fe 3 O 4 suspended in a host fluid,26,27 and entropy-driven demixing in mixtures composed of hard bodies.28,29 Although many studies have examined binary mixtures, relatively few theoretical works exist in which one or both components can form ordered phases. For example, Ising mixtures30 composed of a ferromagnetic and a nonmagnetic compound have been investigated; in other works, mixtures of dipolar spheres of different dipolar coupling strengths31−34 as well as binary mixtures composed of dipolar and apolar spheres have been considered.31,35−37 Properties of these systems and, in particular, their possibility to phase separate have important repercussions for the stability of colloidal suspensions;23 the tendency of binary mixtures of nanoparticles and macromolecules to undergo phase separation is relevant for understanding novel and exciting phenomena such as Casimir forces.38 Because of the timeliness of investigating binary mixtures that are capable of forming ordered liquidlike phases, we are studying here such a mixture where one component consists of Heisenberg particles carrying a magnetic, three-dimensional spin; the other component consists of an apolar fluid composed of hard spheres with superimposed dispersion attraction similar to earlier work by two of us.21 However, unlike in this earlier paper, we consider here the impact of the size of the Heisenberg particles relative to that of the solvent molecules. Varying the size ratio has important consequences for the phase behavior. In particular, the width of the GI two-phase region (and therefore that of the entire phase diagram) is significantly affected. We are discussing this phenomenon in terms of an effective GI critical exponent; the concept of an effective critical exponent was originally introduced by Verschaffelt39 and is still employed in modern theoretical40 and experimental41 studies. As we demonstrate in our work, the effective critical exponent can be linked to highly nontrivial composition changes occurring along the coexistence curve. In this study, we are employing classical DFT in which the three pair correlation functions of the binary mixture are approximated at the level of the second virial coefficient. The advantage of this ansatz is that we obtain an analytic expression for the excess free energy by thermodynamic integration; the disadvantage is that intermolecular correlations are seriously overestimated at low temperatures. However, by invoking this so-called modified mean-field (MMF) approximation for the pair correlation functions, we can eventually derive a set of five equations; simultaneous zeros of these equations correspond to coexisting phases at specific (partial) densities and of a particular degree of (polar) order. The zeros are obtained numerically using a Newton−Raphson scheme.
THE MODEL
We consider a binary mixture of magnetic colloids (component 2) suspended in an apolar solvent (component 1). Molecules of both components are modeled as hard spheres with superimposed (isotropic) dispersion attractions. In general, the molecules of each compound may differ in size. Molecules of component 2 carry an additional three-dimensional, classical “spin” (Heisenberg particles). The interaction potential between a pair of particles in this suspension may therefore be written as (ab) φab(r12, ω1, ω2) = φiso (r12) + φanis(r12, ω1, ω2)δa2δb2 ,
a , b = 1, 2
(1)
where a and b label the respective mixture component, r12 = r1 − r2 is the distance vector connecting the centers of mass of the interacting molecular pair (irrespective of the component to which both particles pertain) located at r1 and r2, respectively, the δ’s are Kronecker symbols, and r12 = |r12|. In eq 1, ω1 and ω2 are sets of Euler angles specifying the orientations of a pair of magnetic nanocolloids. Assuming these to be of uniaxial symmetry, ωi = (ϕi, θi), (i = 1,2) where ϕ and θ refer to azimuthal and polar angles, respectively. The Isotropic Core. Particles of both mixture components have a spherically symmetric core. The interaction between those cores is described by (ab) (ab) φiso (r12) = φLJ(ab)(r12) + φHS (r12)
(2)
where the well-known Lennard−Jones (12,6) potential is given here by the expression
⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ φLJ(ab)(r12) = 4εab⎢⎜ ab ⎟ − ⎜ ab ⎟ ⎥ ⎢⎣⎝ r12 ⎠ ⎝ r12 ⎠ ⎥⎦
(3)
and
⎧∞ , r12 < σab (ab) φHS (r12) = ⎨ ⎩ 0, r12 ≥ σab
(4)
is a hard-sphere potential which we introduce for reasons to be explained below. In eqs 3 and 4, σab and εab are measures of the size of the hard cores and the strength of the interaction between particles belonging to components a and b, respectively. If a ≠ b, σab is taken to be the arithmetic mean of σ11 and σ22. Unlike in our previous paper,21 we consider different volume ratios v ≡ σ113/σ223 where σ11 and σ22 are the respective diameters of the hard cores of solvent molecules and colloids. Moreover, as in our previous study,21 we focus on so-called “pseudomixtures” to limit the dimension of the parameter space governing our model mixtures. These pseudomixtures are characterized by ε11 = ε22 = ε ≠ ε12. Notice that, unlike σ12, ε12 is not computed from the corresponding parameters ε11 and ε22 but is treated as a fixed input parameter chosen according to criteria explained at the end of the Introductory Remarks of the Results section. Anisotropic Attraction. Turning next to the anisotropic contribution φanis, we notice from eq 1 that it depends on r12, ω1, and ω2. For any function of these three variables, it is always possible to expand them according to φanis(r12, ω1, ω2) =
∑ φl1l2l(r12)Φl1l2l(ω1, ω2 , ω) l1l 2l
(5)
where l1, l2, and l are positive, semidefinite integers, ω specifies the orientation of r12 in a space-fixed frame of reference, the set {φl1l2l} are expansion coefficients, and Φl1l2l (ω1, ω2 , ω) ≡ Φl1l2l (·) =
∑ m1m2m
11367
C(l1l 2l ; m1m2m)@l1m1(ω1)@l2m2(ω2)@ *lm(ω)
(6)
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir
dimensionless coupling parameter with ξ = 0 for a (suitably defined) reference system and ξ = 1 for the system of interest, and gab is the orientation dependent pair correlation function. Note that for each and every value of ξ a different pair correlation function is needed. Thus, a suitable ansatz for gab is required. Here, we employ the socalled modified mean-field form of gab which may be cast as
is a rotational invariant.42 In eq 6, C is a Clebsch-Gordan coefficient, pairs of corresponding integers l′ (i.e., l1, l2, or l) and m′ are related through m′ ∈ [−l′, l′], {@l ′ m ′} is the set of spherical harmonics, and the asterisk denotes the complex conjugate. At this stage, we introduce a first approximation by assuming that φanis does not depend on ω. This allows us to recast eq 5 as
∑ (−1)l
φanis(r12 , ω1, ω2) =
l
⎧ 0, r12 < σab ⎪ gab(r12 , ω1, ω2 ; ξ) = ⎨ ⎪ ⎩ exp[− βφab(r12 , ω1, ω2 ; ξ)] r12 ≥ σab
2l + 1 φll0(r12)Pl[u(̂ ω1)· u(̂ ω2)] (4π )3/2 (7)
(12)
indicating that φanis is isotropic for any fixed orientation of a pair of particles of component 2. In eq 7, members of the set {Pl} are Legendre polynomials and û is a unit vector specifying the orientation of a particle. At this stage, we introduce two more assumptions. First, we wish particles of component 2 of our suspension to be polar, in the sense that φanis → − φanis if ω′ = (θ′, ϕ′) → − ω′ = (π − θ′, π + ϕ′) where the prime is used to indicate that the inversion operation should be performed on either ω1 or ω2. This implies that only terms for odd values of l in eq 7 survive. We restrict our model even further by considering in the sum in eq 7 only the leading term (l = 1). Second, we assume the expansion coefficient φ110 to describe dispersion interactions42 similar to the attractive contribution in eq 3. Thus, we finally arrive at ⎛σ ⎞ φanis(r12 , ω1, ω2) = 4εεH⎜ 22 ⎟ u(̂ ω1)· u(̂ ω2) ⎝ r12 ⎠
which becomes exact in the limit of vanishing density. Defining in eq 11 the perturbation potential as (ab) φab(r12 , ω1, ω2 ; ξ) = φHS (r12) + ξφ1(ab)(r12 , ω1, ω2)
taking as φ1 the expression
φ1(ab)(r12 , ω1, ω2) = φLJ(ab)(r12) + φanis(r12 , ω1, ω2)δa2δb2
β Δ- ex = −
To determine properties of nanocolloids suspended in an apolar solvent at thermodynamic equilibrium, we are seeking minima of the grand potential functional
∑ μα ∫ dr dω ρα (r , ω) a
(9)
ρ1(r , ω) =
where - is the free energy functional, ρα is the orientation dependent generic singlet distribution function of component α, and μα is the associated chemical potential. Hence, the expression given in eq 9 may be perceived as a generalized Legendre transform of - . In the absence of external fields, the free-energy functional can be split into several individual contributions according to
∑∫ ab
0
1
dξ
r12 ≥ σab
dr1 dr2
∫ dω1 dω2 ρa (r1, ω1) ρb (r2 , ω2) (15)
ρ1 (16a)
4π
(16b)
where the ρ1 and ρ2 appearing on the right-hand side are the number densities and α is the orientation distribution function. It is normalized according to
(10)
∫ dω α(ω) = 1
(17)
Because of the uniaxial symmetry of the ordered phase, it is possible to expand α in terms of Legendre polynomials. Integrations over orientations can then be carried out analytically as we discuss in some detail in the Appendix. Thus, the excess free energy can eventually be expressed compactly as
β Δ- ex = V
∫ dr1dr2 ∫ dω1dω2φ1(ab)(r12 , ω1, ω2)
× ρa (r1, ω1) ρb (r2 , ω2) gab(r12 , ω1, ω2 ; ξ)
ab
ρ2 (r , ω) = ρ2 α(ω)
where - id is the contribution of an ideal-gas binary mixture, - hs reflects the fact that molecules of both components have spherically symmetric hard cores, and Δ- ex accounts for isotropic and anisotropic “soft” repulsive and attractive interactions between molecules of both components in excess of their “hard” repulsive interactions. The Excess Contribution. We begin by discussiong the excess contribution which may be cast as
β Δ- ex 1 = V 2
∑∫
where fab = exp(φ(ab) 1 ) − 1 is the orientation dependent Mayer ffunction. For suitably chosen values for the sets of model parameters and thermodynamic state variables, the model is capable of forming a gas (G), an isotropic liquid (I), and an ordered (polar) (P) bulk phase of uniaxial symmetry. On account of the homogeneity of all three phases and keeping in mind that the interactions between solvent molecules and those between a solvent molecule and a nanocolloid are purely isotropic, we may decompose the generic singlet distribution functions according to
THEORY
- = - id + - hs + Δ- ex
1 2
× fab (r12 , ω1, ω2)
(8)
where we introduced the dimensionless parameter εH to give us the flexibility to vary the strength of the anisotropic attraction relative to the isotropic one.
Ω[{ρα (r , ω)}] = -[{ρα (r , ω)}] −
(14)
and employing our MMF ansatz given in eq 12, the integration over ξ can be carried out analytically and gives
6
■
(13)
∑ ab
β - ex ab = V
∑ ρa ρb u0(ab) + ρ22 α12u1 ab
(18)
where explicit expressions for the energy parameters u(ab) and u1 are 0 given in eqs 55a and 55b, respectively. Hard-Sphere and Ideal-Gas Contributions. Next, we turn to the hard-sphere contribution - hs in eq 10. Because the spherically symmetric hard cores of particles of both components may generally differ in size, we employ an expression for the free energy of the hardsphere mixture based upon the Boublik−Mansoori−Leland−Carnahan−Starling equation of state.31,46−48 This leads to
(11)
where Σab is shorthand notation for Σa2= 1Σb2= 1. The expression in eq 11 is obtained by thermodynamic integration as explained in detail elsewhere.43 In particular, the relationship between various perturbation theories and thermodynamic integration is discussed in ref 43. However, we emphasize that eq 11 is well-known in the literature44 (see also Chapter 3.4 of ref 45). In eq 11, V denotes volume, ξ is a 11368
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir β - hs ≡ βf hs V ⎡⎛ ζ 3 ⎤ ⎞ ζ23 3ζ1ζ2 ⎥ = ρ⎢⎜ 2 2 − 1⎟ln(1 − ζ3) + + 2 ⎢⎣⎝ ζ0ζ3 − ζ ζ (1 ) ζ ζ − ζ (1 ) ⎠ ⎦ 0 3 0 3 3 ⎥
βΩ = − βP V = − ρ + βf hs − β(ρ1μ1hs + ρ2 μ2hs ) −
ab
−
(19) i + ρ2 σ22 ] (i = 0, ..., 3) and ρ = ρ1 + ρ2. In the limit where ζi ≡ of hard spheres of equal size (i.e., for σ11 = σ22), it is easy to verify that the expression in eq 19 reduces to the well-known expression for f hs that one obtains by integrating the Carnahan−Starling equation of state.49 In addition to the previously discussed excess and hard-sphere contributions to the free energy, there is also an ideal-gas contribution. For the present suspension of colloids in an apolar solvent composed of spherically symmetric particles, the ideal-gas contribution may be cast as
λ(T , ρ2 ) − ρ2 ρ2
⎤
1
∫−1 dx Ψ(x)⎥⎦
(25)
Ψ(x) = exp(− 3ρ2 α1u1x)
(26)
Using these results and the fact that because of eq 25 α̅ is automatically normalized, it is relatively straightforward to show that the entropic term can be cast as
⎡
1
1
⎤
∫−1 dx α̅(x) ln[2α̅(x)] = −ln⎢⎣ 12 ∫−1 dx Ψ(x)⎥⎦ − 2ρ2 α12u1
(20)
where Λ is the thermal de Broglie wavelength, m is the mass of a colloid, and I its moment of inertia. As we are focusing exclusively on thermodynamic states at equilibrium, specific values of these three quantities do not matter. The exponents 3 and 5 in the first and second summands on the right-hand side of eq 20 reflect the fact that solvent molecules possess only three (translational) degrees of freedom whereas the nanocolloids have two additional (rotational) degrees of freedom on account of their spins. Together, the first and second summand on the right-hand side of eq 20 represent the kineticenergy contribution to - id . The third term on the right-hand side of eq 20 is entropic in nature [see eq 43 for a definition of α̅ ]. It accounts for the loss in orientational entropy of the suspension as a whole when the colloids cause the formation of an ordered phase. The extra factor of 2 in the argument has been introduced to make sure that in the isotropic phase the third term on the right-hand side of eq 20 vanishes. This is because 1 in the isotropic phase [see eq 43], α̅ = α0 = 2 . Thermodynamic Stability. Inserting eqs 18−20 into eq 9 allows us to express the grand potential in closed form. Because of eqs 16a and 16b, this general expression turns out to be a function of ρ1 and ρ2 and a f unctional of α̅ . The reader should, however, note that the general expression for Ω applies to all thermodynamic states regardless of whether they are thermodynamically stable, metastable, or unstable. Here, we focus exclusively on thermodynamically stable states satisfying in addition the conditions
(27) Phase Coexistence. After having narrowed the analysis to thermodynamically stable phases in the previous section, we now intend to restrict our considerations even further by focusing exclusively on phase coexistence between thermodynamically stable phases. Consider coexisting phases labeled ′ and ″ where we shall tacitly assume that ρ′ ≤ ρ″; the equality holds at isolated points of the phase diagram such as, for example, critical or tricritial points. Moreover, phase ′′ is assumed to be always isotropic whereas phase ′ may be isotropic or polar (i.e., ordered). From standard thermodynamic considerations, the conditions for coexistence of phases ′ and ″ are given by
T′ = T″
(28a)
P′ = P″
(28b)
μ1′ = μ1″
(28c)
μ2 ′ = μ2″
(28d)
Because of eq 24, eq 28b is equivalent to zeros of the equation s1 = − ρ′ + ρ″ + βf hs′ − βf hs″ − β(ρ1′μ1hs′ + ρ2′μ2hs′) + β(ρ1″μ1hs″ + ρ2″μ2hs″) −
∑ (ρa′ρb′ − ρa″ρb″)u0(ab) + ρ2″ 2 α12u1 ab
(29)
=0
β ⎛ ∂Ω ⎞ ⎟⎟ ⎜⎜ =0 V ⎝ ∂ρ1 ⎠ T ,ρ
(21)
β ⎛ ∂Ω ⎞ ⎟⎟ ⎜⎜ V ⎝ ∂ρ2 ⎠
(22)
Using the expression for βμ1 derived from eq 21, the coexistence condition in eq 28c is equivalent to zeros of the equation ⎛ ρ′ ⎞ s2 = ln⎜⎜ 1 ⎟⎟ + βμ1hs′ − βμ1hs″ + 2u0(11)(ρ1′ − ρ1″) ⎝ ρ1″ ⎠
=0
T , ρ1
⎡1 = − ln⎢ ⎣2
where
1
2
(24) hs
= (∂ f /∂ρi)T,ρj is the hard-sphere where P denotes pressure and chemical potential of component i of the suspension. Solving next eq 23 allows one to express the Lagrangian multiplier as
β - id = ρ1[ln(ρ1Λ3) − 1] + ρ2 [ln(ρ2 Λ5mI ) − 1] V
∫−1 dx α̅(x) ln[2α̅(x)]
ρ22 α12u1 μhs i
π [ρ σ i 6 1 11
+ ρ2
∑ ρa ρb u0(ab)
+ u0(12)(ρ2′ − ρ2″) (30)
=0
β δΩ = λ(T , ρ2 ) V δα̅ (x)
By a similar token, eq 28d is equivalent to zeros of
⎛ ρ′ ⎞ s3 = ln⎜⎜ 2 ⎟⎟ + βμ2hs′ − βμ2hs″ + 2u0(22)(ρ2′ − ρ2″) ⎝ ρ2″ ⎠ ⎡1 1 ⎤ + u0(21)(ρ1′ − ρ1″) + ln⎢ dx Ψ(x)⎥ = 0 ⎣ 2 −1 ⎦
(23)
where λ is a Lagrangian multiplier introduced to guarantee that the solutions of eq 23 automatically satisfy eq 17 and the operator δ indicates a functional derivative. Using the aforementioned general expression for Ω, eqs 21 and 16b can be solved for βμ1 and βμ2, respectively. Inserting these expressions into the general one for Ω then gives
∫
where, of course, of the equation 11369
u(12) 0
=
u(21) 0
(31)
holds. In addition to eqs 29−31, the zero DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir 1
s4 = α1 −
3 ∫−1 dx x Ψ(x) =0 2 ∫ 1 dx Ψ(x) −1
(32)
corresponds to the value of the expansion coefficient α1 [and therefore to that of the order parameter 61, cf., eq 44] at thermodynamic equilibrium. Hence, eqs 29−32 need to be solved simultaneously for the set of unknowns {ρ1′ , ρ1″, ρ2′ , ρ2″, α1} which obviously is impossible in principle because we are dealing with an underdetermined problem. To facilitate a unique solution of eqs 29−32, we therefore introduce
s5 = x 2″ −
ρ2″ ρ1″ + ρ2″
=0 (33)
Figure 1. Phase diagrams in temperature(T)−density(ρ) representation for various volume ratios v; (solid black square) v = 8.000, (solid red circle) v = 1.000, (solid green up triangle) v = 0.512, (solid blue down triangle) v = 0.216, and (purple tilted square) v = 0.008. The linear solid line is the critical line obtained from eq 34. Phase diagrams are obtained for a mole faction x″2 = 0.55 for the polar particles in the denser phase.
and fix the mole fraction of the colloids x″2 in the denser phase (and therefore also the mole fraction x1″ = 1 − x2″ of the solvent in that phase). Alternatively, one could also fix x′2. However, it seems closer to experimental conditions to adjust the initial concentration of nanocolloids in the dense (liquid) phase rather than the composition of the low-density (gas) phase.
■
P phase at sufficiently low temperatures T ≤ Tcep ≃ 1.135 up to the temperature Tcep which demarcates a critical end point. For T > Tcep, one then observes coexistence between a G and an I phase. As T increases, the envelopes of the phase diagram pertaining to G and I phases eventually end at a GI critical point at which both phases become indistinguishable. Phase diagrams of this general topology have been termed “type I” earlier.21 The main effect of varying the particle volume ratio is that, for 0.008 < v ≲ 0.512 ≲ v ≲ 8.000, the two-phase region widens considerably with decreasing v; Tc is affected comparably little by varying v and exhibits a small nonmonotonic dependence on that variable. We observe a decrease of the GI critical density ρc whereas the GI critical temperature Tc exhibits only a small depression as the solvent molecules shrink in size with respect to the polar particles. The shift of ρc is a direct consequence of the expression for - hs in eq 19 which accounts for the different packing of the hard spheres of both mixture components. The critical temperature Tc changes because u(ab) ∝ v, and hence, by 0 varying v, the contribution of the fourth term on the far righthand side of eq 24 varies. Because this term accounts for the attractive mean field exerted on particles of both components, it controls the value of Tc. The shift of Tcep with v also visible in Figure 1 can be rationalized as follows. The line of critical points (critical line) shown in Figure 1 can be obtained from a stability analysis as explained by Range and Klapp,31 later by Giura et al.,11 and by Cattes et al.21 The stability analysis leads to the expression 2 1 ρ2c (T ) = − 3 u1(β) (34)
RESULTS Introductory Remarks. In the following, we shall express all physical quantities in dimensionless (i.e., “reduced”) units. We take σ22 as the basic unit of length, express energy in terms of ε, and temperature in units of ε/kB. Other quantities are given in terms of suitable combinations of these basic units. For instance, density is given in units of σ22−3 and pressure in units of εσ22−3. We are then solving eqs 29−33 iteratively by means of a Newton−Raphson scheme as described in ref 21. To initiate the Newton−Raphson solution of eqs 29−33, it is prudent to begin at sufficiently low T, use as initial guesses ρ1′ = ρ2′ = 10−5, and fix x″2 to a desired value for each phase diagram. Because of this, we take as initial guesses ρ″1 = 1 − x″2 and ρ″2 = x″2 such that ρ″1 + ρ″2 = 1. In addition, we choose 6″1 = 61 = 1.0 where we are dropping the superscript ″ because only the high-density phase can be orientationally ordered. With these initial guesses, the Newton−Raphson method converges in less than 10 iterations. Once convergence has been achieved, we increase the initial temperature by a small increment ΔT = 10−5 and apply the Newton−Raphson algorithm again using the converged solutions from the previous T as new initial guesses. This procedure can be repeated up to a temperature at which the polar phase becomes thermodynamically unstable, that is, where eqs 29−33 can no longer be solved simultaneously. One then reinitiates the Newton−Raphson scheme by setting 61 = 0 and by setting the densities to reasonable values estimated from the low-temperature portion of the phase diagram that has already been obtained. A similar procedure can be applied to phase diagrams which in addition to gasisotropic liquid also exhibit coexistence between isotropic and polar liquid phases. Throughout this work, we employ a spin− spin coupling constant εH = 0.12 and ε12 = 1.20 such that both components of the binary mixtures exhibit a tendency to blend. This choice of the model parameters ε12 and εH is motivated by the observation that for particles of equal size the topology of phase diagrams is richest for these choices (see Figure 8 of ref 21). Phase Diagrams and Concentration Profiles. We begin our presentation of results with the phase diagrams of polar particles dispersed in an apolar solvent. For an intermediate concentration, plots in Figure 1 show that if the particles are of equal size v = 1.000, we observe coexistence between a G and a
for the density along the critical line. The shift of Tcep can then be rationalized because one can verify from eq 55b that u1 is independent of v but is an algebraic function of T because of eq 52b. Moreover, the plots in Figure 1 show that the general topology of the phase diagram changes for sufficiently small v. For example, for v = 0.216, the critical end point vanishes and a tricritical point appears at a density ρtcp = 0.76 and a temperature of Ttcp = 1.33; the location of the tricritical point can also be calculated from eq 34. Thus, in a small part of the phase diagram, one observes (IP) coexistence and a triple point 11370
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir appears at Ttr ≃ 1.26 at which G, I, and P phases coexist. Again, a GI critical point is found for v = 0.216. Phase diagrams of this generic topology are termed “type II”.21 Eventually, for v = 0.008, the GI critical point is submerged in the metastable region of the phase diagram as the plot in Figure 1 reveals. Here, we find coexistence between an isotropic fluid of variable density and a P phase terminating again at a tricritical point. Phase diagrams of this generic topology are referred to as “type III”.21 Because the total densities of coexisting phases vary and because the topology of the phase diagram changes with v, nontrivial effects are anticipated as far as the concentration of polar particles in the participating phases are concerned. Moreover, and because x″2 is a fixed quantity, we focus on the relative mole fraction κ = x2′ /x2″. Starting with the case of v = 8.000, where solvent molecules are bigger than the magnetic particles, we see from Figure 2 that κ < 1.0 for all T ≤ Tcep. This
Changing now to a higher (fixed) mole fraction x″2 of polar particles in the denser phase, we see from Figure 3 that for v =
Figure 3. As in Figure 1, but for x″2 = 0.85; (solid black square) v = 8.000, (solid red circle) v = 3.375, (solid green up triangle) v = 1.000, and (solid blue down triangle) v = 0.216. The solid straight line demarcates a line of critical points obtained from eq 34.
1.000, one observes a phase diagram of type-II topology. Starting from this situation, one eventually obtains phase diagrams of type I by making the polar particles bigger with respect to the solvent molecules (i.e., by increasing v), whereas phase diagrams of topology III are observed if v decreases. The critical density changes for reasons already explained. However, comparing plots in Figure 3 with their counterparts shown in Figure 1 reveals a stronger variation of Tc for the former set of data compared with the latter one. As far as the ratio of mole fractions is concerned, in the case of higher (fixed) concentration of polar particles in the denser phase, plots in Figure 4 reveal very similar trends compared
Figure 2. As in Figure 1, but for the relative mole fraction κ along the phase boundaries shown in Figure 1; (solid black square) v = 8.000, (solid red circle) v = 1.000, (solid green up triangle) v = 0.512, (solid blue down triangle) v = 0.216, and (solid purple tilted square) v = 0.008. In the plot corresponding to v = 0.216, the open blue down triangle refers to IP and the solid blue down triangle to GI phase coexistence, respectively.
indicates that the magnetic particles prefer the P over the G phase and that the G phase is enriched with solvent molecules. For T ≥ Tcep, this situation is inverted. Now, the I phase is depleted of polar particles relative to the G phase as revealed by the fact that κ > 1.0. At T = Tc, κ = 1.0 as it must because G and I phases become indistinguishable. This picture changes if the particles of both components are of equal size (v = 1.000). Here, we observe again κ < 1.0 for all T ≤ Tcep but find κ = 1.0 for all Tcep ≲ T ≤ Tc. If the polar particles become bigger than the solvent molecules (i.e., for v = 0.512), the plot in Figure 2 shows that κ < 1.0 even for temperatures above Tcep, which means that even for GI phase coexistence the I phase is enriched with polar particles compared with the G phase. If the generic topology changes from type I to type II (i.e., for v = 0.216), we observe a discontinuity in the plot of κ in Figure 2. This discontinuity arises at T = Ttr (see Figure 1). Above the triple point, κ consists of two branches. One of these refers to GI coexistence whereas along the other branch I and P phases coexist. In both cases the denser phase (I or P) is depleted of polar particles relative to the less dense one (G or I). Both branches end at the GI critical or IP tricritical points where κ = 1. The discontinuity in the plot of κ versus T in Figure 2 vanishes for v = 0.008 because here the GI critical point is submerged in the two-phase region of the phase diagram presented in Figure 1.
Figure 4. As in Figure 2, but for x2″ = 0.85; (solid black triangle) v = 8.000 (solid red circle) v = 3.375, (solid green up triangle) v = 1.000, and v = 0.216 (solid blue down triangle). In the plot corresponding to v = 1.000, open green up triangle and solid green up triangle refer to IP and GI phase coexistence, respectively.
with the plots for x2″ = 0.55 already presented in Figure 2. A noteworthy exception is the plot for v = 1.000 where a subtle effect can be seen above T ≥ Ttr ≃ 1.37. Compared with the corresponding situation depicted in Figure 2, κ > 1 as far as the branch corresponding to GI coexistence is concerned for v = 1.000. Thus, the G phase is enriched with polar particles relative to the I phase. Along the branch referring to IP coexistence, the situation is different: here, κ < 1 indicating that the P phase is richer in polar particles. The Near-Critical Regime. Revisiting the phase diagrams presented in Figures 1 and 3, one notices that the width of the GI two-phase region in the near-critical regime depends sensitively on the size ratio of the particles of which the binary 11371
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir mixture is composed. To analyze this effect quantitatively, we follow the work of Verschaffelt39 and introduce a (nonuniversal) effective critical exponent
βeff ≡
d ln Δρ d ln τ
(35)
where Δρ ≡ ρ″ − ρ′ along the GI portion of the coexistence curve and τ ≡ 1 − T/Tc is a measure of the “distance” from the GI critical point. Clearly, sufficiently close to the GI critical point (in the limit τ → 0), one has to expect that Δρ ∝ τ β
(36) 1
where β = 2 is the classical (mean-field) exponent.50 This conclusion is based upon the fact that our expression for the excess free energy [cf. eq 12] remains analytic in the nearcritical regime and even directly at the GI critical point.51 Analyzing our data according to eq 35 is numerically rather demanding. This is immediately obvious because eq 35 necessitates that Tc be known with high precision. In principle, Tc can be obtained with arbitrary precision within our DFT approach. However, in practice, numerical problems prevent our solution of eqs 29−33 to go directly to the critical point characterized by Δρ = τ = 0. This is reflected by the fact that around the GI critical point our phase diagrams presented in Figures 1 and 3 exhibit a tiny albeit nonvanishing gap. To circumvent this problem, we use a cubic spline to fit the GI part of the phase diagrams. However, care has to be taken to make sure that in using this approach to extrapolate the highand low-density branches of the GI coexistence curve all the way to the critical point, no artifacts are being generated by the spline. To demonstrate the reliability of our numerical analysis, we show, in Figure 5, a plot of ln Δρ as a function of τ where Tc
Figure 6. Effective critical exponent βeff as a function of the “distance” τ from the GI critical point on a semilogarithmic scale; (solid black square) v = 8.000, (solid light bue tilted square) v = 3.375, (solid red circle) v = 1.000, (solid green up triangle) v = 0.512, and (solid blue down triangle) v = 0.216 (cf. Figure 1).
of ln τ in Figure 6. If particles of both mixture species are of 1 equal size (v = 1.000), we see that βeff = β = 2 if thermodynamic state points in a sufficiently near-critical region are considered. Assuming that this sets the relevant temperature range of the near-critical region irrespective of v, we restrict our analysis to a temperature range 0 < Tc − T ≲ 10−2Tc. On the contrary, if the polar particles are much smaller than the solvent molecules (v = 8.000), βeff is much smaller than the 1 classical mean-field value of β = 2 . Comparing, in Figure 6, the plot for v = 8.000 with the one for a less pronounced volume ratio of v = 3.375, one notices that in the latter case βeff is generally shifted to larger values whereas the shape of both curves looks rather similar. If the size of the polar particles exceeds that of the solvent molecules, we notice that for v = 0.512 the maximum deviation 1 of βeff from β = 2 is smaller compared with the previous two cases. This comports with the fact that in the near-critical regime the widths of the GI portion of the phase diagrams for v = 1.000 and v = 0.512 in Figure 1 look rather similar. Finally, comparing, in Figure 1, the plot for v = 1.000 in the GI nearcritical regime with that for v = 0.216, one realizes that the twophase region has widened considerably suggesting that the envelope of the two-phase region deviates more strongly from that for v = 1.000. However, it is gratifying that, sufficiently close to Tc, βeff 1 assumes the classical mean-field value of β = 2 as one would expect on physical grounds. The precise value of T at which this happens seems to depend on v as plots in Figure 6 suggest. 1 That limτ→ 0 βeff = β = 2 sheds additional credibility on our numerical analysis of the phase behavior in the GI near-critical regime.
Figure 5. Plot of the width of the coexistence curve Δρ as a function of the “distance” τ from the GI critical point on a double-logarithmic scale for v = 0.512 and x2″ = 0.55; (red times sign) DFT results from a numerical solution of eqs 29−33, (solid black line) cubic spline fitted 1 to the (discrete) DFT data, and (solid green line) ln Δρ ∝ 2 ln τ .
is defined as that temperature at which we observe the maximum density of the GI part of the coexistence curve. The plot in Figure 5 shows that the spline provides an excellent representation of the (discrete) DFT results obtained by solving eqs 29−33 simultaneously. Moreover, Figure 5 also shows that the numerical solution of eqs 29−33 allows us to 1 cover the crossover region in which βeff → β = 2 . Thus, we do not expect our extrapolation procedure to generate artificial results in the immediate vicinity of the critical point. Having demonstrated the reliability of our fitting procedure, we show plots of the effective critical exponent βeff as a function 11372
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir Turning next to a larger (fixed) mole fraction of polar particles in the liquid phase, plots, in Figure 7, show that similar
It is then a simple matter to show that Δu0eff (x 2′ , x 2″) = Δu0eff (κ ; x 2″) = [2(1 − κ ) − (1 − κ 2)]x 2″ u0(11) − 2[(1 − κ ) − (1 − κ 2)x 2″]x 2″u0(12) − (1 − κ 2)x 2″2u0(22)
(41)
Clearly, Δueff 0 = 0 for κ = 1 irrespective of x2″, which is a fixed parameter in our treatment anyway. In this limiting case, ueff 0 is the same in the G and I phase and thus β Ω′ β Ω″ − ∝ (ρ′2 − ρ″2 )u0eff V V
(42)
This expression has the same form as the corresponding one in a pure fluid consisting either of the solvent molecules or polar particles. For both of these systems, one can show51 that 1 βeff = β = 2 . Plots in Figures 6b and 7b confirm the above 1
1
analysis. However, whether βeff < β = 2 or βeff > β = 2 outside the regime in which κ = 1 depends on x″2 and v and cannot easily be predicted.
■
DISCUSSION AND CONCLUSIONS In this work, we investigate the phase behavior of a binary mixture in which polar particles carrying a classical, threedimensional spin are dispersed in a solvent composed of apolar molecules. Molecules of both components are generally of different size. We employ classical DFT and derive an analytic expression for the grand potential. To derive this expression, we are treating intermolecular correlations in the low-density (secondvirial coefficient) limit. Packing effects are then solely accounted for by the hard-sphere free energy which we derive from the Boublik−Mansoori−Leland−Carnahan−Starling equation of state. Focusing solely on thermodynamic states at equilibrium that coexist with each other, we obtain five algebraic equations to determine the unknowns ρ′1, ρ′2, ρ″1 , ρ″2 , and α1. The phase diagrams constructed with these solutions reveal three generic topologies that have been observed earlier for pure Heisenberg fluids9,14 and for binary mixtures similar to the ones used in this work but for particles of equal size.21 Here, transformations between the three generic types of phase diagrams are observed dependent on the volume ratio v of polar particles and solvent molecules. Similar transformations were observed by Cattes et al.21 However, in ref 21, changes between the three topologies were controlled by the spin−spin coupling constant εH. This is perhaps not surprising because one can imagine that by varying εH one can create the same effective mean field that one would generate by maintaining εH but varying the volume ratio v and/or the mole fraction x2″. This is an immediate consequence of the mean-field character of the MMF approximation. However, subtle differences arise that can solely be attributed to varying v. One of these is a shift in the critical density controlled by f hs. In addition, the variation of the concentration ratio κ is (besides the attractive mean field) controlled by packing phenomena not accounted for earlier.9,14 Another interesting feature is that the shape of the GI part of the coexistence curve depends in a nontrivial way on v. To quantify this effect, we introduce the concept of an effective critical exponent βeff that deviates from the expected classical 1 mean-field exponent of β = 2 . We rationalize this deviation by
Figure 7. As in Figure 6, but for x″2 = 0.85; (solid black square) v = 8.000, (solid purple tilted square) v = 6.859, (solid light blue down triangle) v = 5.832, (solid red circle) v = 3.375, and (solid green up triangle) v = 1.000.
deviations of βeff from the classical mean-field exponent of 1 β = 2 are observed; unlike in Figure 6 and depending on v, βeff can be larger or smaller than β. In any event, as before in Figure 1 6, βeff → β = 2 as τ → 0 as plots in Figure 7 clearly show. At this point the question arises naturally: What causes the peculiar deviation of βeff from the classical mean-field value of 1 β = 2 ? To rationalize this deviation, we notice from eq 24 that βΩ ∝ V
∑ ρa ρb u0(ab) ab
(37)
which represents the net (attractive) mean field exerted on a reference particle. This net mean field determines GI criticality because along the GI portion of the phase diagram neither the ideal-gas nor the hard-sphere contributions would be sufficient to establish GI phase equilibrium; in addition, α1 = 0 along the GI envelope of the two-phase region. The double sum in the previous expression can easily be rewritten as
∑ ρa ρb u0(ab) = ρ2 [(1 − x2)2 u0(11) + 2(1 − x2)x2u0(12) ab
+ x 22u0(22)] ≡ ρ2 u0eff (x 2)
(38)
where ρ1 = (1 − x2) ρ and ρ2 = x2ρ have also been employed and ueff 0 is an (concentration-dependent) effective mean field. At phase coexistence, we thus have β Ω′ β Ω″ − ∝ ρ′2 u0eff (x 2′) − ρ″2 u0eff (x 2″) (39) V V Because the effective mean fields at mole fractions x2′ and x2″ do not necessarily have to be the same, we introduce
Δu0eff (x 2′ , x 2″) = u0eff (x 2′) − u0eff (x 2″)
(40) 11373
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir
realize that the rotational invariants form a complete orthogonal set of basis functions.43,52 Hence, we can write
introducing an additional, implicitly temperature-dependent (through the temperature dependence of κ) effective mean field Δueff 0 . This additional field vanishes sufficiently close to the GI critical point; at the temperature T where this happens, 1 βeff = β = 2 . Because the shape of the GI coexistence curve is affected by v, one anticipates that also other portions of the phase diagram are affected at least indirectly. This seems plausible because along the envelope of the two-phase region, it is impossible that, for a given ρ, a corresponding T is not defined. In other words, at unique points (ρ,T), different branches of the coexistence curve need to be connected. Thus, if the GI portion of the coexistence curve is affected by v, so will all other portions of the entire phase diagram. In this context, it is interesting to note that for some of the phase diagrams discussed in ref 21, the shape of the GI portion of the phase diagram also looks as if it could be characterized by an effective critical exponent that deviates from the classical mean-field one. However, this effect was missed in that earlier study. The analysis presented in this paper offers a plausible explanation because our argument is based upon the existence of a nonvanishing Δueff 0 which is linked to κ ≠ 1 irrespective of whether or not v ≠ 1. Finally, we wish to emphasize that the variation of the shape of the GI portion of the coexistence curve is solely caused by the size ratio of our particles; this feature, which also manifests itself as a deviation of the effective critical exponent from its classical mean-field value, also has an impact on the rest of the phase diagrams because a physically sensible phase diagram must be represented by a continuous curve. On the contrary, other features can be realized by other changing model parameters. This applies, for example, to the topology of the phase diagram which can be controlled independently by ε12 and/or εH in a similar fashion (see Figure 8 of ref 21). Another example is Tc which could be controlled by varying by ε12 and/ or εH as well as through a variation of v.
f l(labl ) (r12) = 12
f l(labl ) (r12) = 12
2 αl = 2l + 1
12
∫ dω1 dω2 Φl l l(ω1, ω2 , ω) Φ*l ′l ′l′(ω1, ω2 , ω) 12
=
12
2l + 1 δl1l1′δl2l2′δll ′ 4π
(48)
we can conclude that (ab) f 000 (r12) = fab (r12),
a = b = 1, a ≠ b
β Δ- ex ab = −2πρa ρb V
∫σ
∞
2 (ab) dr12 r12 f 000 ,
a = b = 1, a ≠ b
ab
(50)
To derive this expression, we transformed variables according to r1 → r1′ = r1 and r2 → r12 = r1 − r2, performed the one trivial integration over dr′1, expressed dr12 in spherical coordinates, and finally integrated over the orientations ω of r12 analytically. The derivation is slightly more involved in the second case where a = b = 2. Here, the orientation distribution function α̅ arises in the integral over orientations in eq 15. After expanding α̅ according to eq 43, it is possible to perform the integrations over orientations analytically as explained in ref 11. One eventually arrives at
(43)
β Δ- ex 1 22 = − ρ2 2 π V
∞
∑ l=0
1 αl 2 (2l + 1)3/2
∫σ
∞
ab
2 (22) dr12r12 f ll0 (r12)
(51)
l
∫−1 dx α̅(x) Pl(x)
∑ f l(labl ) (r12) Φl l l(ω1, ω2 , ω) l1l 2l
12
(49)
follows trivially. Because of this simple relationship and with the aid of eq 16a, we obtain from eq 15
To obtain closed approximate expressions for the set {f (22) ll0 }, we expand the orientation dependent Mayer f-function in eq 46 in terms of powers of βφanis. It is then again possible to perform the integrations over orientations as we discuss in detail in the Appendix of ref 12. Retaining in the expansion of fab in eq 15 terms up to third order in βφanis and limiting the discussion from now on to the leading terms l = 0, 1 in eq 51, we eventually find
(44)
to quantify order in the various phases. To proceed with the calculation of Δ- ex , we expand the orientation dependent Mayer f-function introduced in eq 15 according to fab (r12 , ω1 , ω2) =
∫ dω1 dω2 Φ*l l l(ω1, ω2 , ω)
where fab now denotes the radial Mayer f-function. Using in the previous expression the fact that
where x ≡ cos θ and the set {αl} are expansion coefficients which allow one to introduce order parameters 0 ≤ 6l ≤ 1 6l =
4π f (r12) 2l + 1 ab
(47)
∞
l=0
12
(46)
APPENDIX: INTEGRATION OVER ORIENTATIONS In this Appendix, we present some mathematical details that allow us to derive eq 18 starting from eq 15. On account of the uniaxial symmetry of our molecules, we begin by expanding the orientation distribution function in terms of Legendre polynomials according to
∑ αlPl(x)
∫ dω1 dω2 fab (r12 , ω1, ω2) Φ*l l l(ω1, ω2 , ω)
The treatment thus far is completely general. However, from now on, we need to distinguish between two different cases. In the first case (a = b = 1 and a ≠ b) the interactions are purely isotropic. However, the second case (a = b = 2) introduces anisotropic contributions to the interactions since the colloids carry a (classical) spin. Let us begin with the simpler of the two cases. Because for a = b = 1 and a ≠ b, φ(ab) = φ(ab) 1 LJ the orientation dependent Mayer f-function in eq 46 can be replaced by its radial counterpart and therefore eq 46 can be rewritten as
■
2πα(ω) = α̅ (x) =
4π 2l + 1
12
(45)
(22) f 000
where the rotational invariants have already been introduced in eq 6. To obtain the expansion coefficients {f (ab) l1l2l }, we have to
(4π )3/2 11374
⎛ β 2A2 ⎞ = exp( −βφLJ(22))⎜1 + ⎟−1 2 ⎠ ⎝
(52a)
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir (22) f110
(4π )3/2
⎛ 3 3 3⎞⎟ = exp( −βφLJ(22))⎜βA + βA ⎝ ⎠ 10
(5) Lomba, E.; Weis, J.-J.; Almarza, N. G.; Bresme, F.; Stell, G. Phase transitions in a continuum model of the classical Heisenberg magnet: The ferromagnetic system. Phys. Rev. E 1994, 49, 5169−5178. (6) Lomba, E.; Weis, J.-J.; Tejero, C. F. Liquid-solid transition of the ferromagnetic Heisenberg fluid: Simulation, density functional, and perturbation theories. Phys. Rev. E 1998, 58, 3426−3435. (7) Lado, F.; Lomba, E. Heisenberg Spin Fluid in an External Magnetic Field. Phys. Rev. Lett. 1998, 80, 3535−3538. (8) Lado, F.; Lomba, E.; Weis, J.-J. Integral equation and simulation studies of the Heisenberg spin fluid in an external magnetic field. Phys. Rev. E 1998, 58, 3478−3489. (9) Tavares, J. M.; Telo da Gama, M. M.; Teixeira, P. I. C.; Weis, J. J.; Nijmeijer, M. J. P. Phase diagram and critical behavior of the ferromagnetic Heisenberg fluid from density-functional theory. Phys. Rev. E 1995, 52, 1915−1929. (10) Oukouiss, A.; Baus, M. Phase diagrams of the classical Heisenberg fluid within the extended van der Waals approximation. Phys. Rev. E 1997, 55, 7242−7252. (11) Giura, S.; Márkus, B. G.; Klapp, S. H. L.; Schoen, M. Isotropicpolar phase transitions in an amphiphilic fluid: Density functional theory versus computer simulations. Phys. Rev. E 2013, 87, 012313. (12) Cattes, S. M.; Gubbins, K. E.; Schoen, M. Mean-field density functional theory of a nanoconfined classical, three-dimensional Heisenberg fluid. I. The role of molecular anchoring. J. Chem. Phys. 2016, 144, 194704. (13) Erdmann, T.; Kröger, M.; Hess, S. Phase behavior and structure of Janus fluids. Phys. Rev. E 2003, 67, 041209. (14) Schoen, M.; Giura, S.; Klapp, S. H. L. Phase behavior of an amphiphilic fluid. Phys. Rev. E 2014, 89, 012310. (15) Wei, D.; Patey, G. N. Orientational order in simple dipolar liquids: computer simulation of a ferroelectric nematic phase. Phys. Rev. Lett. 1992, 68, 2043−2045. (16) Groh, B.; Dietrich, S. Ferroelectric phase in Stockmayer fluids. Phys. Rev. E 1994, 50, 3814−3833. (17) Klapp, S. H. L.; Forstmann, F. Phase transitions in dipolar fluids: An integral equation study. J. Chem. Phys. 1997, 106, 9742−9761. (18) Rovigatti, L.; Russo, J.; Sciortino, F. Structural properties of the dipolar hard-sphere fluid at low temperatures and densities. Soft Matter 2012, 8, 6310−6319. (19) Weeber, R.; Klinkigt, M.; Kantorovich, S.; Holm, C. Microstructure and magnetic properties of magnetic fluids consisting of shifted dipole particles under the influence of an external magnetic field. J. Chem. Phys. 2013, 139, 214901. (20) Schmidle, H.; Jäger, S.; Hall, C. K.; Velev, O. D.; Klapp, S. H. L. Two-dimensional colloidal networks induced by a uni-axial external field. Soft Matter 2013, 9, 2518−2524. (21) Cattes, S. M.; Klapp, S. H. L.; Schoen, M. Condensation, demixing, and orientational ordering of magnetic colloidal suspensions. Phys. Rev. E 2015, 91, 052127. (22) Bramley, A. The Study of Binary Mixtures. Part I. The Densities and Viscosities of Mixtures containing Phenol. J. Chem. Soc., Trans. 1916, 109, 10−45. (23) Rowlinson, J. S.; Swinton, F. L. Liquids and liquid mixtures: Butterworths monographs in chemistry; Butterworth-Heinemann: London, 1982. (24) Loudet, J.-C.; Barois, P.; Poulin, P. Colloidal ordering from phase separation in a liquid-crystalline continuous phase. Nature 2000, 407, 611−613. (25) Staff, R. H.; Rupper, P.; Lieberwirth, I.; Landfester, K.; Crespy, D. Phase behavior of binary mixtures of block copolymers and a nonsolvent in miniemulsion droplets as single and double nanoconfinement. Soft Matter 2011, 7, 10219−10226. (26) Raj, K.; Moskowitz, B.; Casciari, R. Advances in ferrofluid technology. J. Magn. Magn. Mater. 1995, 149, 174−180. (27) Holm, C.; Weis, J.-J. The structure of ferrofluids: A status report. Curr. Opin. Colloid Interface Sci. 2005, 10, 133−140. (28) Dijkstra, M.; Frenkel, D. Evidence for entropy-driven demixing in hard-core fluids. Phys. Rev. Lett. 1994, 72, 298−300.
(52b)
where A is defined as [cf. eq 8] 6 4εεH ⎛ σ22 ⎞ A(r12) ≡ − ⎜ ⎟ 3 ⎝ r12 ⎠
(53)
Last but not least, we introduce (ab)
̃ (r12) = 2πf (ab) (r12) f000 000
1
(22)
̃ (r12) = f000
4 π
a = b = 1, a ≠ b
(22) f 000 (r12)
(54a)
(54b)
for notational convenience. This permits us to define u0(ab) = −
∫σ
∞
(ab)
2 ̃ dr12 r12 f000 (r12),
a , b = 1, 2
ab
u1 =
1 π 33/2
∫σ
∞
2 dr12 r12 f110 (r12)
22
(55a)
(55b)
where we dropped the superscript (22) on both u1 and f110 to limit the notational burden. We also note in passing that the factor 1/4 in eq 54b is a consequence of the normalization of α̅ [see eq 17] which needs to be satisfied even in the isotropic phase. Because in the isotropic phase elements of the set {αl}l≥1 1 vanish, we can conclude that α0 = 2 .
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Stefanie M. Wandrei: 0000-0003-3300-4342 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work is dedicated to Professor Keith E. Gubbins (North Carolina State University) on the occasion of his 80th birthday. One of us (S.M.W.) wishes to express her gratitude for Professor Gubbins’ hospitality during her 6 month research visit to his laboratory. We acknowledge financial support from Deutsche Forschungsgemeinschaf t through funding for the International Graduate Research Training Group 1524. In addition, funding was provided from Deutscher Akademischer Austauschdienst through its Research Internship in Science and Engineering (RISE) program. The grant from RISE made it possible for D.G.M. to spend her research internship as an undergraduate student in Berlin during the summer of 2016.
■
REFERENCES
(1) Omelyan, I.; Mryglod, I.; Folk, R.; Fenz, W. Ising fluids in an external magnetic field: An integral equation approach. Phys. Rev. E 2004, 69, 061506. (2) Omelyan, I.; Fenz, W.; Mryglod, I.; Folk, R. X Y Spin Fluid in an External Magnetic Field. Phys. Rev. Lett. 2005, 94, 045701. (3) Omelyan, I.; Folk, R.; Kovalenko, A.; Fenz, W.; Mryglod, I. Liquid-vapor interfaces in X Y-spin fluids: An inhomogeneous anisotropic integral-equation approach. Phys. Rev. E 2009, 79, 011123. (4) Hemmer, P.; Imbro, D. Ferromagnetic fluids. Phys. Rev. A 1977, 16, 380−386. 11375
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
Article
Langmuir (29) Dijkstra, M.; van Roij, R.; Evans, R. Phase behavior and structure of binary hard-sphere mixtures. Phys. Rev. Lett. 1998, 81, 2268−2271. (30) Fenz, W.; Folk, R. Phase behavior of Ising mixtures. Phys. Rev. E 2005, 71, 046104. (31) Range, G. M.; Klapp, S. H. L. Density functional study of the phase behavior of asymmetric binary dipolar mixtures. Phys. Rev. E 2004, 69, 031201. (32) Range, G. M.; Klapp, S. H. L. Phase behavior of bidisperse ferrocolloids. Phys. Rev. E 2004, 70, 061407. (33) Range, G. M.; Klapp, S. H. L. Density-functional study of model bidisperse ferrocolloids in an external magnetic field. J. Chem. Phys. 2005, 122, 224902. (34) Szalai, I.; Dietrich, S. Global phase diagrams of binary dipolar fluid mixtures. Mol. Phys. 2005, 103, 2873−2895. (35) Díaz-Herrera, E.; Forstmann, F. The density and polarization of an ion-dipole-electrolyte near a charged wall. J. Chem. Phys. 1995, 102, 9005−9017. (36) Brunet, C.; Malherbe, J.; Amokrane, S. Demixing and fieldinduced population inversion in a mixture of neutral and dipolar-hard spheres confined in a slit pore. Mol. Phys. 2012, 110, 1161−1169. (37) Almarza, N. G.; Lomba, E.; Martín, C.; Gallardo, A. Demixing in binary mixtures of apolar and dipolar hard spheres. J. Chem. Phys. 2008, 129, 234504. (38) Hertlein, C.; Helden, L.; Gambassi, A.; Dietrich, S.; Bechinger, C. Direct measurement of critical Casimir forces. Nature 2008, 451, 172−175. (39) Verschaffelt, J. E. On the critical isothermal line and the densities of saturated vapour and liquid in isopentane and carbon dioxide. Proc. Konink. Nederl. Akad. Weten. 1899, 2, 588−592. (40) Vega, L.; de Miguel, E.; Rull, L. F.; Jackson, G.; McLure, I. A. Phase equilibria and critical behavior of square-well fluids of variable width by Gibbs ensemble Monte Carlo simulation. J. Chem. Phys. 1992, 96, 2296−2305. (41) Yin, T.; Lei, Y.; Huang, M.; Chen, Z.; Mao, C.; An, X.; Shen, W. Critical behavior of binary mixture of {xC 6 H 5 CN + (1−x) CH3(CH2)12CH3}: Measurements of coexistence curves, turbidity, and heat capacity. J. Chem. Thermodyn. 2011, 43, 656−663. (42) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids: Fundamentals; Claredon Press: Oxford, 1984; Vol. 1. (43) Schoen, M.; Haslam, A. J.; Jackson, G. Langmuir 2017, 10.1021/ acs.langmuir.7b01849. (44) Evans, R. The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids. Adv. Phys. 1979, 28, 143−200. (45) Hansen, J.-P.; McDonald, I. R. Theory of simple liquids, 4th ed.; Academic Press: New York, 2013. (46) Mansoori, G. A.; Carnahan, N. F.; Starling, K. E.; Leland, T. W., Jr. Equilibrium Thermodynamic Properties of the Mixture of Hard Spheres. J. Chem. Phys. 1971, 54, 1523−1525. (47) Boublík, T. Hard-Sphere Equation of State. J. Chem. Phys. 1970, 53, 471−472. (48) Salacuse, J. J.; Stell, G. Polydisperse systems: Statistical thermodynamics, with applications to several models including hard and permeable spheres. J. Chem. Phys. 1982, 77, 3714−3725. (49) Carnahan, N. F.; Starling, K. E. Equation of State for Nonattracting Rigid Spheres. J. Chem. Phys. 1969, 51, 635−636. (50) Widom, B. Theory of phase equilibrium. J. Phys. Chem. 1996, 100, 13190−13199. (51) Stanley, H. E. Introduction to phase transitions and critical phenomena; Oxford Science Publications: Oxford, 1971; Chapter 10. (52) Frodl, P.; Dietrich, S. Bulk and interfacial properties of polar and molecular fluids. Phys. Rev. A 1992, 45, 7330−7354.
11376
DOI: 10.1021/acs.langmuir.7b01952 Langmuir 2017, 33, 11366−11376
