Confined Quantum Hard Spheres

We present computer simulation and theoretical results for a system of N Quantum Hard Spheres (QHS) particles of diameter σ and mass m at temperature T, confined between parallel hard walls separated by a distance Hσ, within the range 1≤H≤∞. Semiclassical Monte Carlo computer simulations were performed adapted to a confined space, considering effects in terms of the density of particles ρ*=N/V, where V is the accessible volume, the inverse length H−1 and the de Broglie’s thermal wavelength λB=h/2πmkT, where k and h are the Boltzmann’s and Planck’s constants, respectively. For the case of extreme and maximum confinement, 0.5<H−1<1 and H−1=1, respectively, analytical results can be given based on an extension for quantum systems of the Helmholtz free energies for the corresponding classical systems.


Introduction
The behavior of confined fluids have received considerable interest in recent years due to the research on microfluidic applications. Phase transitions that occur in pores are different to the ones observed in bulk systems, as has been determined in several studies of water in nanochanels [1]. Hydrogen in nanopores exhibits important quantum effects even at room temperatures [2]. The study of confined models is of great interest due to their multiple applications. Previous studies have found that the behaviour of a classical system of Lennard-Jones (LJ) particles confined between parallel walls separated by a distance Hσ, where σ is the size parameter of the LJ pair potential u(r) where u(σ) = 0, presents a reduction in their critical point, with a clear asymptotic tendency to a two-dimensional LJ fluid when the separation between walls is reduced to maximum confinement [3]. In the case of the liquid-solid transition, the thermodynamic melting temperature of LJ system confined between parallel walls is higher than for the bulk case (H −1 = 0) [4]. Schmidt and Löwen [5] studied the phase diagram of a hard-spheres (HS) system confined between parallel hard walls, and found a very rich formation of novel phases.
The case of quantum systems confined in slit pores has been studied previously for the HS model by Liu et al. [6]. In previous work [7] we have addressed the behavior of quantum fluids under confinement, considering quantum LJ particles (QLJ) confined between parallel hard walls, and introducing quantum corrections at the level of the Wigner-Kirkwood theory [8,9], with interesting findings in the way that the reduction of dimensionality from bulk (3D) to extreme confinement (2D) gives rise to the coupling of criticality with the pore size and de Broglie's thermal wavelength λ B = h/ √ 2πmkT, where k and h are the Boltzmann's and Planck's constants, respectively, and m is the mass of a particle. From the fundamental solution of the Schrödinger equation, the behaviour of confined D-dimensional hard spheres has been determined for low-density systems [10] based on the assumption that diluted interacting quantum gases can be described as an ideal quantum gas interacting with small excluded spheres randomly distributed, i.e., a kind of an effective Bernoulli model of an ideal gas. A closed related-approach has also been used for classical hard particles [11,12], providing a route to construct equations of state for classical fluids [13]. Effects of confinement in quantum gases have also been studied in detail previously [10,14,15], where contributions to the free energy are functions of the size of confinement, the geometry and the connectivity or number of holes inside the container. What is clear is that the thermodynamics of dense quantum fluids become more elaborate with a complex dependency on the geometry and nature of the container, the thermal wavelength λ B , density, temperature and spin.
The standard way for dealing with the prediction of structural and thermodynamic properties of fluids contained within pores is using density functional theory (DFT). In this paper we follow a different approach, with the motivation of introducing thermodynamic perturbation theory for confined fluids. We present Monte Carlo computer simulation results in the Canonical (NVT) and Isothermal-Isobaric (NPT) ensembles of quantum hard hard spheres (QHS) confined between hard planar walls separated by a distance Hσ, based on a semiclassical approach developed by Yoon and Scheraga for bulk fluids [16] and that we have adapted to confined systems. In the cases of maximum (H −1 = 1) and extreme (0.5 < H −1 < 1) confinement we obtain analytical expressions for the Helmholtz free energy based on the approach given by Franosch et al. [17].

Method
If we consider a system of N quantum hard-sphere particles whose potential energy is U HS (r N ), then the partition function neglecting exchange effects is given by where W(r N ) is the Slater sum defined as [18] W(r N ) = ∑ n e −βE n Ψ * n Ψ n, where E n is the energy for the system's quantum states whose stationary wave function is Ψ n , according to the Schrödinger equation,Ĥ, whereĤ is the Hamiltonian operator. In the classical limit, whereas in the semiclassical approximation λ B < σ and the Slater sum of the system can be expressed as [19,20] where W q (r N ) gives the quantum correction to the classical value. The Slater sum W q (r N ) can be be expressed as a Zwanzig expansion [21] in terms of classical ensemble averages of the quantum Ursell functions U n , considering now the quantum version of these functions, obtained through the solution of the Schrödinger equation. In the two-particle cluster approximation, where with where λ * B = λ B /σ and r * ij =| r i − r j | /σ is the reduced distance between particles i and j with positions r i and r j , respectively. The two-particle Slater function is then given by Using Equations (1), (4) and (5) we can derive a semiclassical partition function for the QHS system, that follows the same mathematical structure as the Zwanzig perturbation theory for classical fluids [22], and introducing the classical HS partition function, where < ... > HS denotes an ensemble average with respect to the classical HS system. The corresponding Helmholtz free energy A QHS = −kT ln Z QHS is then obtained, using Equations (9) and (11), Last expression can be interpreted as the perturbation expansion for a classical system of particles interacting with the pair potential u HS (i, j) − kTU 2 (i, j), where u HS (i, j) is the HS binary potential for particles i and j and −kTU 2 (i, j) < 1 is the perturbation potential. Additionally, where r is the relative position between two particles. Consequently, Equation (12) can be rewritten as that corresponds to the first-order perturbation expansion of the Helmholtz free energy for a QHS system, in terms of its classical counterpart as reference system, with a mean-energy term for the potential −kTU 2 (r). An alternative expression for A QHS was developed by Nordholm and coworkers [23,24], where quantum effects were included in the partition function of N independent HS particles, where packing behaviour was modelled as if one particle is enclosed in a cage induced by the remaining N − 1 particles. This excluded a volume mechanism can be described quantum mechanically using the standard particle-ina-box solution of wave mechanics.
In the MC simulation scheme given by Yoon and Scheraga [16], Equation (14) is used rewritten in terms of the two-body Slater sum w(r), and the approximation U 2 (r) ≈ ln(1 + U 2 (r)) = ln w(r), i.e., Introducing an effective Boltzmann factor, the semiclassical MC simulations are performed with the standard classical MC method, using Γ ij as an equivalent energy and the Slater sum as the quantum-mechanical probability instead of the Boltzmann one. Quantum hard spheres follow Boltzmann statistics since exchange effects decay exponentially to zero in harsh repulsive systems at temperatures where de Broglie thermal wavelengths are comparable to the size of particles, and then quantum diffraction effects are not negligible. In the case of solid phases, previous studies have demonstrated that thermodynamic and structural properties of quantum systems are basically unaltered by modifying the type of statistics (Bose-Einstein, Fermi-Dirac or Boltzmann-Gibbs) [25]. Chandler and Wolynes [26] have demonstrated by a path-integral approach that exchange effects can be described as associating ring molecules, using the exact isomorphism between quantum theory and classical statistical mechanics of polyatomic ring fluids. In a first approximation, considering only dimer association at low temperatures and liquid densities, the hard-core repulsion of a QHS system with λ B < σ reduces the association constant, with respect to its ideal-gas value, by a factor of the order of τ ≈ 2e −2πσ 2 /λ 2 B . In this work we are considering systems with λ B ≤ 0.6σ, that correspond to τ ≤ 2 × 10 −8 , i.e., the majority of the ring molecules will be non-associated since the population of dimers is extremely low, an indication that in the semiclassical approach that we are following here we can safely neglect exchange effects even at high densities.

Results
Computer simulations were obtained using the method described in the previous section, applying the standard Metropolis algorithm [27] with a random sampling technique, for QHS particles confined between parallel hard walls, and considering the particlewall interaction where z is the distance between the fluid particles and the walls, and usual periodic boundary conditions were applied in the x and y coordinates. Simulations were performed using 512 QHS particles, for densities ρ * = 0.1, 0.

Confined QHS Global Behaviour
When we consider the confinement effect due to parallel hard walls, the Helmholtz free energy of the system can be expressed as where A CHS is the Helmholtz free energy of confined classical hard spheres and A CQ is the contribution due to the quantum nature of the system. In the classical approximation, i.e., λ * B = 0, A CQ = 0. On the other hand, for the bulk case, i.e., H −1 = 0, A CHS (ρ * , H −1 = 0) = A HS (ρ * ), and Equation (18) reduces to the Helmholtz free energy of the bulk system, where A HS (ρ * ) can be obtained from the Carnahan-Starling equation of state [28], where A 3D ideal is the ideal contribution and where p 1 = 0.588868 and p 2 = 0.524248 are parameters obtained by fitting to the NVT-MC values. For the general case when H −1 = 0 we can determine A CQ from NVT-MC simulations, as can be seen in Figure 1. We observe that, for all thermal wavelengths, two different ranges are observed where A CQ behaves differently. When 0 ≤ H −1 < 0.4, the free energy contribution due to quantum effects and confinement varies slowly with respect to the bulk phase (H −1 = 0). When 0.4 < H −1 ≤ 1, A CQ has strong variations with H −1 for densities ρ * > 0.4. It is convenient to discriminate a region for 0.5 < H −1 < 1, that will be characterised as of extreme confinement that ends with the case of maximum confinement, H −1 = 1, when the distance between walls is equal to a HS diameter. In this limit, the behaviour of the system can be approximated by a hard disks system. Following the same approach used for bulk QHS, the Helmholtz free energy for quantum hard disks (QHD) is obtained from the equation of state derived by Henderson for classical hard disks [29] as a function of a 2D effective packing fraction γ e , in order to reproduce the NVT-MC simulated values, where and q 1 = 0.35367 and q 2 = 0.224205. For this equation we are considering that ρ 2D = N/A s is a 2D density for N particles distributed along a surface with area A s , and γ is the actual 2D packing fraction of the system,γ = πρ 2D σ 2 /4. The theoretical prediction for the compressibility factor for QHD systems, Z QHD = βP/ρ, is presented in Figure 2, compared with NPT-MC simulation values, obtained from a standard NPT algorithm and also applying the Test Area MC method [30]. From these results it is clear that the introduction of quantum effects increases the values of the pressure with respect to the classical case, and the liquid-solid transition occurs at lower densities when λ B increases. Since it is well known that the quantum correction to hard particles increases its effective size, we can expect this behaviour. However, the nature of the transition as well the presence of a stable hexactic phase requires a more detailed study. The behaviour of the structure factor of QHS at maximum confinement is presented in Figure 3. The transition from the liquid to the solid phases is given for the same density, ρ * = 0.8, varying λ * B . The system evolves from a liquid phase (λ * B = 0) up to a solid phase (λ * B = 0.6) and in between a hexactic-like phase can be observed for λ * B = 0.2. This feature is consistent with the phase diagram for QHD given in Figure 1, since clearly for λ * B ≈ 0.1 the liquid-solid transition boundary has shifted to densities equal to 0.8.
The clear discontinuity observed for A CQ that arises around H −1 = 0.4 at high densities is consistent with the observation of a corresponding strong freezing transition in the phase diagram for confined classical HS around H −1 = 0.5 [5]. By the same effect observed in Figure 2 for hard disks, since the particles swell by increasing the thermal wavelength λ * B we can expect that the corresponding transitions will appear at higher values of H than in the classical state.

Equation of State for Extreme Confinement
For the case of extreme confinement, we follow the same procedure used for QHS and QHD, i.e., to consider an analytical expression for classical HS and then to map onto a quantum expression by considering an effective density.
Franosch et al. [17] derived an expression for the Helmholtz free energy of classical HS confined between two parallel hard walls, for the case 1 ≤ H ≤ 2, defining the confinement parameter L = H − 1, where A 3D ideal is the 3D ideal term, A HD is the free energy for hard disks of diameter σ L = σ 2 − L 2 1/2 , and ∆A is the contribution arising from confinement, given by where V e f f (r; L) = −2kT ln 1 − (σ 2 − r 2 )/L 2 and g(r) are the radial pair distribution function of the hard-disk reference fluid of diameter σ L . The integral involved in Equation (25) can be performed using the approximation g(r) ≈ g(σ L ) [17], and then obtaining The classical expression can be used to describe the contribution to the free energy of the CQHS system, where now the free energy is given in terms of the corresponding quantum expressions, and the free energy of the quantum hard disks of diameter σ L is given by the Henderson equation of state evaluated with an effective packing fraction that reproduces the NVT-MC simulation results for the confined system, where q 3 = (q 31 + q 32 γ)L + q 33 γ 2 L 2 , with q 31 = −0.40144, q 32 = 4.7426 and q 33 = −7.95290. Notice that in the limit L = 0 we recover the expression given by Equations (22) and (23). In Figures 4 and 5 results are given for A * CQHS as functions of density, L and λ B , as obtained from MC simulations and using the quantum version of the free energy of the confined system, Equations (27) and (28).

Discussion
In this work we have presented theoretical and computer simulation results for QHS particles confined between planar hard walls. Since the quantum nature of the system given by the thermal wavelength λ * B introduces a new parameter in the statistical and thermodynamic properties of the confined system, its coupling with the confining parameter H increases the complexity of the phase behaviour already observed for confined classical hard spheres [5]. One of the main features observed in our study is the swelling effect of the hard spheres that modifies the phase diagram of the confined classical system, particularly relevant in the case of maximum confinement, where the liquid-solid transition of QHD occurs at lower densities when compared with classical hard disks. We also considered analytical expressions for the case of extreme confinement, 0.5 < H −1 ≤ 1, based on expressions for the classical system. Similar to what has been already reported for the bulk system [16,31,32], these classical expressions are valid if the packing fraction is modified by considering the swelling effect of the quantum particles and the confinement parameter H. In a future communication we will explore this thermodynamic approach in order to obtain equations of state for a fluid confined in slit pores, based on perturbation theory, following the methods already developed for bulk systems [31,32] .

Data Availability Statement:
The data that support the findings of this article are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.