Density profiles of Ar adsorbed in slits of CO_2: spontaneous symmetry breaking

A recently reported symmetry breaking of density profiles of fluid argon confined by two parallel solid walls of carbon dioxide is studied. The calculations are performed in the framework of a nonlocal density functional theory. It is shown that the existence of such asymmetrical solutions is restricted to a special choice for the adsorption potential, where the attraction of the solid-fluid interaction is reduced by the introduction of a hard-wall repulsion. The behavior as a function of the slit's width is also discussed. All the results are placed in the context of the current knowledge on this matter.


I. INTRODUCTION
In a quite recent paper, Berim and Ruckenstein [1] have reported symmetry breaking of the density profile of fluid argon (Ar) confined in a planar slit with identical walls of carbon dioxide (CO 2 ). These authors claimed that a completely symmetric integral equation provides an asymmetric profile which has a lower free energy than that of the lowest symmetric solution leading to a symmetry breaking phenomenon. It was assumed that the Ar atoms interact via a standard Lennard-Jones (LJ) potential characterized by the strength ε f f and the atomic diameter σ f f . The presented results were obtained from calculations carried out with the smoothed density approximation (SDA) version [2,3] of the nonlocal density functional theory (DFT) in the case of a closed planar slit with an effective width of 15 σ f f . The symmetry breaking was found for temperatures between the experimental triple point for Ar, T t = 83.8 K, and a critical value T sb = 106 K. At each temperature, it was determined a range of average densities ρ sb1 ≤ ρ av ≤ ρ sb2 where the symmetry breaking occurs, outside this range a symmetric profile has the lowest free energy.
As a matter of fact, the adsorption of fluid Ar on a solid substrate of CO 2 was intensively studied for several decades (see, e.g., Refs. 4 and 5). In 1977 Ebner and Saam [6] analyzed phase transitions by assuming that atoms of the fluid interact with the solid wall via a 9-3 van der Walls potential (from here on denoted as ES potential) obtained from the assumption that Ar atoms interact with CO 2 molecules via a LJ interaction with parameters ε sf and σ sf . After this pioneering work, a large amount of work has been devoted to study this system with different numerical and ana-lytical techniques. The attention was focused to analyze features like: the oscillatory behavior of the density profile which leads to a layered structure in the neighborhood of the flat substrate; the thin-to thickfilm transitions; wetting properties and prewetting jumps [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. Berim and Ruckenstein [1] have also adopted the ES potential, however, a hard-wall repulsion was introduced in their calculations. In practice, such a hard wall diminishes the strength of the solid-fluid interaction.
The investigation of symmetry breaking in physical systems is a very exciting issue. This is due to the fact that such a feature may have fundamental theoretical implications. In many fields of physics the discovery of symmetry breaking lead to significant advances in the theory. Therefore it is important to place the results of Ref. [1] in the context of the current knowledge about adsorption into planar slits.
Asymmetric solutions for fluids confined in slits have been previously reported in the literature. About 20 years ago a Dutch Collaboration has carried out calculations on the Delft Molecular Dynamics Processor (DMDP), which was specially designed for Molecular Dynamics (MD) simulations of simple fluids [23,24]. The results were published in a series of papers by Sikkenk et al. [25,26] and Nijmeijer et al. [27,28]. The simulations were performed for a canonical ensemble with two types of particles, 2904 of one type for building a solid substrate and several thousand of the other type for composing the fluid adsorbate. The temperature of the system was kept at T * = k B T /ε f f = 0.9 which is in between the fluid's triple-point temperature T * t ≃ 0.7 and the critical temperature T * c ≃ 1.3. The width of the slit was taken as L = 29.1 σ f f , supposing that this distance be enough large to avoid any capillary effect. Such a sys-tem can support solid-liquid (SL), solid-vapor(SV), and liquid-vapor (LV) interfaces. These authors have studied wetting at LV coexistence by varying over a wide range the relative strength of the solid-fluid and fluidfluid interactions defined by the ratio ε r = ε sf /ε f f of the LJ parameters. The length scale of this interaction was taken as σ sf = 0.941 σ f f . For increasing ε r from ε r ≃ 0.1 towards ε r ≃ 1.0 three cases were observed: (i) at low ε r , symmetric profiles consisting of two SV interfaces and two LV interfaces are obtained, this situation corresponds to a complete wall drying as can be seen in Fig. 1 of [28]; (ii) at intermediate ε r , asymmetric profiles consisting of a SL, a LV, and a SV interface are obtained, here the wall attraction is sufficiently strong to produce a partial wetting, i.e., to support a rather thick film on one wall while a SV interface is present near the other wall, this feature is shown in Fig. 2 of [28]; (iii) for the largest ε r , symmetric profiles consisting of two SL interfaces and two LV interfaces are obtained, now the strength is enough to wet both walls as can be seen in Fig. 3 of [28].
The structure of the profiles mentioned above depends on the balance of the involved surface tensions γ SL , γ LV , and γ SV which are related by the Young's law (see, e.g., Eq. (2.1) in Ref. [29]) where θ is the contact angle. The latter quantity is defined as the angle between the wall and the interface between the liquid and the vapor (see Fig. 1 in Ref. [29]). The transition from (i) to (ii) takes place at the drying point θ = π, whereas the transition from (ii) to (iii) takes place at the wetting point θ = 0. It is worth of notice that Velasco and Tarazona [30] have carried out calculations in the framework of the SDA obtaining density profiles with the same structure to that reported in Refs. [25,26,27,28]. The reader may look at Ref. [31] for a further comparison between MD and DFT results.
It is the aim of the present work to acquire a more accurate picture of the symmetry breaking reported in Ref. 1. In so doing, we explore size effects by comparing the results obtained for slits of widths 15 σ f f and 30 σ f f . Next, we investigate the existence of stable asymmetric solutions for the density profiles when the position of the hard-wall repulsion introduced in Ref.
1 is changed. When the location of this hard wall is moved, the strength of the adsorption potential is varied allowing a connection to the studies described in Refs. [25,26,27,28,30,31]. Several properties of the obtained solutions are discussed.
The paper is organized as follows. In Sec. II we provide a summary of the formalism underlying the present calculations. Special attention is devoted to the adsorption potential. The results and its analysis are reported in Sec. III. Final remarks are given in Sec. IV.

II. THE MODEL
The properties of a fluid adsorbed by an inert solid substrate may be studied by analyzing the grand free energy [32] where F is the Helmholtz free energy, µ the chemical potential, and N the number of particles of the adsorbate Quantity F contains the energy due to the interaction between fluid atoms as well as the energy provided by the confining potential. In a DFT it is expressed in terms of the density profile ρ(r) Here F int [ρ(r)] is the intrinsic Helmholtz free energy functional and U sf (r) is the external potential produced by the slit's walls. This formulation is usually applied to systems described by the grand canonical ensemble, i.e., at constant volume V , temperature T , and chemical potential µ. Such a situation corresponds to an open system in contact with reservoir which fixes T and µ. A minimization of Ω with respect to ρ(r) leads to the Euler-Lagrange equation for the density profile and the number of particles may be evaluated with Eq. (2.2). For a closed system, i.e., a canonical ensemble with fixed N , one should treat µ as an unknown Lagrange multiplier to be determined from the minimization procedure.

A. Density functional theory
Let us now summarize the DFT adopted for F int [ρ(r)]. In the case of inhomogeneous classical fluids at temperature T the intrinsic free energy functional is decomposed into two kind of contributions: (i) the ideal gas term F id [ρ(r)], which is given by the exact expression with Λ = 2 π 2 /m k B T being the thermal de Broglie wavelength of a molecule of mass m; (ii) the excess term F ex [ρ(r)], which accounts for the interparticle interactions is a unique but unknown functional of the local density. For fluids with attractive interactions as the Lennard-Jones (LJ) one, the free energy is decomposed into the repulsive and attractive contributions. The repulsive interactions are then approximated by a hard-sphere functional with a certain choice of the hard-sphere diameter d HS whereas the attractive interactions are treated in most cases in a mean field fashion (2.6) Here Φ attr (r =| r − r′ |) is the attractive part of the LJ potential.
In summary, the intrinsic Helmholtz free energy functional may be expressed as The free energy functional for hard spheres plays a central role in DFT. Expressions for f HS [ρ(r); d HS ] may be taken from the Percus-Yevick [33] or Carnahan-Starling (CS) [34] approximations for the equation of state (EOS) of a uniform non-attractive hard-sphere fluid (see, e.g. Ref. [35]). In a nonlocal DFT this quantity is evaluated as a function of a conveniently averaged densityρ(r).
For the calculations performed in the present work we used the same SDA formalism adopted in the paper of Berim and Ruckenstein [1]. In this approach developed by Tarazona [2,3], the excess of free energy density of hard spheres is written according to the semi-empirical quasi-exact CS expression Hereη =ρ(r) V HS is the packing fraction, where the factor V HS = π d 3 HS /6 is the volume of a hard sphere. The smoothed densityρ(r) is defined as with the following weighting function: The expansion coefficients w 0 (r), w 1 (r), and w 2 (r) are density independent and its expressions as a function of r =| r − r′ | are given in the Appendix of Ref. 3.
To account for the fluid-fluid interaction we adopted, as in Ref. [1], the spherically symmetric L-J potential given in Eq.
where σ f f is the hard-core diameter of the fluid. The authors of Ref. 8 have used this L-J version just for studying the adsorption of Ar on CO 2 and it has been also utilized in several subsequent works on this system. The values of the interaction parameters for Ar are ε f f /k B = 119.76 K and σ f f = 3.405Å.

B. The Euler-Lagrange equation
The equilibrium density profile ρ(r) of the fluid adsorbed in a closed slit is determined by a minimization of the free energy with respect to density variations with the constraint of a fixed number of particles N Here, i.e., for an ensemble with fixed V , T , and N , the Lagrange multiplier µ is an unknown quantity which should be determined from the constraint. It plays a role of a chemical potential but off the liquid-vapor coexistence conditions. Hence, it is not necessarily equal to µ coex of an open slit in equilibrium with a reservoir at temperature T (see, e.g., Ref. 36).
In the case of a planar symmetry where the flat walls exhibit an infinite extent in the x and y directions the profile depends only of the coordinate z perpendicular to the substrate. For this geometry the variation of Eq. (2.12) yields the following Euler-Lagrange (EL) equation For a slit of effective width ℓ w this EL equation may be cast into the form where The number of particles per unit area of one wall of the slit is In order to get solutions for ρ(z) it is useful to rewrite Eq. (2.14) as The relation between µ and N s is obtained by substituting Eq. (2.17) into the constraint given by Eq. (2.16) For the calculations carried out in the present work we set d HS = σ f f as it was done in Ref. 1.
The asymmetry of the density profiles is measured by the parameter (2.20) According to this definition, if the profile is completely asymmetrical about the middle of the slit [ρ(z < ℓ w /2) = 0 and ρ(z ≥ ℓ w /2) = 0] this parameter becomes unity, while for symmetric solutions it vanishes.

C. Adsorption potential
The model van der Waals (9-3) potential proposed by Ebner and Saam [6], i.e. the ES potential, is with ǫ ef f = ε sf ρ s σ 3 sf being the effective strength, was adopted for almost all the abovementioned studies of the adsorption of Ar atoms on a flat wall of solid CO 2 . The exception is the experimental and theoretical investigation performed by Mistura et al. [22], where a more realistic adsorption potential calculated on the basis an ab initio expansion of Marshall et al. [37] was used. The ES expression is obtained when one assumes that Ar atoms interact with CO 2 atoms via a Lennard-Jones (12-6) potential and subsequently integrates this potential over a continuum of CO 2 substrate atoms with a reduced density ρ * s = ρ s σ 3 sf = 0.988. The cross-parameters of the potential are determined by using the Lorentz-Berthelot rules. So that, the van der Waals strength ε sf is the square root of the product of the argon and CO 2 van der Waals strengths, while the hard-core diameter σ sf is the mean of the argon and CO 2 hard-core diameters, while The parameters evaluated in this way are ε sf /k B = 153 K and σ sf = 3.727Å.
Berim and Ruckenstein [1] have investigated the Ar-CO 2 system utilizing, in principle, the ES potential. However, by looking at their paper one realizes that according to Eq. (A5) of the Appendix which accounts for the solid-fluid interaction at one of the walls, a hard-wall repulsion was located at a distance σ sf from the real wall of the slit. In agreement with this assumption, the total confining potential exerted on Ar atoms by the two walls separated by a distance L was expressed as Here the effective width of the slit is This scenario is depicted in Fig. 2 of Ref. 1. In this context, it is interesting to notice that Nilson and Griffiths [38] in order to study the adsorption of a fluid in a planar slit have written in their Eq. (10) the total fluid-solid potential as i.e., locating a hard-wall repulsion at a distance σ sf /2 from the substrate. In this case, the effective width is as it is shown in Fig. 2 of Ref. 38.
In Fig. 1 we compare the potentials outlined in the previous paragraph. The comparison is restricted to the region close to the substrate. The quantity ξ is the perpendicular distance from the real wall being   potential before the minimum be reached. This feature produces important effects on the behavior of the density profiles. In fact, the calculations performed by Berim and Ruckenstein [1] yielded density profiles with ρ(z = 0) and ρ(z = ℓ w ) different from zero indicating that the fluid is in contact with the hard walls, while in the case of Nilson and Griffiths [38] the fluid forms a well defined first layer separated from the wall.
In the present work we shall analyze the evolution of asymmetric solutions when the total adsorption potential is written as

III. NUMERICAL RESULTS AND DISCUSSION
Let us now describe the obtained results. The EL equation (2.14) was solved at fixed ℓ w and T for a given number of particles per unit area N s . The latter quantity determines an average fluid density ρ av = N s /ℓ w . A widely used computational algorithm consisting of a numerical iteration of the coupled Eqs. (2.17)-(2.19) was applied. This procedure yields the density profile ρ(z) and the value of the Lagrange multiplier µ. The convergence of the solutions are measured by the difference between two consecutive profiles where ρ i (z) is the density profile after the i-th iteration, and by the quantity accounting for the deviation from the required N s . In practice, for the calculations it is convenient to use dimensionless variables: z * = z/σ f f for the distance, ρ * = ρ σ 3 f f for the densities, and T * = k B T /ε f f for the temperature. In these units the average density becomes For the numerical task, the region of integration [0, ℓ * w ] was divided into a grid of equal intervals ∆z * = 0.02, i.e., a grid with 50 points per atomic diameter σ f f . It is worthwhile to notice that in the work of Berim and Ruckenstein the number of grid points was taken equal to 10 per atomic diameters. If the obtained profile did not change with increasing precision from δ 1 ≈ 10 −8 to δ 1 ≈ 10 −15 , then it was accepted as a solution of the coupled integral equations.
In a first step, we studied the same systems treated in detail by Berim and Ruckenstein [1]. Hence, we set ν = 1 and solved the EL equation for a slit with an effective width ℓ * w = 15 at T = 87 K (T * = 0.73) for a series of average fluid density ρ * av = N * s /ℓ * w . The ground state solutions yield symmetric density profiles for ρ * av < 0.1 = ρ * sb1 and for ρ * av > 0.514 = ρ * sb2 , while in the range ρ * sb1 ≤ ρ * av ≤ ρ * sb2 the ground-state solutions provide asymmetric density profiles. This is due to the fact that in such a regime the asymmetric solutions have lower free energy than the symmetric ones. The free energies calculated for some selected ρ * av are listed in Table I together with the results obtained by Berim and Ruckenstein [1]. A glance at this table indicates a good agreement between both sets of values. In order to get symmetric solutions in the range 0.1 ≤ ρ * av ≤ 0.513 one must explicitly impose such a condition to Eqs. (2.17)- (2.19).
The Lagrange multiplier µ (equivalent to the chemical potential in the case of open slits) is displayed in Fig. 2 as a function of average density. We show the results for a wider range of ρ * av than it is done in Fig. 9 of Ref. 1. Figure 2 clearly indicates that the asymmetric solutions  occur in the domain where µ is degenerate, namely, the same value of µ corresponds to different ρ * av . That is just the regime where in the case of an open slit the equalarea Maxwell construction should be applied in order to determine the chemical potential (cf., e.g., Fig. 2 in Ref. 36). In addition, it is worthwhile to notice that at ρ * sb2 there is an abrupt jump of µ, we shall come back to this feature below.
If one manages conveniently the EL equation it is possible to get another kind of symmetric solutions in some region of ρ * av . The free energy of such new solutions is included in Table I and the corresponding multiplier µ is plotted in Fig. 2. Figure 3 shows the spacial distribu- tion of all three states listed in Table I for ρ * av = 0.3865. A direct comparison indicates that these three profiles correspond to the drying , the one-wall -wetting, and the two-walls-wetting cases displayed, respectively, in Figs. 1-3 of Ref. [28]. From a glance at Table I one realizes that as long as the film wetting solutions exist the drying one is a symmetric excited state. In the regime ρ * av > ρ * sb2 the drying profile shown in Fig. 3 becomes the capillary condensation (CC) solution. Figure 4 shows a series of asymmetric density profiles ρ * (z) = ρ(z) σ 3 f f . This sequence indicates that for increasing average density, starting from the profile denoted as 1 (ρ * av = 0.1063), the number of oscillations near the left wall as well as its amplitudes increase. This trend continues up to the profile 2 (ρ * av = 0.4638), where the peaks of the oscillating structure begin to decrease. Furthermore, profile 3 corresponds to the biggest asymmetric solution, for larger ρ * av there is a jump to the symmetric profile 4. The latter behavior indicates a transition to the so-called CC phase, i.e., a transition to a situation where the slit is full of liquid argon. The jump of µ addressed in the previous paragraph is also a manifestation of this thick-film to CC transition.
The asymmetry of the profiles for ν = 1 displayed in Fig. 4 was measured by the parameter ∆ N introduced in Eq. (2.20). The results are shown in Fig. 5. These values are essentially the same as that of the equivalent curve displayed in Fig. 5 of Ref. 1. As mentioned in the introduction, in Refs. [25,26,27,28,30,31] it is emphasized that a slit of width ℓ * w = 30 is appropriate to study wetting because it is large enough to avoid any confinement effect. Therefore, it becomes of interest to solve the EL equation for such a big slit and to compare the results with that obtained for ℓ * w = 15. The asymmetry parameter evaluated for the ℓ * w = 30 slit at T = 87 K is also plotted in Fig. 5  For the sake of completeness we also analyzed slits of ℓ * w < 15 at T = 87 K. The asymmetry parameters evaluated for ℓ * w = 12, 7.5 and 6 are included in Fig. 5. It is clear that these curves exhibit everywhere an increasing departure from the ℓ * w = 30 results. Finally, for ℓ * w = 5.5 the parameter ∆ N becomes zero for all ρ * av indicating that the asymmetric solution disappears.
In the next step, we analyzed the change of the properties described above when the parameter ν is taken larger than unity. The variation of ν is performed in such a way that the effective width of the slit is always kept at ℓ * w = 15, hence, only the potential is slightly changed mainly near the walls. At this point the temperature was still kept at T = 87 K. It was found that for increasing values of ν the range of average density ρ * sb1 ≤ ρ * av ≤ ρ * sb2 where there are symmetry breaking decreases. This effect might be observed in successive plots of µ vs. ρ * av , however, we prefer to report directly the evolution of the parameter ∆ N . Figure 6 shows how this parameter decreases with increasing ν. From this figure one may conclude that symmetry breaking persists at most for ρ * av ≈ 0.2. Figure 7 shows the asymmetry parameter ∆ N as a function of the factor ν for the average density ρ * av = 0.1932. These data indicate that the asymmetric solution already disappears for a critical value ν c ≃ 1.18. The evolution of the density profiles from asymmetric to symmetric species is displayed if Fig. 8. In this drawing one may observe the diffusion of fluid argon from the neighborhood of the left wall towards the right one. This process continues until a symmetric density profile is formed for ν ≃ 1.2. The described evolution of ρ * (z) is determined by the behavior of the surface tensions. If one uses the Young's law given in Eq. (1.1) the total surface excess energy of asymmetric profiles may be written as with cos θ = (γ SV −γ SL )/γ LV < 1. By increasing enough the attraction the equality γ SV − γ SL = γ LV is reached yielding cos θ = 1, then the system undergoes to a transition to a symmetric profile with γ tot = 2 γ SL + 2 γ LV . (3.4) In this case both walls of the slit are wet.
Let us now look if something special is going on for the adsorption potential at ν c ≃ 1.18. As a matter of fact, it becomes important to explore where the hardwall repulsion introduced in Ref. 1 is located for such a critical value. The position of the left hard-wall repulsion with respect to the real wall according to Fig. 1 is given by For a slit of effective width ℓ * w = 15 the minimum of the total adsorption potential given by Eq. (2.25) is to a very good approximation determined by the minimum of the ES potential exerted by the left wall. It is located at The ratio of these quantities becomes ξ min ξ c Since this ratio is larger than unity, it indicates that for ν c the minimum of the ES potential is reached inside the slit.
We also analyzed the occurrence of asymmetric solutions for 1 ≤ ν < 2 at temperatures T > 87 K. It was also found that the symmetry breaking disappears at some critical value of ν. Such a behavior may be expected if one takes into account that for temperatures larger than T = 87 K even for ν = 1 the range ρ sb1 ≤ ρ av ≤ ρ sb2 is smaller (see Fig. 5

of Ref. 1).
It is worthwhile to notice that a reliable ab-initio potential utilized by Mistura et al. [22] for investigating the adsorption of Ar on CO 2 exhibits an even stronger attraction than the ES potential. This feature can be observed in Fig. 3 of [22]. Hence, one would not expect any symmetry breaking in the case of such a realistic adsorption potential.

IV. FINAL REMARKS
We reexamined the symmetry breaking found very recently by Berim and Ruckenstein [1] in a study of the adsorption of argon in a closed slit with identical walls of carbon dioxide. It is important to stress that these authors have introduced hard-wall repulsions at distances d BR = σ sf from the real walls of the slit as shown in Fig.  1 (see also Fig. 2 in Ref. [1]), reducing in such a way the strength of the adopted ES potential. Stable asymmetric solutions were obtained in the case of a slit with effective width ℓ * w = 15 for temperatures in the range 87 ≤ T ≤ 106 K.
The present study was mainly devoted to establish how robust are the asymmetric solutions against changes of the adsorption potential. However, in addition, the behavior of the asymmetry parameter was also evaluated for slits of different widths. The calculations have been carried out using the same nonlocal formalism as that adopted in Ref. 1, namely, the SDA density functional theory proposed by Tarazona [2].
Since in a pioneering series of works a Dutch Collaboration have previously found the symmetry breaking by studying adsorption in a slit of width ℓ * w = 30, see Refs. [25,26,27,28], we performed a comparison of the asymmetry parameter ∆ N evaluated for slits of ℓ * w = 15 and 30 at T = 87 K. The difference between the results is rather small. Further calculations showed that the effect due to confinement begins to be important for slits with ℓ * w ≤ 12.
Focusing the analysis on the slit of ℓ * w = 15, in a rather complete plot of µ vs. ρ * av we show clearly the regime where the symmetry breaking occurs and, in addition, we also display the Lagrange multiplier µ for a CC like solution. Furthermore, it is shown that the spacial shape of the obtained solutions correspond to the three sorts of profiles discussed in Refs. [25,26,27,28]. The strength of the attraction determines which one of that profiles is the stable solution.
By shifting the hard-wall repulsions introduced by Berim and Ruckenstein [1] towards the real walls the strength of the attraction is increased producing changes in the balance of surface tensions. It was found that at T = 87 K already for the critical distance d c ≃ d BR /1.2 the asymmetric solutions disappear. At this temperature, close to the triple point, the asymmetry parameter ∆ N reached its largest values in the study of Ref. 1. For higher temperatures the symmetry breaking disappears even more rapidly.
Finally we can state that from the present study it is possible to conclude that the symmetry breaking reported in Ref. 1 can be understood in terms of findings described in Refs. [25,26,27,28]. By locating hard-wall repulsions the authors of Ref. 1 diminish the attraction of the ES adsorption potential exerted on the fluid causing the entrance of the system in a one-wall -wetting regime (asymmetric profiles) corresponding to Fig. 2 of Ref. [28]. By moving the hard-wall repulsions towards the real walls the attraction increases monotonically causing eventually the entrance of the system in the two-walls-wetting regime with symmetric solution of the type displayed in Fig. 3 of Ref. [28].
Furthermore, a reliable ab-initio potential for the interaction of fluid Ar with a structureless smooth wall of CO 2 like that adopted by Mistura et al. [22] is even stronger than the ES one. Therefore, it is possible to infer that no symmetry breaking would be expeted for real Ar/CO 2 systems.