In this work, we employ classical density functional theory (DFT) to investigate for the first time equilibrium properties of a Heisenberg fluid confined to nanoscopic slit pores of variable width. Within DFT pair correlations are treated at modified mean-field level. We consider three types of walls: hard ones, where the fluid-wall potential becomes infinite upon molecular contact but vanishes otherwise, and hard walls with superimposed short-range attraction with and without explicit orientation dependence. To model the distance dependence of the attractions, we employ a Yukawa potential. The orientation dependence is realized through anchoring of molecules at the substrates, i.e., an energetic discrimination of specific molecular orientations. If the walls are hard or attractive without specific anchoring, the results are “quasi-bulk”-like in that they can be linked to a confinement-induced reduction of the bulk mean field. In these cases, the precise nature of the walls is completely irrelevant at coexistence. Only for specific anchoring nontrivial features arise, because then the fluid-wall interaction potential affects the orientation distribution function in a nontrivial way and thus appears explicitly in the Euler-Lagrange equations to be solved for minima of the grand potential of coexisting phases.

If a fluid composed of spherically symmetric molecules has some superimposed vectorial degrees of freedom, ordered liquid-like phases can be formed. These systems are attracting attention from the scientific community for quite some time and continue to do so.1,2 Perhaps the simplest model system in this regard is the three-dimensional, classical Heisenberg fluid3–8 which consists of molecules with isotropic cores and a superimposed classical “spin-spin” interaction such that ordered polar fluid phases can be formed. Therefore, the Heisenberg fluid is a model fluid for so-called ferromagnetic (or antiferromagnetic) fluids.

These fluids are still receiving a lot of interest because a broad range of applications exists.9 These applications include lubrication,10 dynamic sealing,11 mineral processing and water treatment,12 and, more recently, radio-frequency-induced hyperthermia in cancer treatment13 to name only a few examples.

Closely related to the Heisenberg fluid are dipolar fluids composed of (hard or soft) spheres. In these latter model systems, the classical spin is replaced by a point dipole. Whereas the orientation dependence of the dipole-dipole interaction differs from that of the spin-spin interaction, the resulting phase diagrams are topologically rather similar. This is evident by comparing the phase diagrams obtained for the bulk Stockmayer fluid by Groh and Dietrich (see Fig. 15 of Ref. 14) with those presented by Tavares et al. (see Figs. 1–3 of Ref. 4) and later by Schoen et al. (see Figs. 1(a)–1(c) of Ref. 8) for various versions of the Heisenberg fluid.

This is by far not a triviality because it is the long-range character of the dipole-dipole interaction in combination with the orientation dependence that causes the formation of ordered phases whereas short-range spin-spin interactions suffice to cause the formation of polar phases in Heisenberg fluids.

Because bulk properties of Heisenberg fluids3–8 are well understood, it seems interesting to study the impact of confinement. In general, confinement refers to a situation where a bulk phase is exposed to some sort of external field.15 

A special realization of such a field, to which we will restrict the discussion in this paper exclusively, is the solid walls of a nanoscopic slit pore. In this case the interesting physics is caused by an interplay or rather competition between two length scales, namely, the bulk correlation length and the distance between the confining walls of the slit pore which prevents the correlation length from exceeding certain limits in the direction of the confining walls.

In the case of confined dipolar fluids, a rather extensive body of literature already exists which cannot be reviewed exhaustively in a research paper such as the present one. Instead, only a few key publications will be reviewed to give the reader some flavor of the kind of problems that have been studied already.

By means of molecular dynamics (MD) simulations, Lee et al. investigated a confined Stockmayer film in which the dipoles are oriented in a direction parallel to the confining walls if no electric field is present.16 By the same technique, Jordanovic and Klapp studied the competition between the direction of external, homogeneous magnetic fields, and confinement by solid walls.17 Klapp and Schoen focused on the spontaneous ordering of confined dipolar films18 and later on the nature of ferroelectric states in these systems19 by means of Monte Carlo simulations. Ballenegger and Hansen used a combination of linear response theory and MD simulations to study the dielectric response of a confined dipolar fluid.20 Last but not least, Szalai and Dietrich21 focus on phase transitions and orientational effects in confined dipolar fluids based upon density functional theory (DFT) at the modified mean-field level, a version of DFT which is also employed here.

Earlier work on confined dipolar fluids that is particularly relevant to the present one is the study by Gramzow and Klapp.22 First, this is also a DFT study in which pair correlations are treated at the modified mean-field level. Second, these earlier authors are using the same approximation to the singlet density by ignoring packing effects and treating only the orientation distribution function as a local quantity. The difference between our approach and the one by Gramzow and Klapp22 is, however, that in our DFT, we solve the Euler-Lagrange equations under the explicit assumption of phase coexistence, whereas thermodynamic equilibrium states off coexistence are also considered in Ref. 22.

Our motivation to carry out the present study is threefold. First, compared with the rather impressive list of investigations of equilibrium properties of confined dipolar fluids, there does not seem to exist a single study of confined Heisenberg fluids to the best of our knowledge. This is surprising in view of the fact that the treatment of dipolar fluids is much more demanding on account of the long-range nature of the dipole-dipole interaction potential.14,23,24

Second, we will demonstrate here that in previous and related works on confined dipolar fluids, the presence of confining walls causes only a “quasi-bulk” effect at phase coexistence, that is, the walls only impose a reduction of the bulk mean field whereas details of the fluid-wall interaction potential turn out to be irrelevant.

Third, this latter situation only changes if the fluid-wall potential becomes explicitly orientation-dependent. This, in turn, gives rise to rather complex orientational phenomena that have not been explored before.

The remainder of this manuscript is organized as follows. In Sec. II, we introduce our model system. Elements of mean-field DFT are summarized in Sec. III whereas Sec. IV is devoted to a presentation of our results. They are then summarized and discussed in Sec. V. Details of the calculation of the bulk mean field can be found in the  Appendix.

We consider a fluid of N particles interacting in a pairwise additive fashion via the orientation-dependent intermolecular potential function

φr12,ω1,ω2=φisor12+φanisr12,ω1,ω2
(2.1)

which we decompose formally into an isotropic contribution φiso and an anisotropic one represented by φanis. In Eq. (2.1), r12 = r1r2 is the distance vector between the centers of mass of a pair of molecules located at r1 and r2, respectively, r12=r12, and ω1 and ω2 are sets of Euler angles specifying the orientation of both molecules in a space-fixed frame of reference (see also below).

We assume our molecules to have short-range interactions such that

φisor12=4εσr1212σr126=φrepr12+φattr12
(2.2)

is given by the well-known Lennard-Jones potential, where ε is the depth of the attractive well and σ is the van der Waals diameter. In addition, we consider a hard-sphere reference fluid where

φhsr12=,r12<σ0,r12σ
(2.3)

describes the interaction between a pair of hard spheres of diameter σ. The concept of a reference fluid is convenient from the standpoint of DFT because it will allow us to compute the contribution of intermolecular interactions to the free energy from a perturbational approach (see below).

To account for anisotropic interactions, we notice that φanis can be expanded according to25 

φanisr12,ω1,ω2=l1,l2,lφl1l2lr12Φl1l2lω1,ω2,ω,
(2.4)

where Φl1l2l is a rotational invariant, l1, l2, and l are nonnegative integers, ω is a set of Euler angles describing the orientation of  r̂12=r12/r12 in a space-fixed frame of reference, and φl1l2l are expansion coefficients. Here and below we indicate unit vectors by the caret. Focusing exclusively on dispersion interactions φl1l2l = φatt25 irrespective of l1, l2, and l.

To describe the interaction between a pair of classical, three-dimensional (Heisenberg) “spins,” we take

φanisr12,ω1,ω2=εHφattr12cosγ,
(2.5)

where γ is the angle between the spins with orientations specified by ûω1 and ûω2. In Eq. (2.5), the (dimensionless) coupling constant εH > 0 corresponds to the ferromagnetic case.

To rationalize the dependence of φanis on γ, we introduce

Φl1l2lω1,ω2,ω=m1,m2,mCl1l2l;m1m2m×Yl1m1ω1Yl2m2ω2Ylm*ω,
(2.6)

where C is a Clebsch-Gordan coefficient, Ylm is a spherical harmonic, the asterisk denotes its complex conjugate, and corresponding pairs of integers are related through m′ ∈ − l′, …, l′. From this general definition, it then follows that26 

Φ110=14π13m=11Y1m*ω1Y1mω2=34π3/2P1cosγ=34π3/2cosγ,
(2.7)

where P1 is the first Legendre polynomial, and we temporarily dropped the arguments of Φ110 for notational convenience. The relation between Φ110 and P1 follows directly from the addition theorem for spherical harmonics [see Eq. (A.33) of Ref. 25]. From Eqs. (2.5) and (2.7), we then have

φanisr12,ω1,ω2=4π3/23εHφattr12Φ110ω1,ω2,ω.
(2.8)

Two comments apply at this stage. First, it is clear that our molecules have uniaxial symmetry such that ωi=θi,ϕi (i = 1, 2) is a subset of two out of three Euler angles. Second, upon inverting one of the two spins by replacing ωiωi=ωi=πθi,π+ϕi, φanis → −φanis as this replaces γ by γ′ = π + γ in Eq. (2.7). Moreover, φanis is invariant upon interchanging the center-of-mass positions of molecules i and j, that is, if  r̂ijr̂ij=r̂ij. Hence, our model fluid consists of polar molecules and is therefore capable of forming ordered polar phases.

The model fluid introduced in Sec. II A is placed between two planar, structureless walls separated from each other by a distance sz along the z-axis. To describe the interaction of a fluid molecule with the walls, we introduce

φwkz=φhwzζκexpκz±szσ/2z±szσ/2gθϕkω,
(2.9)

where we assume that the lower wall (k = 1, +) to be located at −sz/2 whereas the upper one (k = 2, −) is at sz/2. In Eq. (2.9), the parameter

ζκ=σ2expσκ/2εw
(2.10)

is introduced to guarantee that φw = − εw at z=sz/2 irrespective of the inverse Debye screening length κ and

φhwz=,z<szσ/20,zszσ/2
(2.11)

is a hard-wall background potential. We have chosen the Yukawa potential because it allows us to vary the range of the attraction continuously through a variation of κ and its strength by varying εw.

In Eq. (2.9), 0gθϕkω1 is an orientation-dependent, so-called anchoring function. It serves as a mathematical device to discriminate energetically between desirable and undesirable orientations of the molecule. Experimentally, anchoring can be realized in a number of different ways by either manipulating the walls mechanically (by, for example, rubbing) or by functionalizing them chemically.27,28

In this work, we adopt

gθϕkω=121+ûωêθϕk,
(2.12)

where

êθϕk=sinθkcosϕksinθksinϕkcosθk=cosϕksinϕk0
(2.13)

and θk and ϕk are polar and azimuthal angles, respectively. These angles determine the alignment of molecules at wall k with respect to the wall plane. Throughout this work, we focus exclusively on anchoring scenarios such that the energetically favorable alignment of molecules is in some direction in the xy plane characterized by θ1=θ2=π2. However, in general, ϕ1 may be different from ϕ2 which allows us to implement hybrid anchoring scenarios where homogeneous anchoring (i.e., ϕ1=ϕ2) is included as a special case. A sketch of our model system is presented in Fig. 1.

FIG. 1.

Sketch of a Heisenberg fluid confined to a slit pore with explicit hybrid anchoring (indicated by the big red and blue arrows pointing up and down, respectively). Arrows attached to the spherically symmetric particles represent their classical spins; (a) isotropic phase, (b) polar phase.

FIG. 1.

Sketch of a Heisenberg fluid confined to a slit pore with explicit hybrid anchoring (indicated by the big red and blue arrows pointing up and down, respectively). Arrows attached to the spherically symmetric particles represent their classical spins; (a) isotropic phase, (b) polar phase.

Close modal

For the model introduced in Sec. II, we can associate thermodynamic states at equilibrium with minima of the grand potential29 

Ωρr,ω=Fρr,ωμdrdωρr,ω,
(3.1)

where F is the intrinsic Helmholtz free energy and μ is the chemical potential. Both Ω and F are functionals of the orientation-dependent local singlet density ρr,ω. In Sec. III A, we will briefly sketch how these thermodynamic equilibrium states can be identified.

As a first approximation, we immediately introduce

ρr,ω=ραr,ω=ραz,ω,
(3.2)

where ρ is the mean number density and α is the local orientation distribution function normalized according to

dωαz,ω=1.
(3.3)

Because of confinement in the z-direction, properties of the fluid are translationally invariant in each of the x-y planes stacked along the z-axis. Therefore, α turns out to depend only on z and ω. In the isotropic phase, α=1/4π for molecules of uniaxial symmetry.

Strictly speaking, however, the approximation introduced in Eq. (3.2) ignores packing effects that arise on account of the fluid-wall interaction. In reality, these packing effects manifest themselves as oscillations in the local density. The periodicity of the oscillations is closely related to σ. They are damped as one moves away from the walls depending on the range of the fluid-wall potential. Eventually, at sufficiently large distances from the walls, the local density will approach ρ. Hence, Eq. (3.2) is only approximately valid for sufficiently large sz and in regions located at distances from the walls that are large compared with 1/κ.

On the one hand, as was shown earlier by Gramzow and Klapp in their study of confinement effects in dipolar fluids,22 Eq. (3.2) allows one to capture confinement-induced orientational phenomena sufficiently realistically. On the other hand, the approximation introduced in Eq. (3.2) is advantageous in several ways. This is because it allows us to decompose F according to

F=For+Fw+Fhs+ΔFex,
(3.4)

where all members of the equation are functions of ρ and functionals only of α. The first term on the right side of Eq. (3.4) is an entropic (ideal-gas) contribution that can be cast as

βForV=ρszsz/2sz/2dzdωαz,ωln4παz,ω,
(3.5)

where β = 1/kBT (kB is the Boltzmann’s constant and T is the temperature), and V denotes volume. The expression in Eq. (3.5) reflects the loss in orientational entropy when a polar fluid phase forms. The factor 4π in the argument of the logarithm is introduced to make sure that For vanishes in the isotropic phase (see above).

The next term in Eq. (3.4) represents the contribution of the external field represented by the confining walls. It can be cast as

βFwV=ρβszk=12sz/2sz/2dzdωαz,ωφwkz,ω,
(3.6)

where φwk given in Eqs. (2.9) and (3.2) has also been used.

Because of the approximation of nonlocality of the density of the confined fluid [see Eq. (3.2)], we can describe the hard-sphere reference fluid by a suitable equation of state. Here we adopt the one proposed by Carnahan and Starling.30 Including also the kinetic (i.e., ideal gas) contribution, we may thus write

βFhsV=βfhs=ρlnρΛ5m/I1+ρ4η3η21η2
(3.7)

for the free-energy density fhs of the hard-sphere reference system, where Λ is the thermal de Broglie wavelength, m the molecular mass, I the moment of inertia, and η = πρσ3/6 is the hard-sphere packing fraction. The exponent of 5 reflects the number of degrees of freedom (three translational and two rotational) of our uniaxial molecules.

Last but not least, the contribution of φ to F can be written in a perturbative spirit as31,32

ΔFex=1201dλdr1dr2dω1dω2ρr1,ω1ρr2,ω2×gr1,r2,ω1,ω2;λλ[φhsr12+λφr12,ω1,ω2],
(3.8)

where the dimensionless parameter λ allows us to specify a one-dimensional path of integration that takes us from the hard-sphere reference system to the system of interest. In Eq. (3.8), g is the orientation-dependent pair correlation function. In its most general form, g depends on r1 and r2 separately, because of the presence of the walls. However, because of the approximation introduced in Eq. (3.2) (namely, homogeneity of the density of the confined fluid), g depends only on r12.

We approximate g at the so-called modified mean-field level4,8,14,23,24,26,33,34 as

gr12,ω1,ω2;λ=0,r12<σexpβφhs+λφ,r12σ,
(3.9)

where we temporarily dropped the arguments of both φhs and φ for notational convenience. This expression can be viewed as the leading term in a density expansion of g and therefore becomes exact in the limit of vanishing ρ. Using Eq. (3.9) the integration over λ in Eq. (3.8) can be performed analytically.

Next, we remind ourselves that the confined fluid has cylindrical symmetry [as reflected, for example, by Eq. (3.2)]. To proceed it is therefore advisable to introduce cylindrical coordinates ri=ai,zi, where ai=ai,γi is the projection of ri onto the xy plane and γi is the polar angle. Thus, it is useful to transform variables such that a1,γ1a1,γ1=a1,γ1 and a2,γ2a12,γ12.

One can then perform the trivial integration over da1 and dγ1 in Eq. (3.8) analytically and finally obtain

βΔFexV=ρ22szsz/2sz/2dz1dz2a12da12Θr12σ×dω1dω202πdϕαz1,ω1αz2,ω2×fr12,ω1,ω2,
(3.10)

where Θ is the Heaviside function and

fr12,ω1,ω2=expβφr12,ω1,ω21
(3.11)

is the orientation-dependent Mayer f-function arising from the analytic integration over dλ in Eq. (3.8) where Eqs. (3.2) and (3.9) have also been employed. In addition, we used the fact that γ12 = ϕ, where ϕ is the azimuth describing the orientation of  r̂12 in spherical coordinates.

The integrations over dω1, dω2, and dϕ in Eq. (3.10) can be carried out analytically. To achieve this we first expand the orientation distribution functions in Eq. (3.10) according to

αzi,ωi=li,miαlimiziYlimiωi,i=1,2.
(3.12)

This expansion introduces a fourfold summation over l1,m1 and l2,m2 into Eq. (3.10) where corresponding pairs of these integers l,m are related as already explained below Eq. (2.6).

Next, we expand the Mayer f-function [Eq. (3.10)] in terms of rotational invariants according to

fr12,ω1,ω2=l1,l2,lfl1l2lr12Φl1l2lω1,ω2,ω.
(3.13)

The expansion coefficients in Eq. (3.13) can be obtained from

fL1L2Lr12=4π2L+1dω1dω2dωfr12,ω1,ω2×ΦL1L2L*ω1,ω2,ω
(3.14)

because Φl1l2l is a complete set of orthogonal basis functions, that is,

dω1dω2dωΦl1l2lΦL1L2L*=2L+14πδl1L1δl2L2δlL,
(3.15)

where the δ’s are Kronecker symbols.

Replacing now in Eq. (3.10), f via Eq. (3.13) and using Eq. (2.6), it is clear that this introduces six more summations over integer pairs L1,M1, L2,M2, and L,M. Thus, we now have to deal with a tenfold summation over such integer pairs. However, this seemingly too cumbersome problem can be simplified substantially.

Focusing only on orientation-dependent terms in Eq. (3.10), we notice that

L1,L2,LM1,M2,MCL1L2L;M1M2M02πdϕYLM*ωdω1Yl1m1YL1M1dω2Yl2m2YL2M2=2πL,LMCLLL;MM¯0YL0θdω1Yl1m1YLMdω2Yl2m2YLM¯=2πL,LMCLLL;MM¯0YL0θδl1Lδm1Mδl2Lδm2M,
(3.16)

where θ is the polar angle specifying the orientation of r̂12 in spherical coordinates and M¯=M.

To rationalize Eq. (3.16), we first realize that

02πdϕYLM*ω02πdϕexpiMπ=2πδM0
(3.17)

so that only YL0θ remains. Second, the selection rule M1 + M2 = M for nonzero Clebsch-Gordan coefficients has been used from which M1 = − M2 = M follows. However, this implies also L1 = L2 = L′. Third, we utilize [see Eq. (A.3) of Ref. 25]

YLM¯ω=1MYLM*ω
(3.18)

which allows us to carry out the integrations over dω1 and dω2 analytically and from which the last line of Eq. (3.16) emerges [see Eq. (A.39) of Ref. 25].

Using now Eq. (3.16) together with Eqs. (3.12) and (3.13), Eq. (3.10) may be recast as

βΔFexV=ρ2szL,LMsz/2sz/2dz1dz2αLM¯z1αLMz2×uLLMz12,
(3.19)

where the coefficients

uLLMz12=πCLLL;MM¯0×a12da12YL0θa12fLLLr12×Θr12σ
(3.20)

can only depend on z12=z1z2 because both r12=a122+z122 and cosθ = a12/r12 depend on z12.

Finally, we seek to eliminate the sum on L from the expressions given above. To that end we expand the Mayer f-function in Eq. (3.11) as a power series in terms of βφanis where we assume this quantity to be less than one. From Eq. (2.5), we realize that φanisφatt. We notice also that φatt is smallest at r12 = σ. The temperature range will be such that kBT/ε is of the order of one. Hence, we anticipate βφanis to be of the order of 10−2 to 10−1 given that εH = 0.06 employed here throughout. This suggests that the expansion of the Mayer f-function might actually be converging rather rapidly. Inserting the expansion now into Eq. (3.14), one can show that

uLMz12=π4π1L+M2L+1a12da12fLL0r12×Θr12σ,
(3.21)

where the countably infinite set uLM constitutes the bulk mean field. Details of this derivation can be found in the  Appendix. Thus, we finally rewrite Eq. (3.10) as

βΔFexV=ρ2szL,Msz/2sz/2dz1dz2αLM¯z1αLMz2uLMz12.
(3.22)

Collecting the individual contributions to the free energy functional given in Eqs. (3.5)(3.7) and (3.22), the grand potential functional Ω can be cast as

βΩV=βfhs+ρszsz/2sz/2dzdωαz,ωln4παz,ω+ρβszk=12sz/2sz/2dzdωαz,ωφwkz,ω+ρ2szL,Msz/2sz/2dz1dz2αLM¯z1αLMz2×uLMz12βρμ
(3.23)

indicating that Ω depends on α through the second, third, and fourth terms on the right side. At fixed T, thermodynamic equilibrium states correspond to minima of Ω with respect to ρ and α. Hence, we are seeking simultaneous solutions of the equations

βVΩρ=0,
(3.24a)
βVδΩδαz,ω=λz;T,ρ,
(3.24b)

where the operator δ indicates a functional derivative and λ is a Lagrangian multiplier introduced to ensure that Eq. (3.3) is always satisfied locally (i.e., at every z).

It is then a simple matter to verify that from Eqs. (3.23) and (3.24a), we have

0=βμhs+1szsz/2sz/2dzdωαz,ωln4παz,ω+βszk=12sz/2sz/2dzdωαz,ωφwkz,ω+2ρszL,Msz/2sz/2dz1dz2αLM¯z1αLMz2×uLMz12βμ,
(3.25)

where

βμhs=lnρΛ5m/I+8η9η2+3η31η3
(3.26)

is the chemical potential of the hard-sphere reference fluid obtained from Eq. (3.7). Solving Eq. (3.25) for βμ and replacing this term in Eq. (3.23) by the resulting expression immediately give

βΩV=βfhsβρμhsρ2szL,Msz/2sz/2dz1dz2αLM¯z1αLMz2uLMz12.
(3.27)

At this point a couple of comments seem to be useful. First, one realizes from Eq. (3.27) that at thermodynamic equilibrium, Ω is independent of the a priori unknown chemical potential μ which is a rather gratifying feature of the above approach.

Second, a term ρlnρΛ5m/I arises in the expressions for βfhs [see Eq. (3.7)] and βρμhs [see Eq. (3.26)] which consequently disappears when the latter is subtracted from the former. Hence, at thermodynamic equilibrium, Ω is independent of any material-specific constants as it must.

Third, the last term in Eq. (3.27) describes a contribution to the grand potential arising from a coupling between local order parameters αLM. The strength of this coupling is controlled by the bulk mean field uLM. Interestingly, uLM is solely determined by intermolecular interactions and is independent of the interaction between molecules and the confining walls as an inspection of Eqs. (3.11), (3.13), and (3.21) reveals. Consequently, uLM is a bulk property. However, if both molecules are sufficiently close to the walls, the bulk mean field will be reduced somewhat on account of the Heaviside function appearing in the integrand in Eq. (3.21). We shall refer to features of the confined fluid that can solely be ascribed to the wall-induced reduction of the bulk mean field as “quasi-bulk-like.”

Using now Eq. (3.23) together with the second equilibrium condition Eq. (3.24b) and Eqs. (3.3) and (3.12), one finds that

szλz;T,ρρρ=ρ2πsz/2sz/2dzu00zzlnΨzω,
(3.28)

where the angle-averaged orientation distribution function is given by

Ψzω=14πdωΨz,ω
(3.29)

and

Ψz,ω=exp[2ρL=1Msz/2sz/2dzuLMzz×αLMzYLM¯*ω]expβk=12φwkz,ω.
(3.30)

Notice that the first term on the right side of Eq. (3.28) arises because we have treated the isotropic contribution (L = M = 0) in Eq. (3.30) separately.

Having derived in Eq. (3.27) an expression for the grand-potential density at thermodynamic equilibrium, we are now in a position to focus on coexisting equilibrium phases. Let ′ and ″ be two equilibrium phases with densities ρ′ < ρ″ coexisting at temperature T. Because of the presence of the walls, both phases are characterized by sets of local order parameters αLM and αLM in general.

Because ′ and ″ are at thermodynamic equilibrium, we have from Eq. (3.24a)

βVΩρρ,αLM=βVΩρρ,αLM=0.
(3.31)

In addition, at coexistence ′ and ″ must satisfy

βΩV=βΩV.
(3.32)

Combining this second coexistence condition with Eq. (3.27), it is evident that phase coexistence is associated with zeros of the function

s1=βfhsρβρμhsρβfhsρ+βρμhsρρ2ρ24π0szdz12u00z12+βΔFexVβΔFexV,
(3.33)

where from Eq. (3.22)

βΔFex,V=ρ,2szL=1Msz/2sz/2dz1dz2αLM¯,z1αLM,z2×uLMz12.
(3.34)

Notice that the expression on the second line of Eq. (3.33) arises by considering in Eq. (3.22) the term for L = M = 0 separately because it is always present regardless of whether we have a disordered or an ordered phase. Here we used the fact that because of Eqs. (3.3) and (3.12), α00=1/4π changed variables according to z1z1=z1 and z2z12 and performed the one trivial integration over dz1 analytically.

To implement the first coexistence condition Eq. (3.31), we take α from Eq. (3.24b) and insert it into Eq. (3.25). Replacing in the resulting expression the Lagrangian multiplier λ via Eq. (3.28) eventually gives

βμ=βμhs+ρ2π0szdz12u00z121szsz/2sz/2dzlnΨzω.
(3.35)

Thus, at coexistence (μ′ = μ″), we are seeking zeros of the function

s2=βμhsρβμhsρ1szsz/2sz/2dzlnΨzωΨzω+ρρ2π0szdz12u00z12
(3.36)

which is again independent of the unknown chemical potential and a term proportional to lnΛ5m/I [see Eq. (3.26)].

Finally, we realize that Eq. (3.30) involves L2L+1 unknown order parameters because of the relation between L and M. These order parameters can be determined from zeros of the functions

sl+2,z=αLMzdωΨ,z,ωYLM*ωdωΨ,z,ω.
(3.37)

To solve Eqs. (3.33) and (3.36), we employ a Newton-Raphson scheme described in Ref. 8 whereas Eq. (3.37) is solved by an algorithm proposed by Ng.35 

Before turning to a discussion of our results, we obviously need to truncate the infinite sums over L and M in Eqs. (3.34) or (3.30) at some sensible level. Indeed, we can limit ourselves to L = 1, M = 0, ± 1 because for these values of L and M, nonzero order parameters αLM would indicate a certain degree of polar order present in the confined Heisenberg fluid. We follow Gramzow and Klapp22 and introduce the polarization Pz=Pxzêx+Pyzêy+Pzzêz such that the total polarization Pz=Pz1. Because of the relation of Cartesian components of unit vectors to spherical harmonics [see Eq. (A.63) of Ref. 25], it is easy to verify that

Pxz=dωαz,ωsinθcosϕ=4π3121/2dωαz,ωY11ω+Y11*ω=8π3Reα11z,
(4.1a)
Pyz=dωαz,ωsinθsinϕ=8π3Imα11z,
(4.1b)
Pzz=dωαz,ωcosθ=4π3α10z,
(4.1c)

where “Re ” and “Im ” refer to real and imaginary parts, respectively.

Henceforth, we shall present all quantities in the customary dimensionless (i.e., “reduced”) units. For instance, length will be given in units of σ, density in units of σ−3, and temperature in units of ε/kB.

We begin by presenting results for the simplest possible situation, namely, that of hard substrates for which εw = 0 [see Eq. (2.9)]. Clearly, in this case, anchoring at the walls is irrelevant. From Eqs. (2.9) and (2.10), it is then clear that φwkz=φhwz. One then realizes from Eq. (3.30) that for hard walls, the somewhat simpler expression

Ψz,ω=exp[2ρL=1Msz/2sz/2dzuLMzz×αLMzYLM¯*ωΘzszσ/2]
(4.2)

is obtained.

For the relatively small Heisenberg coupling constant εH = 0.06, one anticipates a phase diagram of type I in the bulk according to the classification scheme of Schoen et al.8 Our rationale for this choice is that in a type I phase diagram all three phases, namely, gas (G), isotropic liquid (I), and a polar phase (P) are thermodynamically stable depending on T and ρ. At the same time, the topology of the phase diagram is sufficiently simple.

Plots in Fig. 2 support the notion of a phase diagram of type I in the bulk and in confinement. Below T = 1.042, one observes coexistence between a gas (G) and a polar phase (P) in the bulk; above that temperature, a G phase coexists with an isotropic liquid (I). At Tc ≃ 1.238 and ρc ≃ 0.250, a GI critical point is found in agreement with our earlier study.8 Hence, T ≃ 1.042 marks a critical end point at which a critical line would start that we omit for reasons of clarity (see, for example, Fig. 1(a) of Ref. 8).

FIG. 2.

(a) Phase diagrams in Tρ representation for a Heisenberg fluid confined to hard-wall slit pores of variable widths; (black square) bulk, (red circle) sz = 15.0, (green triangle) sz = 7.5, and (blue triangle down) sz = 5.0.

FIG. 2.

(a) Phase diagrams in Tρ representation for a Heisenberg fluid confined to hard-wall slit pores of variable widths; (black square) bulk, (red circle) sz = 15.0, (green triangle) sz = 7.5, and (blue triangle down) sz = 5.0.

Close modal

Confining the Heisenberg fluid to slit pores of variable pore widths sz does not affect the overall topology of the phase diagram. It merely serves to shift the critical temperature down monotonically with decreasing sz such that at the same T the densities of G and I or P phases become more alike because at coexistence the density of the G phase is somewhat less affected by the magnitude of sz compared with those of the I or P phases. This shift of the phase diagrams is to be expected because a fluid molecule is exposed to a reduced bulk mean field the closer the walls get. These attractive interactions affect more the high- than the low-density side of the phase diagram according to one’s intuition.

For four thermodynamic state points identified by symbols in the phase diagrams presented in Fig. 2 we present the corresponding profiles of the total polarization Pz in Fig. 3 for different pore widths. Because of the shift of the phase diagrams with sz we can only maintain ρ but have to adjust T to generate the curves presented in Fig. 3. We note in passing that plots similar to the ones shown in Figs. 2 and 3 can also be generated for thermostate points along an isotherm but varying densities.

FIG. 3.

Plots of the total polarization Pz across slit pores of variable widths; (black square) bulk, (red circle) sz = 15.0, (green triangle) sz = 7.5, and (blue triangle down) sz = 5.0. Thermodynamic state points T,ρ are chosen as marked in Fig. 2.

FIG. 3.

Plots of the total polarization Pz across slit pores of variable widths; (black square) bulk, (red circle) sz = 15.0, (green triangle) sz = 7.5, and (blue triangle down) sz = 5.0. Thermodynamic state points T,ρ are chosen as marked in Fig. 2.

Close modal

As anticipated, Pz0.57 does not depend on z for the bulk system. In a wide pore of sz = 15.0, one notices the effect of confinement: Pz declines strongly in the vicinity of the walls (and will eventually go to zero directly at the walls). However, this pore is still wide enough such that a bulk-like region exists in which Pz is independent of z. This is not so for smaller pores where the curves corresponding to sz = 7.5 and sz = 5.0 show that Pz depends on z across the entire pores and passing through maxima at z = 0.

The decline of Pz as zsz/2 can again be described to the quasi-bulk effect. Because the phase diagrams in Fig. 2 approach the bulk limit as sz → ∞, a similar shift is expected for Pz. Plots in Fig. 3 confirm this expectation.

Next we turn to a slightly more complicated case, namely, that of attractive walls with nonspecific anchoring which we realize by choosing εw > 0 and gθϕk=1 in Eq. (2.9). It is then clear that in Eq. (3.30), φwkz,ωφwkz. Hence, when performing integrations over dω in Eq. (3.37), expβk=12φwk [see Eq. (3.30)] cancels between denominator and numerator of the expression on the right side.

A similar cancellation occurs in Eq. (3.36) between denominator and numerator of the argument of the logarithmic term. Because at the same time Eq. (3.33) is independent of the fluid-wall attraction one anticipates that for fixed sz neither the strength of this attraction nor its range has any effect on the phase diagram of the confined Heisenberg fluid. Plots of phase diagrams and of the corresponding polarization profiles in Fig. 4 confirm this expectation.

FIG. 4.

(a) Phase diagrams in T-ρ representation for sz = 5.0 and εw = 0.0 (black solid line), εw = 1.0 (red square), and εw = 2.0 (green circle) [κ = 1.0, see Eq. (2.9)]. Also shown are results for the inverse Debye lengths κ = 0.2 (blue triangle) and κ = 2.0 (magenta triangle down) (εw = 2.0). (b) Total polarization Pz across a slit pore where symbols and line refer to the same values of εw and κ as in part (a). Data are shown for T = 0.872.

FIG. 4.

(a) Phase diagrams in T-ρ representation for sz = 5.0 and εw = 0.0 (black solid line), εw = 1.0 (red square), and εw = 2.0 (green circle) [κ = 1.0, see Eq. (2.9)]. Also shown are results for the inverse Debye lengths κ = 0.2 (blue triangle) and κ = 2.0 (magenta triangle down) (εw = 2.0). (b) Total polarization Pz across a slit pore where symbols and line refer to the same values of εw and κ as in part (a). Data are shown for T = 0.872.

Close modal

Before closing this section, we note in passing that different strengths of fluid-wall attraction have also been considered by Marguta et al. in their study of a Lebwohl-Lasher (lattice) model of a confined liquid crystal (see Sec. V).36 

1. Homogeneous anchoring

The above situation only changes if molecules are anchored at the walls in specific ways. Consider, for example, the case of homogeneous anchoring characterized by gθϕ1=gθϕ2=gθϕ. More specifically, we take

gθϕω=ûωêx=8π3ReY11ω,
(4.3)

where we used Eq. (A.63) of Ref. 25. Thus, molecules align themselves with the positive x-axis at both walls. Because the anchoring function causes φwk to be explicitly orientation-dependent, no cancellation of expβk=12φwk in Eq. (3.37) occurs and hence solutions of Eqs. (3.36) and (3.37) depend explicitly on the fluid-wall interaction.

This is illustrated by plots of phase diagrams in Fig. 5 where one notices that they are affected most significantly over a temperature range around the critical end point which, however, arises only in the case of confinement by hard walls (see Fig. 3) or for nonspecific anchoring [see Fig. 4(a)]. At the critical end point, the dT/dρ changes discontinuously.

FIG. 5.

Phase diagrams for a pore width sz = 5.0, κ = 0.2 and homogeneous anchoring of molecules along the x-axis at both walls; εw = 0.0 (black square), εw = 1.0 (red circle), and εw = 2.0 (green triangle).

FIG. 5.

Phase diagrams for a pore width sz = 5.0, κ = 0.2 and homogeneous anchoring of molecules along the x-axis at both walls; εw = 0.0 (black square), εw = 1.0 (red circle), and εw = 2.0 (green triangle).

Close modal

For εw > 0 and homogeneous anchoring, the critical end point is transformed into a shoulder which becomes increasingly weaker with increasing εw. Thus, for εw > 0, dT/dρ varies continuously with ρ. Compared with the temperature range around the critical end point, other parts of the phase diagram are far less affected by the strength of the fluid-wall interaction as plots in Fig. 5 reveal.

The impact of the fluid-wall attraction with specific anchoring is very much akin to that of an external, spatially homogeneous magnetic field H.37 There it was found that even an infinitesimally small magnetic field destroys a similar critical end point and with it the associated critical line for the same type I phase diagram.

The similarity between the impact of H and that of a wall with specific anchoring becomes clear if one realizes that the external magnetic field is proportional to the angle γ between H and the molecular orientation û. Hence, applying Eq. (3.24b) to the magnetic-field contribution to Ω causes H to appear explicitly in the orientation distribution function [see, for instance, Eq. (3.22) of Ref. 37] in very much the same way as an anchoring-dependent fluid-wall interaction does in the present case.

Anchoring at the walls renders the confined Heisenberg fluid strongly anisotropic and orientationally inhomogeneous. This is illustrated by plots in Fig. 6 which show that Px varies nonmonotonically across the pore whereas Py = Pz ≃ 0. In other words, gθϕ in Eq. (4.3) completely dominates the direction of P if the pore is small enough.

FIG. 6.

Plots of components Pαz (α = x, y, or z) of the polarization vector Pz as functions of position z for a slit pore of width sz = 5.0, εw = 2.0, and κ = 0.2. Data are shown for the P phase at coexistence with the G phase and for a temperature T = 0.815; α = x (red square), α = y (green circle), α = z (blue triangle).

FIG. 6.

Plots of components Pαz (α = x, y, or z) of the polarization vector Pz as functions of position z for a slit pore of width sz = 5.0, εw = 2.0, and κ = 0.2. Data are shown for the P phase at coexistence with the G phase and for a temperature T = 0.815; α = x (red square), α = y (green circle), α = z (blue triangle).

Close modal

The data plotted in Fig. 6 also illustrate the robustness of our numerical procedure to solve Eqs. (3.33), (3.36), and (3.37) simultaneously. To initiate the iterative solution of these equations, we initially set P = 1 for the P phase with a direction of P such that Px=Py=Pz=1/3 where Eqs. (4.1a), (4.1b), and (4.1c) are used to compute initial values of αLM. The simultaneous iterative solution of Eqs. (3.33), (3.36), and (3.37) lowers P = Px ≃ 0.80 somewhat and causes P=Pxêx to point in the x-direction because of our current choice of gθϕ [see Eq. (4.3)]. We have also tried a couple of other initial choices of P without observing any significant impact on the plots shown in Fig. 6.

Moreover, plots in Fig. 7 show that varying the strength of the fluid-wall attraction changes the orientation distribution, unlike in the case of nonspecific anchoring (see Fig. 4). First, the plots in Fig. 7 illustrate that even the G phase now possesses a certain degree of order in the vicinity of the walls. Towards the center of the pore, this order diminishes and becomes smallest around the middle of the pore because here the fluid-wall attraction is weakest. In the P phase order near the walls still decreases on account of the reduction of the bulk mean field.

FIG. 7.

As Fig. 6, but for Pxz as a function of position z between the solid walls (εw = 2.0) and T = 0.85; open symbols G phase, filled symbols P phase.

FIG. 7.

As Fig. 6, but for Pxz as a function of position z between the solid walls (εw = 2.0) and T = 0.85; open symbols G phase, filled symbols P phase.

Close modal

However, because of anchoring this reduction is compensated for to some extent by the strength of the fluid-wall attraction. Consequently, one observes an upward shift of the curves shown in Fig. 7 with increasing εw that is most pronounced near the walls. If εw is sufficiently large, the plot of Px may even exhibit a nonmonotonic variation as zsz/2.

2. Hybrid anchoring

After discussing the impact of homogeneous anchoring in Sec. IV B 1, we now focus on hybrid anchoring for which the anchoring of molecules at both walls differs. We realize hybrid anchoring in two ways. First, we consider anchoring such that molecules at the lower wall tend to align with the +x-axis whereas those near the upper wall align preferentially with the −x-direction. We realize this hybrid anchoring scenario via the anchoring functions

gθϕ1ω=8π3ReY11ω,
(4.4a)
gθϕ2ω=8π3ReY11ω.
(4.4b)

In the second hybrid anchoring scenario, we replace Eq. (4.4b) by

gθϕ2=8π3ImY11ω
(4.5)

so that at the upper wall, molecules align preferentially with the +y-axis.

Plots in Fig. 8 illustrate the impact of hybrid anchoring on the phase behavior of the confined fluid. In general, hybrid anchoring causes a depression of the GI critical point compared with a homogeneous anchoring scenario. This reflects that for hybrid anchoring the net attraction experienced in a fluid phase is generally lower. This effect turns out to be more pronounced for antiparallel [see Eqs. (4.4)] as opposed to orthogonal hybrid anchoring [see Eqs. (4.4a) and (4.5) and Fig. 8]. It is then not surprising that the fluid confined between hard walls has the lowest GI critical temperature of all the phase diagrams plotted in Fig. 8.

FIG. 8.

As Fig. 4, but for hybrid anchoring; (blue triangle down) antiparallel anchoring [see Eqs. (4.4)], (green triangle) orthogonal anchoring [see Eqs. (4.4a) and (4.5)]. For comparison, we also show phase diagrams for hard walls (black square) and homogeneous anchoring (red circle) [see Eq. (4.3)]; sz = 5.0, εw = 2.0, and κ = 0.2. Horizontal lines demarcate two isotherms (see Fig. 9).

FIG. 8.

As Fig. 4, but for hybrid anchoring; (blue triangle down) antiparallel anchoring [see Eqs. (4.4)], (green triangle) orthogonal anchoring [see Eqs. (4.4a) and (4.5)]. For comparison, we also show phase diagrams for hard walls (black square) and homogeneous anchoring (red circle) [see Eq. (4.3)]; sz = 5.0, εw = 2.0, and κ = 0.2. Horizontal lines demarcate two isotherms (see Fig. 9).

Close modal

For reasons already emphasized in Sec. IV B 1 only the fluid confined between hard walls exhibits a critical end point. Nevertheless, that part of the coexistence curve corresponding to the P phase and antiparallel hybrid anchoring exhibits again a shoulder as a vestige of the critical end point. The phase diagrams for the other two anchoring scenarios presented in Fig. 8 lack a similar shoulder. Together with the depression of the GI critical point, we conclude that the walls stabilize the confinement in the direction no anchoring (hard wall) → antiparallel hybrid anchoring → orthogonal hybrid anchoring → homogeneous anchoring.

Corresponding to the phase diagrams presented in Fig. 8 we present plots of the total local polarization P in Fig. 9. We selected two isotherms T = 0.815 and 0.930. For the former we have GP phase coexistence in all the four cases considered in Fig. 9(a) whereas in the latter case, where we are just slightly above the temperature of the hard-wall critical end point, we have GI phase coexistence for the fluid confined between hard walls. From the phase diagrams in Fig. 8 the situation is not so clear as far as both homogeneous and heterogeneous anchorings are concerned.

FIG. 9.

As Fig. 6, but for the total polarization Pz as a function of position z. Data correspond to the phase diagrams presented in Fig. 8 for temperatures T = 0.815 (a) and T = 0.930 (b).

FIG. 9.

As Fig. 6, but for the total polarization Pz as a function of position z. Data correspond to the phase diagrams presented in Fig. 8 for temperatures T = 0.815 (a) and T = 0.930 (b).

Close modal

In Fig. 9(a) the fluid between hard walls P is a monotonically decreasing function as z changes from 0 to sz/2. In the three other cases, the walls compensate the confinement-induced reduction of the bulk mean field which causes P to rise in the immediate vicinity of the walls. Order is highest for homogeneous anchoring followed by orthogonal hybrid anchoring.

Of the two hybrid anchoring scenarios, the antiparallel one exhibits the smaller degree of order. This can be ascribed to the fact that in this case anchoring at the walls requires molecules to change their orientation by an angle of π as one moves from one wall to the other whereas for the orthogonal hybrid anchoring only a change of the molecular orientation by π/2 is required. Given the size of the pore orientational frustration of the molecules is therefore larger for antiparallel hybrid anchoring and consequently P is lower than for orthogonal hybrid anchoring.

At the higher T = 0.930, P for the fluid between hard walls vanishes. This can be understood from the corresponding phase diagram plotted in Fig. 8 where one can see that T = 0.930 is above the temperature of the critical end point. For antiparallel hybrid anchoring, a critical end point no longer exists. However, T = 0.930 is in the temperature range where the phase diagram in Fig. 8 exhibits a weak shoulder on the higher density side as a vestige of the former critical end point. Therefore, it seems consistent that P for antiparallel hybrid anchoring is still substantial but decays as one moves towards z = 0 [see Fig. 9(b)]. The cusp in the plot of P at z = 0 can be understood because here P=Px and Px changes sign at z = 0 (see below). As one can also see from Fig. 9(b), P for homogeneous and orthogonal hybrid anchoring remains nearly constant across the entire pore.

To analyze variations in local order in greater detail, we now focus on individual components of P and begin with antiparallel hybrid anchoring [see Eqs. (4.4)]. As one can see from Fig. 10, Px decreases monotonically from positive (alignment with +x-axis) to eventually negative values (alignment with −x-axis). Moreover, Px is point symmetric with respect to the pore’s midpoint at z = 0 as the plot in Fig. 10 shows, that is, Pxz=Pxz.

FIG. 10.

Plots of Pαz [α = x (red square), α = y (green circle), α = z (blue triangle)], and Pz (black triangle down) as functions of position z between the slit-pore walls for antiparallel anchoring (sw = 5.0, εw = 2.0, κ = 0.2).

FIG. 10.

Plots of Pαz [α = x (red square), α = y (green circle), α = z (blue triangle)], and Pz (black triangle down) as functions of position z between the slit-pore walls for antiparallel anchoring (sw = 5.0, εw = 2.0, κ = 0.2).

Close modal

The other two components of the polarization vector P exhibit a completely different behavior. First, as one can see from Fig. 10, Py = Pz. Both components exhibit axial symmetry with respect to z = 0 unlike their counterpart Px, that is, Pαz=Pαz (α = y or z). Second, Py and Pz pass through a maximum at z = 0 and are smallest at the walls. Hence, in changing the alignment from the +x- to the −x-direction, the fluid’s orientation “escapes” to the y–z plane to a certain degree. This effect is ascribed to orientational frustration of molecules around the pore’s midsection.

A third noteworthy feature shown in Fig. 10 is that despite the strong dependence of all three individual components of P on position, P itself appears to be remarkably constant across the pore and exhibits only very minute variations. Together, P and the three components of P indicate that the confined spin fluid is highly ordered but in a strongly inhomogeneous and anisotropic fashion.

Because of the variation of components of P across the pore illustrated by plots in Fig. 10, it seems interesting to take a closer look at the impact of other fine details of the fluid-wall interaction on the local order in a confined Heisenberg fluid. For example, if one reduces the range of the fluid-wall interaction by making the inverse Debye length larger a comparison of the plots in Figs. 10 and 11(a) clearly shows that qualitatively not much changes. Even quantitatively the data are almost unaffected by the range of the fluid-wall potential. Apparently only a small film near the walls needs to be ordered and then this order propagates into portions of the confined fluid that are more remote from the walls without them directly interfering any more. The thickness of this film, controlled by κ, apparently does not matter too much.

FIG. 11.

As Fig. 10, but for κ = 2.0 (a) and sz = 15.0 (b).

FIG. 11.

As Fig. 10, but for κ = 2.0 (a) and sz = 15.0 (b).

Close modal

On the contrary, if one maintains the relatively long-range fluid-wall potential as in Fig. 10 but makes the pore three times wider as in Fig. 11(b) the plots exhibit some subtle changes. First, one notices a more sigmoidal variation of Px as one moves across the pore, whereas Pα (α = y or z) exhibit a bell-like shape.

Second, near the walls Pα is much smaller for the larger sz indicating that in the vicinity of the walls orientational frustration of molecules is lower if the pore is wider. This is because in a sufficiently wide pore molecules in the vicinity of one wall interact predominantly with that wall and not with the other one; for smaller wall separation the total external field exerted on a molecule is due to the combined effect of both walls. As in Fig. 10, however, P is nearly constant in Fig. 11(b) but its value is marginally larger.

A completely different orientational structure is observed if one replaces Eq. (4.4b) by Eq. (4.5), that is, for orthogonal hybrid anchoring. Accordingly, Px is relatively large at the lower substrate and then decreases monotonically as zsz/2; Py shows the opposite behavior (see Fig. 12). At the same time, Pz vanishes identically everywhere. Again one notices symmetry between the curves plotted in Fig. 12 which can be expressed through the relation

Pxz=Pyz,szzsz.
(4.6)

In addition, P is roughly constant across the pore very similar to the antiparallel wall anchoring shown in Fig. 10.

FIG. 12.

As Fig. 10, but for anchoring along the +x-axis at the lower and along the +y-axis at the upper substrate.

FIG. 12.

As Fig. 10, but for anchoring along the +x-axis at the lower and along the +y-axis at the upper substrate.

Close modal

The variation of Px and Py illustrated in Fig. 12 indicates that in the P phase molecules are rotating in the xy plane. Thus, it is conceivable that this rotation may lead to helicoidal structures evolving around the z-axis under favorable conditions. We defer an in-depth exploration of these wall-induced, “quasi-cholesteric” structures in nanoconfined Heisenberg fluids to part II of this work.

In this work we develop a DFT for a confined fluid composed of classical, three-dimensional Heisenberg “spins” with ferromagnetic spin-spin coupling. Within our DFT pair correlations are treated at the modified mean-field level at which they are approximated through an orientation-dependent Mayer f-function. Even though this approach overestimates the correlation length at low temperatures, it is generally believed4,8,14,22–24,26,33,34 to be an improvement over a simple mean-field treatment in which pair correlations are ignored altogether.

The fluid is placed between two planar walls that have a hard core plus a superimposed soft attraction. The attractive part of the fluid-wall interaction is modeled by a Yukawa potential. This particular form has been chosen because it allows us to control the range of the attraction via the Debye screening length. For the confined fluid we ignore packing effects, that is, we take the singlet density to be composed of a constant number density times a position-dependent orientation distribution function following the work of Gramzow and Klapp.22 

We focus on model parameters for which the bulk phase diagram is of type I according to the classification scheme of Schoen et al. (see Fig. 1 of Ref. 8). This is the simplest topography of a phase diagram where besides a G phase, an I as well as a P phase forms.

Even though our model is three-dimensional, there is an interesting link to a study of Blöte et al.38 of a two-dimensional Heisenberg fluid. Their work should be closely related to our’s but only in the limit of very narrow slit pores. Blöte et al.38 observe that in two dimensions the classical Heisenberg fluid does not exhibit any order-disorder phase transition for T > 0 in agreement with the celebrated Mermin-Wagner theorem.39 However, if the Hamiltonian is a sufficiently nonlinear function of the spin-spin interaction, a discontinuous order-disorder phase transition appears and preempts a Kosterlitz-Thouless transition. Nonetheless, in the ordered phase, order remains short-range.

Another study, even more closely related to the present one, has been carried out by Marguta et al.36 By means of computer simulations and a mean-field theory these authors investigated the Lebwohl-Lasher (lattice) model of a liquid crystal confined to a slit pore. In this model, the orientation dependence of the intermolecular interaction potential is, however, quite different from ours in that it is of the Mayer-Saupe form and therefore invariant with respect to the transformation ûû [cf. Eq. (2.5)]. This invariance reflects the “head-tail” symmetry of many liquid-crystalline systems which, at the Mayer-Saupe level, causes the leading (nontrivial) term in the expansion Eq. (2.4) to be Φ22026 rather than Φ110 as in this work. Interestingly, Marguta et al.36 allow for hybrid anchoring in the sense that the strength of the interaction between a liquid-crystal molecule and at one wall of the slit pore may differ from that at the other. However, within the framework of the present study, this would again most likely cause the rather trivial quasi-bulk effects discussed above.

In our work, confinement effects arise in two different ways. First, because of the mere presence of the walls and regardless of whether they are hard or soft the mean-field exerted on a molecule near the walls is always reduced. This causes the orientation distribution function to decline near the walls indicating that portions of the fluid in their vicinity are always less ordered than its inner parts. This has been observed in a qualitatively similar fashion by Gramzow and Klapp for their confined Stockmayer fluid (see Fig. 9 of Ref. 22).

However, Gramzow and Klapp also observe very weak maxima in the x-component of their polarization vector that are not observed here.22 Presumably, these maxima are caused by the fact that Gramzow and Klapp have chosen a state point deep in the ferroelectric polar phase whereas we are exclusively working under coexistence conditions.

If the walls have a superimposed attractive tail the soft attraction of the fluid-wall interaction potential is completely irrelevant at coexistence as long as the fluid-wall interaction is not explicitly orientation-dependent. This is because then the fluid-wall interaction potential becomes a (position-dependent) prefactor in the orientation distribution function that vanishes in the Euler-Lagrange equations to be solved under the explicit assumption of phase coexistence. This fine point was apparently not noted explicitly by Gramzow and Klapp.22 

The fluid-wall interaction enters the picture only in a less trivial fashion if the fluid-wall potential is explicitly orientation-dependent which we realize here through anchoring functions. One notices immediately that an orientation-dependent fluid-wall potential remains in the integrand in Eq. (3.29). Hence, no trivial cancellation occurs between denominator and numerator of the logarithmic term in the integrand in Eq. (3.36) because this weighting is different between G and P phases on account of the inequality ρGρP. The same is true for the integrands in Eq. (3.37). In this sense an orientation-dependent fluid-wall potential has the same formal status within our DFT approach as that of a nonlocal external magnetic field imposed on a bulk Heisenberg fluid.37 

Last but not least, we notice that quite a bit of experimental work on systems related to the present one has already been carried out as pointed out in the review by Cabuil.40 Experimentally one is mostly concerned with thermodynamic properties of nanoparticle suspensions in confinement whereas the present work is devoted to the phase behavior of pure magnetic nanoparticles.

In the future we also plan to abolish the approximation to the singlet density used here by a more realistic one that takes into account a locally varying density where the locality is a reflection of short-range correlations between hard spheres in the reference fluid. To that end one needs to give up the Carnahan-Starling equation of state in favor of a fundamental measure-theoretical treatment of Fhs in the spirit of the work by Szalai and Dietrich.21 

We are grateful to the Deutsche Forschungsgemeinschaft for financial support under the auspices of the International Graduate Research Training Group 1524 “Self-assembled soft matter nanostructures at interfaces” and to Professor Dr. Sabine H. L. Klapp for discussions. In addition, K.E.G. thanks the U.S. National Science Foundation for support of this work through Grant No. CBET-1160151. M.S. is grateful for support through the grant Scho 525 9-1.

In this appendix, we will give a brief account of the calculation of the bulk mean field uLM [see Eq. (3.21)]. A key ingredient of this calculation is the coefficients fLL0 in the expansion of the Mayer f-function. We begin by expanding the anisotropic part of the Mayer f-function in Eq. (3.14) and obtain

fL1L2Lr12=4π2L+1expβφisor12{I0βAr12I1+β2A2r122!I2β3A3r123!I3±}1.
(A1)

We truncate this expansion after the third-order term. In Eq. (A1) [cf. Eqs. (2.7) and (2.8)]

Ar12=4π3/23εHφattr12,
(A2a)
In=dω1dω2dωΦ110nΦL1L2L*.
(A2b)

The first few members of In can be computed with moderate effort. For example, using Y00=1/4π25 and Eqs. (2.6) and (3.15) one can easily verify that

I0=4π3/2dω1dω2dωΦ000ΦL1L2L*=4πδL10δL20δL0.
(A3)

Moreover, Eq. (3.15) leads directly to

I1=14πδL11δL21δL0.
(A4)

The calculation of I2 turns out to be somewhat more demanding as it involves the product of three rotational invariants in the integrand in Eq. (A2b). Using the product rule for rotational invariants [see, for example, Eq. (B8) of Ref. 14] this triple product can be reduced to a product of just two rotational invariants which then allows one to apply Eq. (3.15). We shall sketch this analysis below. We begin by considering (see Appendices A.3 and A.5 of Ref. 25)

ΦL1L2LΦ110=1L1+L2+L+1+1+04π3/2μ1μ2μ[2L+1220+122L1+12L2+1×21+121+12μ1+12μ2+1]1/2L11μ1000L11μ2000L0μ000L1L2L110μ1μ2μΦμ1μ2μ,
(A5)

where the terms in parentheses on the fourth line are Wigner 3j symbols whereas the term in braces denotes a Wigner 9j symbol.25 Inserting this lengthy expression into Eq. (A2b) (n = 2) one can employ Eq. (3.15) to realize that this leads to μ1 = μ2 = 1 and μ = 0. Because (at least) one entry in the resulting 9j symbol is zero, we can employ Eq. (A.292) of Ref. 25 to reduce the 9j to a 6j symbol which can then be further reduced to a 3j symbol by means of Eq. (A.285). From this scheme we finally arrive at

I2=4π5/232L+11/2L110002δL0δL1L2,
(A6)

where L1 = L2 = L′.

Notice now that from Eq. (A2b), I3 involves a quadruple product of rotational invariants. This product can be reduced by applying exactly the same manipulations that finally led to the expression in Eq. (A6), only this time the product rule has to be applied twice. The final result is

I3=4π433μ2L+12μ+121/2×L1μ000μ110002δL0δL1L2.
(A7)

At this stage it is useful to notice that the Wigner 3j symbols in Eqs. (A6) and (A7) are proportional to Clebsch-Gordan coefficients [see Eq. (A.139) of Ref. 25]. From the parity selection rule of the Clebsch-Gordan coefficients [see Eq. (A.155) of Ref. 25] this implies that the 3j symbol in Eq. (A6) vanishes unless L′ is zero or even. Similarly, the second 3j symbol in Eq. (A7) vanishes unless μ′ vanishes or is even. This, in turn, implies that the first 3j symbol in Eq. (A7) is nonzero only if here L′ is odd.

Because all the expressions in Eqs. (A3), (A4), (A6), and (A7) contain the Kronecker symbol δL0, fL1L2LfLL0 where we replaced L′ again by L to simplify the notation. Hence, from Eq. (A1) we finally have

fLL0=4π3/2expβφiso{δL0+βAδL1+3β2A22!2L+11/2L110002+33β3A33!2L+11/2μ2μ+1L1μ0002×μ110002}4π3/2δL0,
(A8)

where A=A4π3/2. With Eq. (A8) we can in principle compute the full mean field uLM from Eq. (3.21). However, in the main body of the paper we restrict the treatment to the terms L = 0 and L = 1. In this case only the much simpler expansion coefficients

f0004π3/2=expβφiso1+β2A221,
(A9)
f1104π3/2=expβφisoβA+310β3A3
(A10)

are needed. To arrive at Eqs. (A9) and (A10) we used again the parity selection rule [see Eq. (A.155) of Ref. 25]. It suggests that for L = 0 the 3j symbol involved in the second-order term proportional to β2A2 in Eq. (A8) is nonzero. However, μ appearing in the third-order term has to be odd and even at the same time for both 3j symbols to be nonzero simultaneously. Obviously, this is a contradiction so that this term vanishes for L = 0.

For L = 1, however, the parity selection rule suggests that the 3j symbol in the second-order term vanishes and that μ = 2n (n ∈ ℕ) in the third-order term, where ℕ is the set of nonnegative integers. Applying in addition the triangle inequality [see Eq. (A.131) of Ref. 25] one then realizes that only the summands for μ = 0, 2 survive in Eq. (A8).

Finally, upon insertion of Eqs. (A9) and (A10) into Eq. (3.21) the remaining integrations over da12 are carried out numerically using a simple trapezoidal rule.

1.
P. C.
Hemmer
and
D.
Imbro
,
Phys. Rev. A
16
,
380
(
1977
).
2.
I. P.
Omelyan
,
W.
Fenz
,
I. M.
Mryglod
, and
R.
Folk
,
Phys. Rev. Lett.
94
,
045701
(
2005
).
3.
E.
Lomba
,
J. J.
Weis
, and
G.
Stell
,
Phys. Rev. E
50
,
3853
(
1994
).
4.
J. M.
Tavares
,
M. M. T.
da Gama
,
P. I. C.
Teixeira
,
J. J.
Weis
, and
M. J. P.
Nijmeijer
,
Phys. Rev. E
52
,
1915
(
1995
).
5.
A.
Oukouiss
and
M.
Baus
,
Phys. Rev. E
55
,
7242
(
1997
).
6.
E.
Lomba
,
J. J.
Weis
, and
C. F.
Tejero
,
Phys. Rev. E
58
,
3426
(
1998
).
7.
F.
Lado
,
E.
Lomba
, and
J. J.
Weis
,
Phys. Rev. E
58
,
3478
(
1998
).
8.
M.
Schoen
,
S.
Giura
, and
S. H. L.
Klapp
,
Phys. Rev. E
89
,
012310
(
2014
).
9.
V.
Bashtovoy
and
B. M.
Berkovsky
, in
Magnetic Fluids and Applications
(
Begell House
,
New York
,
1996
), Chap. 4.4.1.
10.
S.
Miyake
and
S.
Takahashi
,
ASLE Trans.
28
,
461
(
1985
).
12.
Š.
Jakabský
,
M.
Lovs
,
A.
Mockovčiaková
, and
S.
Hredzk
,
J. Radioanal. Nucl. Chem.
246
,
543
(
2000
).
13.
N. A.
Brusentsov
,
V. V.
Gogosov
,
T. N.
Brusentsova
,
A. V.
Sergeev
,
N. Y.
Jurchenko
,
A. A.
Kuznetsov
,
O. A.
Kuznetsov
, and
L. I.
Shumakov
,
J. Magn. Magn. Mater.
225
,
113
(
2001
).
14.
B.
Groh
and
S.
Dietrich
,
Phys. Rev. E
50
,
3814
(
1994
).
15.
S. H. L.
Klapp
,
J. Phys.: Condens. Matter
17
,
R525
(
2005
).
16.
S. H.
Lee
,
J. C.
Rasaiah
, and
J. B.
Hubbard
,
J. Chem. Phys.
85
,
5232
(
1986
).
17.
J.
Jordanovic
and
S. H. L.
Klapp
,
Phys. Rev. Lett.
101
,
038302
(
2008
).
18.
S. H. L.
Klapp
and
M.
Schoen
,
J. Chem. Phys.
117
,
8050
(
2002
).
19.
S. H. L.
Klapp
and
M.
Schoen
,
J. Mol. Liq.
109
,
55
(
2004
).
20.
V.
Ballenegger
and
J.-P.
Hansen
,
J. Chem. Phys.
122
,
114711
(
2005
).
21.
I.
Szalai
and
S.
Dietrich
,
Eur. Phys. J. E
28
,
347
(
2009
).
22.
M.
Gramzow
and
S. H. L.
Klapp
,
Phys. Rev. E
75
,
011605
(
2007
).
23.
P.
Frodl
and
S.
Dietrich
,
Phys. Rev. A
45
,
7330
(
1992
).
24.
P.
Frodl
and
S.
Dietrich
,
Phys. Rev. E
48
,
3203
(
1993
).
25.
C. G.
Gray
and
K. E.
Gubbins
,
Theory of Molecular Fluids
(
Clarendon Press
,
Oxford
,
1984
), Vol.
1
.
26.
S.
Giura
and
M.
Schoen
,
Phys. Rev. E
90
,
022507
(
2014
).
27.
J. A.
Castellano
,
Mol. Cryst. Liq. Cryst.
94
,
33
(
1983
).
28.
A. A.
Sonin
,
The Surface Physics of Liquid Crystals
(
Gordon and Breach
,
Amsterdam
,
1995
).
29.
C. G.
Gray
,
K. E.
Gubbins
, and
C. G.
Joslin
,
Theory of Molecular Fluids
(
Oxford University Press
,
Oxford
,
2011
), Vol.
2
.
30.
N. F.
Carnahan
and
K. E.
Starling
,
J. Chem. Phys.
51
,
635
(
1969
).
32.
J.-P.
Hansen
and
I. R.
McDonald
, in
Theory of Simple Liquids
, 3rd ed. (
Elsevier
,
Burlington
,
2006
), Chap. 3.4.
33.
G. M.
Range
and
S. H. L.
Klapp
,
Phys. Rev. E
69
,
041201
(
2004
).
34.
S.
Giura
,
B. G.
Márkus
,
S. H. L.
Klapp
, and
M.
Schoen
,
Phys. Rev. E
87
,
012313
(
2013
).
35.
K.-C.
Ng
,
J. Chem. Phys.
61
,
2680
(
1974
).
36.
R. G.
Marguta
,
Y.
Martínez-Ratón
,
N. G.
Almarza
, and
E.
Velasco
,
Phys. Rev. E
83
,
041701
(
2011
).
37.
S. M.
Cattes
,
S. H. L.
Klapp
, and
M.
Schoen
,
Phys. Rev. E
91
,
052127
(
2015
).
38.
H. W.
Blöte
,
W.
Guo
, and
H. J.
Hilhorst
,
Phys. Rev. Lett.
88
,
047203
(
2002
).
39.
N. D.
Mermin
and
H.
Wagner
,
Phys. Rev. Lett.
17
,
1133
(
1966
).
40.
V.
Cabuil
,
Curr. Opin. Colloid Interface Sci.
5
,
44
(
2000
).