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.

## I. INTRODUCTION

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 fluid^{3–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 treatment^{13} 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 fluids^{3–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 films^{18} and later on the nature of ferroelectric states in these systems^{19} 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 Dietrich^{21} 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 Klapp^{22} 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.

## II. MODEL

### A. Fluid phase

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

which we decompose formally into an isotropic contribution *φ*_{iso} and an anisotropic one represented by *φ*_{anis}. In Eq. (2.1), *r*_{12} = *r*_{1} − *r*_{2} is the distance vector between the centers of mass of a pair of molecules located at *r*_{1} and *r*_{2}, 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

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

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 to^{25}

where Φ_{l1l2l} is a rotational invariant, *l*_{1}, *l*_{2}, and *l* are nonnegative integers, *ω* is a set of Euler angles describing the orientation of $r\u030212=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} = *φ*_{att}^{25} irrespective of *l*_{1}, *l*_{2}, and *l*.

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

where *γ* is the angle between the spins with orientations specified by $u\u0302\omega 1$ and $u\u0302\omega 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

where *C* is a Clebsch-Gordan coefficient, $Yl\u2032m\u2032$ 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 that^{26}

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

Two comments apply at this stage. First, it is clear that our molecules have uniaxial symmetry such that $\omega i=\theta i,\varphi i$ (*i* = 1, 2) is a subset of two out of three Euler angles. Second, upon inverting one of the two spins by replacing $\omega i\u2192\omega i\u2032=\u2212\omega i=\pi \u2212\theta i,\pi +\varphi 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\u0302ij\u2192r\u0302ij\u2032=\u2212r\u0302ij$. Hence, our model fluid consists of polar molecules and is therefore capable of forming ordered polar phases.

### B. Confinement

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

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

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

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), $0\u2264g\theta \varphi k\omega \u22641$ 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

where

and $\theta k$ and $\varphi 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 *x*–*y* plane characterized by $\theta 1=\theta 2=\pi 2$. However, in general, $\varphi 1$ may be different from $\varphi 2$ which allows us to implement hybrid anchoring scenarios where homogeneous anchoring (i.e., $\varphi 1=\varphi 2$) is included as a special case. A sketch of our model system is presented in Fig. 1.

## III. ELEMENTS OF DENSITY FUNCTIONAL THEORY

For the model introduced in Sec. II, we can associate thermodynamic states at equilibrium with minima of the grand potential^{29}

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 $\rho r,\omega $. In Sec. III A, we will briefly sketch how these thermodynamic equilibrium states can be identified.

### A. Basic expressions and approximations

As a first approximation, we immediately introduce

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

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, $\alpha =1/4\pi $ 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 *s*_{z} 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

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

where *β* = 1/*k*_{B}*T* (*k*_{B} 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

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

for the free-energy density *f*^{hs} 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 as^{31,32}

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 *r*_{1} and *r*_{2} 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 *r*_{12}.

We approximate *g* at the so-called modified mean-field level^{4,8,14,23,24,26,33,34} as

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,\gamma i$ is the projection of *r*_{i} onto the *x*–*y* plane and *γ _{i}* is the polar angle. Thus, it is useful to transform variables such that $a1,\gamma 1\u2192a1\u2032,\gamma 1\u2032=a1,\gamma 1$ and $a2,\gamma 2\u2192a12,\gamma 12$.

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

where Θ is the Heaviside function and

### B. Notes on integrations over orientations

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

This expansion introduces a fourfold summation over $l1,m1$ and $l2,m2$ into Eq. (3.10) where corresponding pairs of these integers $l\u2032,m\u2032$ 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

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

because $\Phi l1l2l$ is a complete set of orthogonal basis functions, that is,

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

where *θ* is the polar angle specifying the orientation of $r\u030212$ in spherical coordinates and $M\xaf=\u2212M$.

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

so that only $YL0\theta $ remains. Second, the selection rule *M*_{1} + *M*_{2} = *M* for nonzero Clebsch-Gordan coefficients has been used from which *M*_{1} = − *M*_{2} = *M* follows. However, this implies also *L*_{1} = *L*_{2} = *L*′. Third, we utilize [see Eq. (A.3) of Ref. 25]

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].

where the coefficients

can only depend on $z12=z1\u2212z2$ because both $r12=a122+z122$ and cos*θ* = *a*_{12}/*r*_{12} depend on *z*_{12}.

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 *r*_{12} = *σ*. The temperature range will be such that *k*_{B}*T*/ε 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

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

### C. Thermodynamic equilibrium

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

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

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*).

where

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

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 $\rho ln\rho \Lambda 5m/I$ arises in the expressions for *βf*^{hs} [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,

*u*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.”

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

where the angle-averaged orientation distribution function is given by

and

### D. Phase coexistence

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 $\alpha LM\u2032$ and $\alpha LM\u2033$ in general.

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

In addition, at coexistence ′ and ″ must satisfy

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

where from Eq. (3.22)

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), $\alpha 00=1/4\pi $ changed variables according to $z1\u2192z1\u2032=z1$ and *z*_{2} → *z*_{12} and performed the one trivial integration over $dz1\u2032$ 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

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

which is again independent of the unknown chemical potential and a term proportional to $ln\Lambda 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

## IV. RESULTS

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 Klapp

^{22}and introduce the polarization $Pz=Pxze\u0302x+Pyze\u0302y+Pzze\u0302z$ such that the total polarization $Pz=Pz\u22641$. 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

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

### A. Non-specific anchoring

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 ε/*k*_{B}.

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 $\phi wkz=\phi hwz$. One then realizes from Eq. (3.30) that for hard walls, the somewhat simpler expression

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 *T*_{c} ≃ 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).

Confining the Heisenberg fluid to slit pores of variable pore widths *s*_{z} does not affect the overall topology of the phase diagram. It merely serves to shift the critical temperature down monotonically with decreasing *s*_{z} 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 *s*_{z} 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 *s*_{z} 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.

As anticipated, $Pz\u22430.57$ does not depend on *z* for the bulk system. In a wide pore of *s*_{z} = 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 *s*_{z} = 7.5 and *s*_{z} = 5.0 show that $Pz$ depends on *z* across the entire pores and passing through maxima at *z* = 0.

The decline of $Pz$ as $z\u2192sz/2$ can again be described to the quasi-bulk effect. Because the phase diagrams in Fig. 2 approach the bulk limit as *s*_{z} → ∞, 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\theta \varphi k=1$ in Eq. (2.9). It is then clear that in Eq. (3.30), $\phi wkz,\omega \u2192\phi wkz$. Hence, when performing integrations over d*ω* in Eq. (3.37), $exp\u2212\beta \u2211k=12\phi 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 *s*_{z} 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.

### B. Explicit anchoring

#### 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\theta \varphi 1=g\theta \varphi 2=g\theta \varphi $. More specifically, we take

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 $\phi wk$ to be explicitly orientation-dependent, no cancellation of $exp\u2212\beta \u2211k=12\phi 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 d*T*/d*ρ* changes discontinuously.

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, d*T*/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

**and the molecular orientation $u\u0302$. 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.**

*H*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 *P*_{x} varies nonmonotonically across the pore whereas *P*_{y} = *P*_{z} ≃ 0. In other words, *g*_{θϕ} in Eq. (4.3) completely dominates the direction of ** P** if the pore is small enough.

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 $\alpha LM$. The simultaneous iterative solution of Eqs. (3.33), (3.36), and (3.37) lowers

*P*=

*P*

_{x}≃ 0.80 somewhat and causes $P=Pxe\u0302x$ 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

**without observing any significant impact on the plots shown in Fig. 6.**

*P*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.

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 *P*_{x} may even exhibit a nonmonotonic variation as $z\u2192sz/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

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

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.

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.

In Fig. 9(a) the fluid between hard walls *P* is a monotonically decreasing function as $z$ changes from 0 to *s*_{z}/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 *P*_{x} 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,

*P*

_{x}decreases monotonically from positive (alignment with +

*x*-axis) to eventually negative values (alignment with −

*x*-axis). Moreover,

*P*

_{x}is point symmetric with respect to the pore’s midpoint at

*z*= 0 as the plot in Fig. 10 shows, that is, $Pxz=\u2212Px\u2212z$.

The other two components of the polarization vector ** P** exhibit a completely different behavior. First, as one can see from Fig. 10,

*P*

_{y}=

*P*

_{z}. Both components exhibit axial symmetry with respect to

*z*= 0 unlike their counterpart

*P*

_{x}, that is, $P\alpha z=P\alpha \u2212z$ (

*α*= y or z). Second,

*P*

_{y}and

*P*

_{z}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

**indicate that the confined spin fluid is highly ordered but in a strongly inhomogeneous and anisotropic fashion.**

*P*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.

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 *P*_{x} 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 *s*_{z} 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, *P*_{x} is relatively large at the lower substrate and then decreases monotonically as *z* → *s*_{z}/2; *P*_{y} shows the opposite behavior (see Fig. 12). At the same time, *P*_{z} vanishes identically everywhere. Again one notices symmetry between the curves plotted in Fig. 12 which can be expressed through the relation

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

The variation of *P*_{x} and *P*_{y} illustrated in Fig. 12 indicates that in the P phase molecules are rotating in the *x*–*y* 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.

## V. DISCUSSION AND CONCLUSIONS

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 believed^{4,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 $u\u0302\u2192\u2212u\u0302$ [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 Φ_{220}^{26} 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}*≠

*ρ*. 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.

^{P}^{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}

## Acknowledgments

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.

### APPENDIX: EXPANSION COEFFICIENTS *f*_{LL0} OF THE ORIENTATION-DEPENDENT MAYER *f*-FUNCTION

*f*

*f*

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 *f*_{LL0} 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

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

Moreover, Eq. (3.15) leads directly to

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)

where the terms in parentheses on the fourth line are Wigner 3*j* symbols whereas the term in braces denotes a Wigner 9*j* 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 9*j* symbol is zero, we can employ Eq. (A.292) of Ref. 25 to reduce the 9*j* to a 6*j* symbol which can then be further reduced to a 3*j* symbol by means of Eq. (A.285). From this scheme we finally arrive at

where *L*_{1} = *L*_{2} = *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

At this stage it is useful to notice that the Wigner 3*j* 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 3*j* symbol in Eq. (A6) vanishes unless *L*′ is zero or even. Similarly, the second 3*j* symbol in Eq. (A7) vanishes unless *μ*′ vanishes or is even. This, in turn, implies that the first 3*j* 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}, *f*_{L1L2L} → *f*_{LL0} where we replaced *L*′ again by *L* to simplify the notation. Hence, from Eq. (A1) we finally have

where $A\u2032=A4\pi \u22123/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

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 3*j* symbol involved in the second-order term proportional to *β*^{2}*A*′^{2} in Eq. (A8) is nonzero. However, *μ* appearing in the third-order term has to be odd and even at the same time for both 3*j* 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 3*j* symbol in the second-order term vanishes and that *μ* = 2*n* (*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).