Phase transitions and ordering of confined dipolar fluids

We apply a modified mean-field density functional theory to determine the phase behavior of Stockmayer fluids in slitlike pores formed by two walls with identical substrate potentials. Based on the Carnahan-Starling equation of state, a fundamental-measure theory is employed to incorporate the effects of short-ranged hard sphere - like correlations while the long-ranged contributions to the fluid interaction potential are treated perturbatively. The liquid-vapor, ferromagnetic liquid - vapor, and ferromagnetic liquid - isotropic liquid first-order phase separations are investigated. The local orientational structure of the anisotropic and inhomogeneous ferromagnetic liquid phase is also studied. We discuss how the phase diagrams are shifted and distorted upon varying the pore width.


I. INTRODUCTION
The structure and thermodynamics of confined dipolar fluids (such as molecular liquids, ferrofluids, and electrorheological fluids) have recently attracted considerable attention. Theoretical 1,2,3 , computer simulation 4,5,6,7,8,9 , and experimental 10,11 studies have provided significant information about the orientational and spatial arrangement of dipolar particles in the vicinity of solid walls. The strongly anisotropic and long-ranged character of dipolar forces causes particular difficulties for the theoretical description of these systems. As long as one is not aiming for a quantitative description of a specific system but for general phenomena and trends, the so-called Stockmayer model has turned out to be rather useful 12 . It considers spherical particles interacting with Lennard-Jones (LJ) pair potentials plus pointlike permanent dipoles at their centre. In this paper we adopt the magnetic language, assuming that the particles carry magnetic dipole moments of strength m because the main applications we have in mind are ferrofluids. (The results are identical for ferroelectric particles.) Bulk dipolar fluids such as the Stockmayer model system may exhibit three distinct fluid phases: isotropic vapor, isotropic liquid, and ferromagnetic liquid 12,13 . (In the present study we focus on sufficiently high temperatures so that freezing is not a concern 14 .) Accordingly the following first-order phase transitions can occur: isotropic vapor -isotropic liquid, isotropic liquidferromagnetic liquid, and isotropic vapor -ferromagnetic liquid. The phase diagrams exhibit the usual liquid-gas critical point, a triple point, a tricritical point, and a line of critical points corresponding to second-order phase transitions between the isotropic and the ferromagnetic liquid. For bulk dipolar fluids these phase transitions and critical points and lines have been detected in various theoretical studies 12,13,14,15,16 .
However, to the best of our knowledge there is no theoretical study which predicts the complete global phase diagrams and the structure of confined dipolar fluid phases. 3 have studied the phase behavior and orientational structure of Stockmayer fluids confined to slit pores. Within their modified mean-field (MMF) density functional theory (DFT) they considered only the subspace of spatially homogeneous local densities throughout the slit pore and focused on the orientational structure only. This allowed them to study the influence of the confinement on the phase behavior but neglecting the inevitable spatial inhomogeneities of the phases under study. Our results will show that in the case of confined dipolar liquids, as for other fluids, there are strong spatial inhomogeneities in the vicinity of the walls. This inhomogeneous character of the liquid phases becomes more pronounced upon increasing of the density. These phenomena require a more sophisticated theory in order to describe the structure and the phase equilibria of confined dipolar systems.

Recently Gramzow and Klapp
In the present study we also use the modified meanfield (MMF) DFT approximation. However, for the description of the hard sphere reference system, instead of the homogeneous bulk Carnahan-Starling free energy expression, a fundamental measure theory free energy functional 17 is applied in order to capture the shortranged hard sphere -like correlations. Our results for the confinement induced shifts and the distortions of the phase boundaries relative to the bulk ones are compared with the corresponding Kelvin equation.
The paper is organized as follows. In Sects. 2 and 3 the model and the extension of the MMF density functional theory for confined dipolar fluids are presented and developed, respectively. In Sec. 4 the calculation of the orientational order and the physical meaning of the order parameters are discussed. Sections 5 and 6 present the details of the Euler-Lagrange and of the Kelvin equations, respectively. The results of our calculations and a discussion are given in Sec. 7. Certain important computational details are summarized in the Appendices A, B, and C.
Thus the total pair potential is given by u S (r 12 , ω 1 , ω 2 ) = u r (r 12 ) + u exc (r 12 , ω 1 , ω 2 ). (5) This decomposition lends itself to choose a hard-sphere reference system characterized by u r (r 12 ) = u hs (r 12 ) for which reliable approximations for the free energy are known. The excess part u exc (r 12 , ω 1 , ω 2 ) will be treated perturbatively in an appropriate way.
For the fluid-wall interaction we consider two different potentials. In the case of a purely hard wall we use the hard-wall potential where z is the coordinate normal to the walls and L is the distance between the surfaces of the hard, parallel walls, which we call repulsive walls. For the study of attractive walls we choose the external substrate potential where ǫ w sets the energy scale. This equation can be obtained from the assumption that both walls consist of LJ particles interacting with the fluid particles of the same size also via the LJ potential with the energy scale ǫ w . Averaging the LJ interactions over all positions of the wall particles and approximating the averaged repulsive part of the LJ potential by a hard core potential one obtains the fluid particle -wall interaction potential given by Eq. (7).

III. MODIFIED MEAN-FIELD DENSITY FUNCTIONAL THEORY FOR ONE-COMPONENT FLUIDS
For a one-component nonuniform fluid configuration ρ(r, ω) denotes the number density of dipolar particles at a point r = (x, y, z) with an orientation ω = (θ, φ) relative to a spatially fixed coordinate system. The total number density of particles, independent of orientation, is given aŝ This allows one to splitρ(r, ω) into the total number densityρ(r) and a normalized space-and angle-dependent orientational distribution functionα(r, ω): ρ(r, ω) =ρ(r)α(r, ω), dωα(r, ω) = 1. (9) Within density functional theory 18 the equilibrium density distribution ρ(r, ω; T, µ) of an inhomogeneous fluid in the presence of an external potential u w minimizes the grand canonical potential (10) and thus satisfies the Euler-Lagrange equation Here F is the Helmholtz free energy density functional of the system and µ is the chemical potential. The last term in Eq. (10) takes into account the normalization condition (Eq. (9)) for the orientational distribution function α(r, ω) by a Lagrange parameter κ(r). Due to the two distinct contributions to the particle interaction potential (Eq. (5)) the free energy functional F decomposes into an ideal gas term F id and two corresponding parts: where F ref is the reference and F exc is the excess free energy density functional. In order to describe the planar wall -fluid interfaces, we consider a slab-shaped macroscopic system with the surfaces of the slab parallel to the x − y plane. The distance between the parallel surfaces is L. Thus, apart from spontaneous symmetry breaking like freezing, the equilibrium configuration of the system is inhomogeneous only in the z direction and translationally invariant in the x and y directions. Under these circumstances the total number density ρ(r) = ρ(z) is a function of z only and the orientational profile α(r, ω) = α(z, ω) depends on z and the polar angle ω. In a coordinate system fixed in space the actual orientational distribution function can be expanded in terms of the spherical harmonics Y lm (ω) where α 00 = 1/ √ 4π due to the normalization and α * lm = (−1) m α lm because α is real. (Here and in the following m ≡ −m.) The expansion coefficients α lm (z), which can be interpreted as orientational order parameters, are related to the full distribution function via As presented below, for a slab-shaped system of transversal size L one obtains explicit expressions for the three terms of the Helmholtz free energy density functional (see Eq. (12)).
A. The free energy functional for the ideal gas For the slablike shape of the sample under consideration the ideal gas contribution has the form where A is the lateral cross-sectional area of the system, Λ is the de Broglie wavelength, and β = 1/(k B T ) is the inverse temperature.

B. The reference free energy functional
For the description of the hard sphere reference free energy functional we adopt the fundamental measure theory, which was initially proposed by Rosenfeld 19,20 : the function Φ and the weighted densities n α (z) are defined in Appendix A.

C. The excess free energy functional
We approximate the excess free energy functional in terms of the modified mean-field density functional theory 21,22,23 . Using Eq. (13) for the present slablike geometry one obtains the following expression: Here and in the following the summations over l i , m i are taken as In order to extract the explicit proportionality to the large cross-sectional area A we introduce the sum and the difference of the lateral coordinates (x 1 , y 1 ) and (x 2 , y 2 ): In terms of these new variables the integration of any function f (r 12 ) can be written as . (20) In order to proceed we replace the rectangular coordinates (x 12 , y 12 ) by polar coordinates (R 12 , φ 12 ) and take into account Eqs. (17) and (20) for the excess free energy functional: where In Eq. (22) f M denotes the Mayer function of the excess potential: where in accordance with the planar polar coordinates used for the slab-geometry Via the vector r 12 the rotationally invariant function D(ω 1 , ω 2 , ω 12 ) (Eq. (2)) depends on ϑ 12 as well. There-fore in Eq. (22) the substitutions sin ϑ 12 = R 12 have to be carried out prior to the integration over the variable R 12 . In our calculation we have expanded the exponential term in Eq. (23), containing the dipole-dipole interaction, into a Taylor series: In accordance with Eqs. (22) and (26) u l1m1l2m2 is given as a Taylor series, too: with the coefficients In Eq. (28) the coefficients A (i) l1m1l2m2 can be expressed as In Appendix B we present the nonzero coefficients A (i) l1m1l2m2 for i, l 1 , l 2 = 0, 1, 2. In the bulk limit the expressions for u l1m1l2m2 reduce to the corresponding ones in Ref. 12 : where the coefficients u (i) l1 , for l 1 = 0, 1, 2, 3, 4, are given by Eqs. (4.9)-(4.13) in Ref. 12 . It is important to note that among the coefficients u l1m1l2m2 there are terms which are nonzero in a slablike geometry, but which vanish in the bulk limit.

IV. MAGNETIZATION AND ORDER PARAMETERS
Focusing on the structure of confined dipolar fluids, here we are particularly interested in the occurrence of a spontaneous magnetization. In the slab geometry, which can be considered as the limiting case of an oblate ellipsoidal sample with vanishing aspect ratio, the demagnetization field vanishes if the fluid sample is magnetized along a spontaneously chosen direction within the xy plane. For such a type of magnetization the formation of various domains can be excluded. Therefore the confined fluid with the magnetization within the xy plane is comparable to a single domain bulk fluid in a needlelike volume with longitudinal magnetization. The free energies of these fluid samples can be mapped onto each other 12 . The definition of the local magnetization is By inserting the expansion in terms of the spherical harmonics and performing the integration we obtain For M z = 0, i.e., α 10 (z) = 0, the magnetization within the xy plane is measured by the order parameter so that the modulus M xy of the in-plane magnetization is The higher coefficients of the orientational distribution function α(z, ω) (Eq. (13)) may also be considered as local orientational order parameters. Nonzero coefficients for l = 2 indicate a certain type of ordering of the dipoles. An especially interesting quantity is the order parameter α 20 (z) that describes the orientation of dipole-axes relative to the z axis. This is related to the zz component of the local quadrupole tensor as where P 2 is the second order Legendre polynomial. Q zz = 1 corresponds to the perfect ordering of dipole axes normal to the walls while Q zz = −1/2 indicates the perfect ordering along directions in the xy plane. Randomly oriented dipoles correspond to Q zz = 0.

V. EULER-LAGRANGE EQUATION
The equilibrium configuration characterized by ρ(z) and α(z, ω) for given values of T and µ follows from minimizing the total grand canonical functional with respect toρ(z) and the functionα(z, ω): where µ ref is the chemical potential functional (first order direct correlation function) of the hard sphere reference system (see Appendix A). The corresponding set of integral equations is discussed in Ap-pendix C. As solutions of Eqs. (36) and (37) ρ(z; T, µ) and α(z, ω; T, µ) become functions of T and µ.

VI. PHASE EQUILIBRIA AND KELVIN EQUATION
The phase coexistence curves µ(T ) and the coexisting densities and orientational order parameters follow from requiring the equality of the grand potentials for the coexisting phases I and II: The functions ρ (i) (z; T, µ) and α (i) (z, ω; T, µ) (i = I, II) denote the corresponding equilibrium density and orientational distribution functions obtained from Eqs. (36) and (37). Equation (38) requires simultaneous solutions of the Euler-Lagrange equations (36) and (37) for a wide range of temperatures and chemical potentials. As mentioned in Sec. 1, we consider three kinds of two-phase equilibria: isotropic liquid -isotropic vapor, isotropic vapor -ferromagnetic liquid, and isotropic liquid -ferromagnetic liquid. In order to estimate in leading order of 1/L the location of the phase coexistence curves, the Kelvin equation (see, e.g. Ref. 24 and references therein) is used to obtain the chemical potential difference between the bulk phase points and the corresponding ones for slabs of thickness L at a given temperature and for two identical confining walls: where ρ i (i = I, II) are the bulk densities of the coexisting phases, γ wi are the corresponding wall-fluid interfacial tensions, γ I,II is the surface tension between the coexisting fluid phases I and II, and Θ I = arccos[(γ wII − γ wI )/γ I,II ] is the contact angle of a drop of fluid I with the wall (Young's equation). In the slab-shaped system, for a fluid which is in contact with two identical walls, the surface tension γ wi follows from the equilibrium density ρ(z; T, µ) and oriental distribution α(z, ω; T, µ) profiles as where p is the bulk pressure and Ω L is the grand potential per cross-sectional area A of the slablike system of transversal size L.
For calculating the expansion of the dipole-dipole interaction Mayer function we used a second order approximation. This choice restricts the expansion coefficients α lm (z) of the orientational distribution function to those with l ≤ 2 and thus |m| ≤ 2.

A. Phase diagrams
As a first step, on the basis of the corresponding bulk theory we have determined the bulk phase diagrams of Stockmayer fluids for different dipole moments. For comparison and as a reference the corresponding bulk phase diagram is depicted in all figures discussing the phase diagrams in slabs. We note that the bulk phase diagrams presented here differ numerically from those in Ref. 12 because there a fourth-order expansion with respect to the corresponding Mayer functions has been used as well as a temperature dependent hard sphere diameter for the reference system. Extending our present analysis of spatially inhomogeneous system also to fourth order in that expansion is in principle possible but requires unreasonably large technical efforts. The bulk phase diagrams within the present approximation exhibit only minor quantitative differences from those in Ref. 12 .
For the dipole moment m * = 1.5 and the wall distance L/σ = 10 Fig. 1 presents our numerical results for the phase diagrams in the chemical potential -temperature and the density -temperature planes. Figure 1 compares the bulk phase diagrams with those of the slab for repulsive or attractive walls. For describing the attractive walls (Eq. (7)) here and in the following we choose ǫ * w = 1/2. Figure 1(a) shows that below the triple point temperature T T P the isotropic vapor coexists with the ferromagnetic liquid. The confinement by repulsive walls shifts the chemical potential -temperature phase diagrams towards higher chemical potentials relative to the bulk data. In the case of attractive walls the shift is smaller and in the opposite direction. This figure also shows that for the chosen parameters the Kelvin equation provides a good estimate of these shifts. Between the temperatures T T P and T CP of the triple point temperature and of the liquid-vapor critical point, respectively, there are three possible phases: the isotropic vapor, the isotropic liquid, and the ferromagnetic fluid. The firstorder phase transition between the isotropic liquid and ferromagnetic fluid turns into a second-order phase transition above the tricritical temperature T T CP . (For large µ * and at high ρ * orientationally disordered or ferromagnetic solid phases appear 14 , which are not considered here.) We find that, compared with the bulk, the temperature range for liquid-vapor and liquid -ferromagnetic fluid coexistence is narrowed by the confinement, both for repulsive and attractive walls. According to 1(b) at low temperatures (T < T T P ) the high density ferromagnetic fluid coexists with the low density isotropic vapor. Above the triple point temperature T T P , at low and medium densities two isotropic fluids coexist, becoming identical above the liquid-vapor critical temperature T CP . At higher densities, the isotropic liquid (with the lower density) and the ferromagnetic fluid (with the higher density) are separated by first-order phase transitions which turn into second-order phase transitions above the tricritical temperature T T CP . The critical line CL of second-order phase transitions divides the liquidlike thermodynamic states into an isotropic liquid phase and into a ferromagnetic fluid phase. In comparison with the bulk coexistence curves the shift of the first-order coexistence region to lower densities is more significant for repulsive than for attractive walls. The tricritical liquid densities are lowered by the confinement in comparison with the corresponding bulk ones whereas the tricritical vapor densities are increased. The liquid-vapor critical density is lowered due to the confinement by repulsive walls and slightly increased in the case of attractive walls. At a given temperature both confinements promote the formation of the ferromagnetic phase at lower densities. These findings are in qualitative agreement with the results of Gramzow and Klapp 3 , but the shift of the coexis- tence curves of confined systems relative to the bulk ones is more significant than predicted by their homogeneous local density theory. As shown in Fig. 2, a stronger confinement (L/σ = 4) does not change the topology of the phase diagram. However, in this case the shift of the phase boundaries relative to the bulk ones are more pronounced. The temperature interval T CP − T T P , within which liquid-vapor coexistence is thermodynamically stable, becomes narrower with increasing confinement. All critical temperatures (T CP , T T CP , T T P ) are lowered upon increasing the confinement. (The lowering of the liquidvapor critical temperature for Stockmayer fluids with increasing confinement is in a qualitative aggrement with the Gibbs ensemble Monte Carlo simulation results of Richardi et al 4 .) For attractive walls the critical density of the liquid-vapor equilibria does not change significantly whereas it is reduced for repulsive walls. Comparing our results for m * = 1.5 and L/σ = 4 with those of Gramzow and Klapp 3 we note that, in contrast to their findings, in the chemical potential -temperature plane the first-order phase boundaries of the bulk and of the confined systems do not intersect.
For L/σ = 10, increasing the dipole moment from m * = 1.5 to m * = 2 changes the corresponding phase diagrams shown in Fig. 1 considerably and leads to the phase diagram displayed in Fig 3. This increase of the dipole moment causes the disappearance of the liquidvapor transitions for both the bulk and the confined systems. According to Fig. 3(a) for this dipole moment the phase boundaries of the bulk and of the confined system with attractive walls do intersect. On the basis of the Kelvin equation (Eq. (39)) this multiple crossing of the ferromagnetic fluidisotropic fluid (I = f f , II = if ) phase boundaries for the bulk and for the confined system implies a non-monotonic temperature dependence of the surface tension difference γ w,f f − γ w,if and vice versa; according to the full line in Fig. 3 Since γ f f,if is expected to decrease monotonically upon increasing temperature, these crossings amount to a nonmonotonic temperature dependence of the contact angle of the ferromagnetic fluid around 90 • .
(For Θ f f = 90 • the leading behavior of ∆µ is given by the contribution In the case of confinement with repulsive walls the ferromagnetic fluid -isotropic fluid phase boundary is generally shifted to higher values of the chemical potential as compared with its bulk counterpart. Confinement lowers the tricritical temperatures T T CP and there is no significant difference between the tricritical temperatures for systems with attractive and repulsive walls. As for m * = 1.5 ( Fig. 1(a)), also for m * = 2 and L/σ = 10 the estimates obtained from the Kelvin equation are in good agreement with the full phase equilibrium calculations. 1-3. In (a) the crossings of the shifted phase boundary for the attractive wall confinement with the bulk phase boundary is more clearly visible than for L/σ = 10 shown in Fig. 3(a). peared in the bulk phase diagram, this holds also for the confined system, independent of the character of the wall. For m * = 2, upon decreasing the wall separation from L/σ = 10 to L/σ = 4 the topology of the phase diagrams does not change (compare Figs. 3 and 4). Figure  4(a) shows that for repulsive walls the stronger confinement leads to a larger shift of the phase boundary to higher values of the chemical potentials. For attractive walls this shift remains much smaller even for this narrow slab with L/σ = 4. Different from the case m * = 1.5 (see Fig. 2(a)), for m * = 2 two crossings of the bulk coexistence curve with the phase boundary for attractive walls are clearly visible. Their interpretation in terms of the Kelvin equation is less valid than for Fig. 3(a) because due to the small value L/σ = 4 subdominant terms in the expression for ∆µ become relevant. Fig-ure 4(b) shows that enhancing the confinement causes a further shrinking of the two-phase coexistence region. The tricritical temperatures T T CP and densities ρ T CP decrease upon reducing the wall separation. This means that for smaller distances between the walls the ferromagnetic fluid -isotropic fluid phase coexistence occurs within a narrower range of the thermodynamic variables.
In a recent publication Trasca and Klapp 7 have studied inter alia the second-order phase transition of strongly coupled dipolar fluids confined to narrow slit pores. They performed Monte Carlo simulations using purely repulsive wall confinements and the dipolar soft sphere model for the particles. In confined systems, they have found that, at a given temperature, with the decrease of the wall separation the paramagnetic-ferromagnetic phase transition density (averaged across the slit pore) decreases relative to the corresponding bulk one. They found that the direction of this shift is inconsistent with their very simple mean-field theoretical predictions (see Eq. (3a) in Ref. 7 ). Analyzing the shifts of the second-order transition lines for the repulsive wall confinements (see Figs. 3(b) and 4(b)) we find that within the framework of our DFT at a given temperature the critical averaged density decreases upon the decrease of the wall separation, which is in qualitative aggrement with the simulation findings of Trasca and Klapp 7 . (There is no possibility for a quantitative comparison because these Monte Carlo simulations were carried out for a reduced dipole moment m * = 3; for such large values of m * the quantitative reliability of the present DFT would be reduced anyhow.)

B. Structural properties
In the following the structural properties of the coexisting ferromagnetic fluid and isotropic gas phases are discussed. Figure 5 shows the corresponding results for reduced dipole moment m * = 1.5 and wall separation L/σ = 10 at the reduced temperature T * = 1.15. At this low temperature the ferromagnetic fluid can be denoted as a liquid (see Fig. 1(b)). (As a hint we note that for a given temperature the chemical potential at two-phase coexistence is different for attractive and repulsive walls, see Fig 1(a).) Figure 5(a) shows that the confined ferromagnetic liquid phase is strongly structured both for attractive and repulsive walls. In the case of confinement by attractive walls the contact value ρ * w ≃ 19.8 of the reduced density (i.e., at |z| = (L − σ)/2) is very high (not displayed on the scale of the figure) which refers to a strong adsorption on the walls. The density profiles of the coexisting isotropic gases are displayed in Fig.  5(b). For the attractive wall case the contact value of the density is ρ * w ≃ 0.55 which indicates a visible spatial inhomogeneity in the vicinity of the walls. However, in the case of repulsive walls, at this low temperature the density of the vapor phase does not show any detectable inhomogenity. The behavior of the preferential orientation α 20 (z) shown in Fig. 5(c) tells that in the liquid phase the dipole axes are preferentially oriented parallel to the walls throughout the pore for both types of confinements. The ferromagnetic ordering of the particles in the liquid phase is confirmed by the variation of the order parameter α xy (z) (Eq. (33)) which can be seen in Fig.  5(e) which describes the net orientation of the dipoles parallel to the walls, keeping in mind that for the orientational order parameter we find α 10 (z) = 0 (M z = 0, see Eq. (32)) throughout of the pore. We note, that for all confined fluid phases studied here we found that α 10 (z) = 0, i.e., throughout the pores the z component of the magnetization M z is zero. Figure 5(f) displays the corresponding magnetization M xy (z) (Eq. (34)), which shows strong spatial inhomogeneities for both types of confinements induced by the structure of ρ(z). (For the attractive walls the contact value of the magnetization is M xy,w ≃ 25.7.) For the vapor phase we find α xy (z) = 0 for both attractive and repulsive confinements, which together with M z = 0 implies that the vapor phase is an isotropic phase. However, in Fig. 5(d) the order parameter α 20 (z) shows that close to the walls there is a slight ordering in the vapor phases as well. There the ordering is more pronounced for the confinement by attrac- tive walls than by repulsive walls. In the vapor phase, close to the attractive walls the dipole axes are preferentially ordered parallel to the walls (α 20 < 0), but the amplitude of the orientation is much smaller than in the coexisting liquid phase (Fig. 5(d)). Between this layer of parallel orientation and the orientationally disordered interior of the pore (i.e., α 20 (z) is close to zero), there is a layer with preferential orientation perpendicular to the wall (α 20 > 0, Fig. 5(d)). With respect to the orientational order, our DFT results are consistent with those of Gramzow and Klapp 3 but they show stronger spatial inhomogenities which is borne out by using a more sophisticated free energy functional for the hard sphere reference system. For m * = 1.5 and T * = 1.15, upon decreasing the wall separation from L/σ = 10 to L/σ = 4 the density and orientational distributions change significantly. Figure 6(a) shows that for both attractive and repulsive confinements the ferromagnetic liquid phases become more structured. Moreover, the confinement by attractive walls induces more ordered structures than by the corresponding repulsive ones. The density profiles for the coexisting isotropic vapor phase are displayed in Fig.  6(b), telling that only the attractive walls induce a spatial inhomogeneity in the density distribution of the isotropic vapor. Figure 6(c) shows that in the liquid phase the extent of orientation of the dipole axes parallel to the walls is comparable with that in the wider (L/σ = 10) pore but it varies more strongly across the pore. According to Fig. 6(d) the same conclusion holds for the structural properties of the confined vapor phase; since the vapor phase is isotropic this shows that the preferential ordering parallel and perpendicular to the walls does not hinge on the formation of ferromagnetic order. From Fig. 6(e) one infers that the in-plane ferromagnetic order parameter of the liquid phase for L/σ = 4 is comparable with that for L/σ = 10. However, the thinner pore exhibits a stronger magnetization and more pronounced spatial variations ( Fig. 6(f)). Figure 7 displays the structural properties for the wide pore L/σ = 10 at the higher temperature T * = 1.6 and for an increased dipole moment m * = 2. Compared with Fig. 5 (T * = 1.15, m * = 1.5) there are no qualitative differences. The quantitative differences concern the vapor density, which is lower due to the higher temperature ( Fig. 7(b)), and the magnetization of the liquid phase, which is higher due to the increased dipole strength (Fig.  7(f)). The occurrence of preferential ordering parallel to the walls in the vapor phase (α 20 (z) < 0, Fig. 7(d)) is in line with the findings of Gramzow and Klapp 3 ; however their approach does not render the preferential ordering perpendicular to the walls (α 20 (z) > 0) in the subsequent layer towards the interior of the pore. But this preferential orientation is about a factor of ten smaller (Fig.  7(d)) as compared with the case shown in Fig. 5(d).
For m * = 2 and T * = 1.6, decreasing the wall separation from L/σ = 10 to L/σ = 4 leads to a significant increase of both the liquid and the vapor densities (Figs. 8(a) and (b)) as well as to more pronounced density oscillation. The preferential orientional ordering is similar for pore. The magnetization is stronger in the more confined pore (Fig. 8(f)); this and the more pronounced spatial variation is caused by the corresponding behavior of the underlying number density ( Fig. 8(a)). The vapor phase does not acquire a ferromagnetic order even in the narrow pore, i.e., α xy (z) = 0 for the vapor phase.

VIII. SUMMARY
In the present study of the phase behavior and structural ordering of confined Stockmayer fluids (Eq. (1)) the following main results have been obtained. 1. We have applied an extension of the modified meanfield density functional theory to determine the phase behavior and the structure of Stockmayer fluids in slitlike pores formed by either attractive (Eq. (7)) or repulsive (Eq. (6)) walls at wall separations L/σ = 10 and L/σ = 4, for reduced dipole moments m * = 1.5 and m * = 2, and for a reduced wall strength ǫ * w = 1/2. This system exhibits three distinct fluid phases: isotropic vapor, isotropic liquid, and ferromagnetic fluid. 2. For the reduced dipole moment m * = 1.5 we have found isotropic liquid -isotropic vapor, ferromagnetic fluid -isotropic vapor, and ferromagnetic fluid -isotropic liquid first-order phase separations. The confinement by repulsive walls shifts the phase boundaries in the chemical potential -temperature plane towards higher chemical potentials while attractive walls lead to a shift to-wards lower chemical potentials relative to the corresponding bulk data. The decrease of the wall separation from L/σ = 10 to L/σ = 4 does not change the topology of the phase diagrams, but in this case the shift of the phase boundaries relative to the bulk ones are more pronounced (see Figs. 1 and 2). The increase of the dipole moment from m * = 1.5 to m * = 2 erases the isotropic liquid -isotropic vapor phase transitions for both the confined and the bulk systems. We have found that for this stronger dipole moment only the phase boundary of the confined system with repulsive walls is significantly shifted relative to the bulk one (see Figs. 2 and 3). The phase boundaries of the bulk and of the confined system with attractive walls do intersect without a significant shift of the phase boundaries of the confined system relative to the bulk one. 3. We have shown that in the temperature -chemical potential plane the Kelvin equation provides a good estimate of the confinement induced shift of the chemical potential for isotropic vapor -ferromagnetic fluid ( Fig.  1(a)) and isotropic fluid -ferromagnetic fluid ( Fig. 3(a)) phase coexistence. The multiple crossings of the isotropic fluid -ferromagnetic fluid phase boundaries in the temperature -chemical potential plane correspond to a nonmonotonic temperature dependence of the contact angle of a drop of ferromagnetic fluid on a single wall. 4. On the basis of our DFT approach we have found that the number density ρ(z) of the confined liquidlike phases (coexisting with the corresponding vaporlike phases) for both dipole moments (m * = 1.5 and m * = 2) and for both confinements (L/σ = 10 and L/σ = 4), with either attractive or repulsive walls, are strongly structured (see Figs. 5(a)-8(a)). The figures for the order parameter α 20 (z) (Eq. (14)) for preferential orientation of the dipole axes show that in the liquidlike phases the dipole axes are preferentially oriented parallel to the walls throughout the pore for both attractive and repulsive walls as well as for both sizes of confinements (see Figs. 5(c)-8(c)). The ferromagnetic ordering of the dipoles in the confined liquidlike phases is displayed by the variation of the ferromagnetic order parameter α xy (see Eq. (33) and Figs. 5(e)-8(e)) which, together with the result M z (z) ≡ 0, describes the net orientation of the dipoles parallel to the walls. The corresponding magnetizations M xy (z) = mρ(z)α xy (z) are displayed in Figs. 5(f)-8(f). Since α xy (z) = 0 and α 10 (z) = 0 throughout the pore for the vaporlike phases (at coexistence with the ferromagnetic liquidlike phases), they are completely isotropic. Attractive wall confinements give rise to strong adsorption at the walls even in the vaporlike phases (see Figs. 5(b)-8(b)). In the vaporlike phases, close to the attractive walls a slight preferential parallel ordering has been found. Between the layer of this preferential parallel orientation and the orientationally disordered interior of the pore a layer with preferential perpendicular orientation has been detected (see Figs. 5(d)-8(d)). One of the main limitations of the original version of fundamental-measure theory is that the underlying bulk fluid equation of state reduces to the Percus-Yevick compressibility equation. As it is well known, within this approximation the contact value of the density profile of a one-component hard-sphere fluid at a planar wall is significantly overestimated at high bulk densities. In order to improve the quantitative accuracy and to maintain internal consistency, following Roth et al. 17 and Tarazona 25 here we apply a fundamental-measure theory based on the Carnahan-Starling equation of state. Within this approach the hard-sphere excess free energy density functional is given in terms of weighted densities as where Φ = −n 0 ln(1 − n 3 ) + n 1 n 2 − n 1 ·n 2 1 − n 3 +(n 3 2 − 3n 2 n 2 ·n 2 ) In Ref. 26 it is shown that the weighted densities in the present slab geometry are given as where the reduced weight functions w (i) as functions of z are The chemical potential functional (first-order direct correlation function) of the reference hard sphere system is In the summation (Eq. (A46)) α = 4, 5 referes to the vectorial weight functions w (1) and w (2) , respectively.   (29)). By using Eq. (29) we have calculated these coefficients and in the following we provide all coefficients which are nonzero for i, l 1 , l 2 = 0, 1, 2: so that the functions are symmetric, too, which simplifies Eqs. (36) and (37). We note that the integral in Eq. (28) splits into two parts due to the Heaviside function Θ(R 2 + z 2 − σ 2 ): which fulfills the normalization requirement in Eq. (9). Using the expansion of α(z, ω) in terms of spherical harmonics (Eq. (13)) one obtains for the coefficients α lm (z) (Eq. (14)): Equations (C53) and (C55) are a set of coupled integral equations which are solved iteratively for the density profile ρ(z) and the coefficients α lm (z) of the orientational distribution. To this end, in Eq. (C53) the orientational entropy term must be expressed in terms of the coefficients α lm (z). Using Eq. (C54), after some calculations one arrives at dωα(z, ω) ln[4πα(z, ω)] = ln[4π/C(z)] − 2β l1,m1 l2,m2 α l1m1 (z) where C(z) = dω exp   −2β