Mean properties and Free Energy of a few hard spheres confined in a spherical cavity

We use analytical calculations and event-driven molecular dynamics simulations to study a small number of hard sphere particles in a spherical cavity. The cavity is taken also as the thermal bath so that the system thermalizes by collisions with the wall. In that way, these systems of two, three and four particles, are considered in the canonical ensemble. We characterize various mean and thermal properties for a wide range of number densities. We study the density profiles, the components of the local pressure tensor, the interface tension, and the adsorption at the wall. This spans from the ideal gas limit at low densities to the high-packing limit in which there are significant regions of the cavity for which the particles have no access, due the conjunction of excluded volume and confinement. The contact density and the pressure on the wall are obtained by simulations and compared to exact analytical results. We also obtain the excess free energy for N=4, by using a simulated-assisted approach in which we combine simulation results with the knowledge of the exact partition function for two and three particles in a spherical cavity.


I. INTRODUCTION
The properties of a few particles confined in a spherical cavity has become attractive recently, due to the possibility of synthesizing the so-called silica nano-rattles [1] and their promising wide range of applications. An important biomedical application of these mesoporous silica nanoparticles is in the area of drug delivery, targeted to disease diagnosis and therapy [2]. Within the field of core-shell systems, hollow silica nano-particles containing smaller particles or molecules were synthesized in different conditions and packed as mesoporous materials [2,3]. Single silica nano-containers confining gold nano-particles have been proposed as an option for drug delivery vehicles in cancer therapies [4] and also as imaging probes with fluorescent dies [5]. Hollow nanospheres with mesoporous shells and metal cores have also interesting optical and electrical properties, which are relevant in imaging and catalysis applications. Hindering the interaction between metal nanoparticles, by encapsulating them by mesoporous shells, keep the metallic particles active for longer times. These systems are known as yolk-shell nanoparticles and are also promising as nanoreactors, in bio-medicine applications and photo-catalysis [1,6].
From a theoretical point of view, a few particles confined in a cavity can be regarded as the few-bodies limit of the liquid or fluid states and also as a limiting case in which the curvature of the confining wall is very significant. Systems of many particles interacting through simple potentials, that form bulk fluid phases, confined in spherical cavities produce a rich set of structural and thermal properties and were studied thoroughly [7][8][9].
The spherically confined HS system has been studied using molecular dynamics (MD) [16], Monte Carlo simulations (MC) and density functional theories (DFT) [17,18], mainly in the context of liquid state theory. However, the limit of few particles confined in a spherical cavity has received much less attention. One of us performed analytical calculations for two and three HS confined in a spherical cavity [19][20][21][22] and in pores of different shapes [21]. Huang and co-workers studied recently, by molecular dynamics simulation the structure and equation of state of hard sphere confined in a small spherical geometry, with a focus in a freezing-like transition at high densities [9]. Based on the HS interaction, Winkler et al. analyzed phase transitions of a colloid-polymer mixture inside a hard spherical cavity [12]. New insights on the curvature dependence of the surface free energy of a fluid in contact with a curved wall were also revealed, by studying the HS fluid [23,24].
In this work, we study the free energy and the structural and mean properties of few-body HS confined in a spherical cavity, as a function of the number of particles. For the sake of clarity we give operational expressions for configuration integral Z 4 and the pressure on the wall P w as a function of the effective cavity radius. Two different approaches are used: we perform event-driven molecular dynamics [25] with a fixed-temperature spherical wall and compare the results with exact or quasi-exact theoretical calculations, obtained from analytical expressions of the free energy. Our minimal model utilizes the simple HS interaction potential between particles and a spherical hard boundary, chosen to capture the purely geometric effect of the spherical confinement and the behavior of the translational entropy of the particles, for a wide range of number densities. The number of particles N and T are kept fixed. A typical configuration of the system with N = 4 particles, as obtained from computer simulations, and the relevant length scales of the system are shown in Ê Figure 1. (color online) A configuration of four hard spheres (yellow) confined in a spherical cavity with a number density ρ = 0.5σ −3 , as obtained from an event-driven molecular dynamics simulation. The sphere of radius Ro = 1.74σ, is shown with transparency. The diameter σ of the HS particle is indicated in one of them, and the effective radius R eff ≡ Ro − σ 2 of the cavity is presented in dashed line. The snapshot was generated with the VMD program [26]. As the statistical mechanical ensembles equivalence does not apply for few-particle systems, the approach is based on the canonical ensemble partition function (CPF). We go through the theoretical details in section II and devote section III to describe the simulation techniques and details. The combination of analytical and simulation results to obtain the expressions for Z 4 and the free energy is described in section IV. The mean properties of the system and the free energy for 2, 3 and 4 HS particles in a spherical cavity are given in section V. Final remarks and conclusions are presented in section VI.

II. THEORY
The system under study is a set of N HS-particles with hard repulsion distance σ, confined in a spherical cavity with radius R o . The particles interact with the cavity through a hard wall potential which prevents that particles escape to the outside. Thus, the effective radius of the cavity is R eff ≡ R o − σ/2, which represents the maximum possible distance between the center of the cavity and the center of each HS. The temperature of the fluid is determined by the wall temperature T , which is fixed.
The statistical mechanical and thermodynamic properties of the system are obtained from its CPF, Q N . One actually works with the configuration integral (CI), given that kinetic degrees of freedom integrate trivially. Thus, the CPF reads Here Λ = h/(2πm k B T ) 1/2 is the thermal de Broglie wavelength, while m is the mass of each particle, h is the Planck's constant and k B is the Boltzmann constant. In addition, Z N is the CI of the system that accounts for its non-trivial properties here e jk = exp [−βφ (r jk )] = Θ (r jk − σ) is the Boltzmann factor, φ is the HS pair potential, r jk the distance between the j and k particles, and the Heaviside function is Θ (x) = 1 if x ≥ 0 and Θ (x) = 0 otherwise. The inverse temperature is defined as β = 1/T (for simplicity we fix k B ≡ 1). The spatial domain of the integral for each particle is constrained to the spherical cavity. Naturally, Z N is a function of the cavity's radius R o and σ.
The CI for N = 1 is identified with the volume of the system through Z 1 = V (V = 4πR 3 eff /3), and the next two terms Z 2 and Z 3 were evaluated and discussed in Ref. [19,22]. In general, the structure of Z N is related to the Z i 's with i < N through a well known polynomial form (see p.135 in Ref. [27]), which also involves the inhomogeneous system cluster integral of N particles, τ N , as constant term. Alternatively, Z N can also be expressed as a polynomial in τ i 's with i ≤ N . τ N was introduced by Mayer in his cluster theory of homogeneous system with a different notation, and was latter used to develop the diluted gas equation of state (EOS), known as Virial series. A general form for the cluster integrals of the HS system confined in a spherical cavity is [22,28] where A = 4πR 2 eff , J = 8πR eff , K = 4π are the surface area and the mean and Gaussian extensive curvatures, respectively. S N includes the higher order terms in R −1 eff , while ∆τ N is the non-analytic core, being ∆τ N = 0 for R eff > (N − 1) σ. The constant coefficients {b N , a N , c N,J , c N,K } are the volume, area and curvature components of τ N . The known components up to order four are presented in Table I. For N = 2, 3 and 4 the CIs are [29] Two limiting regimes characterize the features of Z N independently of the number of particles in the cavity. For very large values of R eff (very small number densities) the CI scales as Z N ≃ V N and the system behaves basically as an ideal gas. In the opposite limit, for each N there exist a small enough cutoff value of R eff , the packing radius R m , such that R eff < R m implies that the hard spheres do not fit in the cavity. Thus Z N = 0. Close to this high-density limit, at R eff > ∼ R m , the behavior of Z N is characterized by its root, while τ N is finite and the pressure on the wall diverges as will be shown in the next section.
A. Free energy of a few body system The standard relations between CPF and thermal quantities, that are usually applied to systems with a large number of particles, also apply to few-body systems [30]. Of course, this statement is valid provided that any assumption about the large number of particles involved and the so-called thermodynamic limit are avoided [21,28]. Thus, for a system in stationary conditions with fixed N and T the Helmholtz's free energy (F ) is where U is the system energy and S, its entropy. The energy term of HS particles is simply βU = 3N/2. The reversible work done at constant temperature, to change the cavity radius between states a and b is equal to (minus) the change of free energy with dV = A dR eff . Here, P w is the pressure on the spherical wall which is an EOS of the system. The derivative of Eq. (7) at constant T gives which meets the exact relation known as contact theorem [30,31], In this ideal gas-like relation, ρ R eff is the value that takes the density profile at contact with the wall. This is an extended version of the contact theorem for planar walls. In this case, it applies to curved walls of constant curvature (spheres and cylinders), for both open and closed systems. Note that Eq. (9) follows directly from Eqs. (1,2,5,7) and (8), where in Eq. (2) it is convenient to explicitly write the Boltzmann factors for the external potential, i.e. Θ(R eff − r i ).
As it was mentioned above, the use of Eq. (1) to Eq. (9) is well established for many-body systems. Yet, these Eqs. also apply to the much less studied case of few-body systems. As this might be regarded as controversial, we provide here a few arguments to support it [21,28,32]. An important point is that statistical mechanics theory of ensembles does not exclude few-body systems [33]. In addition, independently of the considered number of particles, we relate the system properties with the corresponding mean-ensemble and thermal quantities such as energy, free energy, entropy. Thus, we follow the basic assumptions of statistical mechanics for many-body systems, except those that explicitly consider a large number of particles, to postulate the validity of Eqs. (5) and (6) for few-body systems. In particular, we avoid the ensembles equivalence, that is valid for large number of particles but does not apply to few-body systems. Finally, the ergodicity of the system is also assumed to claim the equivalence between the mean ensemble value of a magnitude and the simulation time average. Ultimately, the comparison of theoretical results with simulation or experimental outcome serves to validate our approach.

III. SIMULATION TECHNIQUES
As it was already mentioned, in few-body systems, the statistical mechanical equivalence between different ensembles does not apply. Therefore, simulation and statistical mechanical approches should correspond in detail to the same physical conditions, to ensure that the obtained results are comparable. In this work, we focus on a system at constant temperature. From the statistical mechanical theory point of view, it corresponds to a canonical ensemble, and thus one assumes a Maxwell-Boltzmann velocity distribution. On the other hand, the simulations should include a thermostating mechanism that ensures constant temperature and Maxwell-Boltzmann velocity distribution. The effect of utilizing different ensembles on the properties of few, spherically confined, hard-spheres was analyzed in Ref. [34].
Molecular dynamics simulations of few hard spheres confined in a spherical cavity were performed with a standard event-driven algorithm (EDMD) [25]. Constant temperature was achieved by using a thermal wallsthermostat, which changes the velocity of the particle colliding with the wall by means of a velocity distribution compatible with canonical ensemble for a given temperature T . We emphasize that this way of thermostating the system is unusual for bulk simulations. For these, local thermostats acting on each particle, like the Langevin thermostat, or a global thermostat that act on the system as a whole, by means of extra degrees of freedom and dynamical equations (i.e. Nosé-Hoover) are usually employed [35]. In our case, the particles are thermalized by collisions with the wall, in the same way in which a real system exchange heat with a confining medium, which acts also as thermal bath [36]. We consider that this procedure is very representative of the experimental case and, as we deal with very small systems, the collisions with the wall are very frequent which produce an efficient thermalization. We notice also, that the HSsystem is usually considered athermal because the interacting potential does not introduce an energy scale. For example, the isotherms P (T ) vs ρ collapse in a unique curve when plotted βP vs. ρ. In the simulation, however we consider the HS system in contact with a thermal bath which imposes its temperature to the system. In this line of thinking, we consider that the HS-system is in thermal equilibrium at a fixed temperature.
In this way, we account for two types of events, namely, particle-particle collision and particle-wall collision. For the particle-particle collision, the usual EDMD algorithm was used, in which the particle moves with rectilinear and constant-velocity dynamics, between particle collisions [25].
The time to collision of a particle i approaching a particle j is calculated as: where r ij ≡ r i − r j and v ij ≡ v i − v j are the relative positions and velocities of the particle pair, respectively. The parameter b ij ≡ r ij · v ij must be negative if the particles are approaching each other and positive otherwise, and we consider only the positive values of t ij . Here, it is assumed that the spheres are not already overlapping. The Eq. (10) is obtained by simply asking that at collision time t ij , the condition |r ij + v ij t ij | = σ must be fulfilled. An ordered list of events, with increasing collision times t ij is generated. Between collisions times, the particles move with r i = v i t. Once a collision occurs, the new velocities of the pair of particles involved in the collision are obtained as: [25]. In our case, however, it must be also considered the time for which each particle collides with the wall t w i . This time is calculated by the condition The nearest next event is chosen as the minimum of the next particle and wall events: min(min(t ij ), min(t w i )). If the particle-wall collision is the next event, the system is evolved until the particle reaches the wall. At this point, the thermal-walls thermostat acts on the particle by imposing it a new velocity, which is chosen stochastically from the probability distribution densities: here n and t stands for the normal and tangential components of the velocities, that lie in directions −r and v old − (v old ·r)r, respectively. The thermal walls described by Eq. (11) fix the temperature of the system. They were tested for HS confined by planar walls and produce a velocity distribution compatible with that of Maxwell-Boltzmann [37]. In the present study we have verified the same behavior for the temperature profile and the distribution of the velocities for curved walls. It should be noted that it is not correct to simply use a Gaussian distribution for the normal component of the velocity bouncing in the wall (Eq. (11)). This can be rationalized by understanding that the flux of particles leaving the wall is what must be sampled and not the distribution of velocities of particles in the wall [37,38].
In Fig. 2 we provide an example of temperature profile, for a number density ρ = 0.1, N = 3 and T = 2. The local temperature across the cavity agrees with the imposed value, and higher fluctuations and deviations are observed only close to the center of the cavity, due to poorer statistics. This feature is typical of all the profiles of local quantities for the spherical geometry. The Inset shows the distribution of velocities of the particles, giving an excellent fit to a Gaussian distribution of T = 2 (see caption). This behavior was verified for all the densities and number of particles considered in this work. All the simulation results presented throughout this work correspond to T = 2. Different number of particles were simulated using typically 4 × 10 7 particle-particle colliding events for each physical condition. The densities were varied in a broad range from very low densities to high densities, in the vicinity of the closed-packing condition.
During the simulation we evaluate various thermal quantities and spherically averaged functions of position. The studied structure and strain position functions are: the one body distribution function ρ (r) and the two independent components of the local pressure tensor: the normal P n (r) and tangent P t (r) components, respectively. Naturally, the calculated profiles for ρ (r), P t (r) and P n (r) are mean values over a discrete domain, obtained on a grid of spherical shells during the elapsed simulation time. The grid takes steps of ∆r = R eff /200. The definition for the local pressure tensor is that of Irving and Kirkwood (IK), i.e. a straight line path for the contour joining two centers of force [39,40]. The IK pressure tensor is well behaved in spherical coordinates and can be written as [41] P αβ (r) = ρ (r) k B T δ αβ + 1 2 where δ αβ is the Kronecker delta and α, β run over the local orthogonal unit vectors in spherical coordinatesr, θ andφ. The configurational contribution to the tensor involves the ensemble average · · · and each pair of particles contributes with where δ (x) is the Dirac delta function. During the simulation of the HS system, each collision contributes to the pressure with This contribution is distributed in M points located at r l (l = 1, · · · , M ), that are equally spaced along the line of length σ, between the pair of colliding particles. Each point r l , that lies in the k-th spherical slab of the grid with mean radius, contributes to P αβ (r k ) with Πr α ijr β ij /M where the projections on the spherical coordinates arer α ij =r ij .α (r l ). The two independent components of the diagonal tensor are P n = P rr and P t = (P θθ + P φφ ) /2. The discretized domain produces average local values that can be slightly different from the exact ones. Particularly, in the near-wall and central regions, the binning of different volume produce a poorer statistics for smaller r, which can be noticed in the profiles.

IV. SIMULATION-ASSISTED QUASI-EXACT ANALYTICAL APPROACH
The CI of the system of four HS confined in a spherical cavity Z 4 (R eff ) reduces to an integral in nine variables, as it is also the case of the fourth cluster integral τ 4 (R eff ) of the system. A taste of the difficulty of solving these integrals can be gained by recognizing that, for the case of homogeneous HS system where τ N = N !b N V , the highest order cluster integral exactly known is the fourth (the constant b 4 reduces to an integral in six variables). [42] On the other hand, the higher order cluster integral known for an inhomogeneous HS system is τ 3 (R eff ) for spherical confinement. This integral also reduces to a six-variable integral [22]. In this work, we develop a simulationassisted procedure that uses the mean simulation value of wall pressure of the 4-HS in a hard-wall spherical cavity to obtain, based on fitting, an analytic expressions for both τ 4 and Z 4 . These expressions accurately describe both functions along the complete range of cavity sizes R eff > 0.6124. In this section, we set σ ≡ 1 to deal with simpler notation.
Cluster integrals of HS system in a spherical confinement are non analytic functions of R eff . If we analyze from the infinitely diluted configuration by decreasing the cavity radius, the first analytic branch of τ N extends from to R eff → ∞ to R eff = (N − 1) /2. There, a discontinuity in a higher order derivative may appear. For example, at R eff = 1 τ 3 R eff has a discontinuous fifth order derivative. [22] For even smaller R eff , a sequence of non-analytic points appear in τ N until the highest density packing is attained at R eff = R m , where Z N collapses to zero and τ N takes a finite value. For each N , the number density ρ 0 = N/V reaches a different maximum driven by the geometry of the packing configuration, that exceeds the packing density of the bulk system ρ 0 = 1.4142. The system with N = 2 attains its maximum density ρ 0 = 3.8197 at R m = 0.5, for N = 3 it is ρ 0 = 3.7215 at R m = 0.5773, while for N = 4 it corresponds to ρ 0 = 4.1584 at R m = 0.6124.
By using the analytical expressions (4), we obtained Z 4 by fitting simulation results of P w in the low and high density branches in a process that we call simulationassisted calculation. A detailed explanation of this procedure is provided in Appendix A.
A peculiarity of the N = 4 system is that an ergodicity break is produced when the cavity radius decreases below R eff = 0.7071, that corresponds to ρ 0 = 2.70095. This purely classical effect, that is absent for N = 2 and 3, appears for distinguishable particles, but is absent for indis-tinguishable ones. In the ergodicity break, the available region of the spatial domain of the particles separates into two identical disjointed regions. This produces a discontinuous step in the CI that decreases by a factor 1/2. This curious behavior appears because, for R eff < 0.7071, it is impossible the exchange of positions between two particles, keeping the other two fixed. The spatial constraint forces the 4 HS to move as it were a tetrahedron, that can only rotate and have small deformations. Alternatively, once the labeled particles are constrained to such a small cavity, there is no path that allows to go from a given available state to another one, obtained from its speculate reflection. Systems of 2 HS confined in cylindrical and cuboidal pores have also shown an ergodicity break behavior [21].

V. MEAN PROPERTIES AND FREE ENERGY OF CAVITIES WITH 2, 3 AND 4 PARTICLES
Even at this very small number of particles in the studied system, we define a bulk-like behavior to compare the mean and thermal quantities of the system, with their macroscopic homogeneous fluid counterparts. These bulk-like behavior is found typically at intermediate densities in the center of the cavity, where an approximate constant density profile is observed. In many properties, we found a qualitatively different behavior, depending upon the number of particles, whose root is a different (geometrical) packing with varying number of particles. We will discuss these differences while presenting the results. Figure 3 shows the density profiles across the spherical cavity from the center (r = 0) towards the wall of the cavity for three distinctive total number densities ρ 0 . The density profiles are scaled with the number density of the system ρ 0 and the r coordinate, with the effective radius of the cavity R eff to allow comparison between different number of particles. For 2, 3 and 4 particles, the overall features are similar. At low number density (ρ 0 = 0.01), the profile is constant (homogeneous) across the pore, indicating ideal gas behavior. At intermediate densities, represented here by ρ 0 = 0.35, a clear homogeneity in the density profile shows-up, with an increase of density close to the wall. A central region, which we call plateau, of approximately constant density is still present, but its extension decreases towards the center of the cavity for higher number of particles. At high number densities (ρ 0 = 1.40σ −3 ) the inhomogeneous profile persists, with the distinctive feature of vanishing density in the central region. In this high packing regime, the particles cannot occupy the center of the cavity, due to the presence of the other ones. In the case of N = 2, a direct comparison between theory (open circles) and simulation for ρ(r)/ρ 0 along the entire range of r is presented. The agreement between simulation and theory is excellent. This strengthens the confidence in the thermostating scheme used in simulations and the overall theoretical and numerical approaches. We also checked temperature conservation with the temperature profiles (see Fig. 2) across the cavity, obtained by the application of the equipartition theorem to the local kinetic energy of the particles. In the case of N = 3, only ρ center and ρ(R eff ) were found theoretically. For N = 4 the simulation-assisted approach and Eq. (9) enables to obtain ρ(R eff ).
In Figure 4, we present details of the scaled density profiles for an intermediate value of number density (ρ 0 = 0.35σ −3 ) and N = 2, 3, 4. Open circles for N = 2 and for N = 3 at r = 0, correspond to exact theoretical results. Diamonds show the scaled wall pressure βP w /ρ 0 , that should be identical to the scaled density ρ(r)/ρ 0 at the wall in accord with contact theorem [Eq. approach. The figure also shows, with open diamonds, the scaled wall pressure βP w /ρ 0 obtained with MD. They were shifted to the left and a horizontal line was draw to guide the eyes. The agreement with the contact theorem and between theory and simulation are, again, both excellent. These results are very important to validate the simulation-assisted fitting, used for quasi-exact calculations in the case N = 4. The profiles have similar features, but they have significant differences as changes the number of particles. A main common feature is the enhancement of density towards the wall. This can be understood as an effective attraction produced by the translational entropy maximization. The slightly enhanced presence of particles close to the wall, reduces the mean excluded volume, which increases the number of possible configurations, and therefore the entropy of the system. This can be also considered as a special case of the depletion interaction [43]. Another interesting feature, noticeable for N = 3 and 4 is that the density profile presents a local minimum, and also structure in the central region of the cavity. This features in the density profile were also observed by White et. al. for Monte-Carlo and density functional theory calculations of N = 6, 8, 10 [44].
In Figure 5, we show pressure on the wall for N =2,3,4 particles as a function of mean number density. For N = 2, 3 the high degree of coincidence between simulation (circles and squares) and theoretical results (continuous and dotted lines) can be observed. In the case of four particles, one can observe that the simulation results (diamonds) are well described by that found using the simulation-assisted approach (dashed line).
In Figure 6 we show an estimator of the system- This quantity can be regarded as a measure of the interface tension for these systems of few particles, far from the thermodynamic limit. Symbols show simulation results and lines between symbols are a guide for the eyes. For comparison, the curve in full line (blue) shows the surface tension for the bulk HS system in contact with a planar wall.
substrate surface tension γ est ≡ (P center − P w ) R eff /2, as a function of number density, based on the Laplace equation [45]. At low densities γ est is close to 0 for all the number of particles considered, showing a vanishingly small energetic cost to confine the system. The inset in Fig. 6 is a zoom showing how the curves depart from 0 upon increasing confinement. A local minimum of γ est is observed for N = 2. At higher number densities, there is an interesting important difference among the cases N = 2 and N = 3, 4. This qualitatively different behavior is a consequence of the different packing features, when the system becomes strongly constrained. There are limiting number densities ρ c 0 = 1.10266 (N = 3) and ρ c 0 = 1.47021 (N = 4) such that P center vanishes for ρ 0 > ρ c 0 . In these cases, there are no possible configurations of the system such that the segment joining the center of two colliding particles, also passes through the center of the cavity. This geometrical constraint becomes evident at lower densities ρ 0 < ∼ ρ c 0 , where there are few configurations that produce non-vanishing pressure in the center. At higher densities, the pressure difference diverges as −P w . The case of two particles does not has this geometrical restriction and then P center diverges with increasing density. In this case, the opposite mechanism takes place: when the caging is high and a lower total number of configurations is possible, a progressively higher number of them have inter-particle segment going through the center of the cavity. This divergence of P center is stronger than that of P w and thus γ est diverges to +∞ for N = 2. For comparison, the system-substrate surface tension of the infinitely extended HS system in contact with a planar hard-wall is also shown [24]. It is negative and its behavior follows that found for N = 2, 3, 4 at very low densities. The attained maximum density is in this case ρ 0 = 1.4142, the maximum possible density of the bulk-solid system.
We also studied the pressure tensor profile across the spherical cavity. As an example, Figure 7 presents the profiles of the normal and tangential components of the pressure tensor for N = 3. In the center of the cavity r ∼ 0 for low (ρ 0 = 0.01) and intermediate (ρ 0 = 0.35) densities, the behavior P t = P n = 0 is representative of a zone of bulk-like behavior. In the interior of an homogeneous and isotropic liquid all the components of the pressure tensor have equal magnitude. Upon increasing number density, the bulk-liquid region is reduced towards the center of the cavity and the region in which there are important differences between the normal and tangential components of the pressure tensor increases. This can be regarded as the interfacial-like region of the cavity in which the differences in the components of the pressure give rise to a non-vanishing contribution to the surface tension as was shown in Fig. 6.
At high densities the pressure in the central region is usually zero (P t = P n = 0), as it is observed in the panel for ρ 0 = 1.40 in Fig. 7. In line with the surface tension, this can be rationalized by thinking in an event of colliding particles. For a given pair of colliding particles, the contribution to the local pressure is assigned along the straight line which connects the pair. If, due to geometrical restrictions of confinement and excluded volume, all of these lines lay out of the central region, the contribution to the local pressure tensor components vanishes. At these high densities, there is no bulk behavior at all, and the center of the particles are not allowed to move through the center of the spherical cavity. It can also be the case (not shown) in which the density in the central region is 0 and the local pressure have a non-vanishing value. This corresponds to a high-packing case, in which the center of particles cannot be close to r ∼ 0, but the lines connecting a colliding pair can pass through the  central region. This happens, for example, in the case of N = 2 and R eff < 1 or number density ρ 0 > 6 4π ≃ 0.48. In Figure 8 the adsorption Γ = (N − ρ center V ) /A is shown. The data from simulation is depicted with symbols and the theoretical results with lines. This quantity measures the excess of number density in the surface, as compared to a gas with mean number density ρ center . The curves present theoretical results for N = 2, 3 and 4. For N = 4, adsorption is drawn from the ρ 0 value such that that ρ center = 0. An excellent agreement between theory and simulation is observed. At low densities the adsorption tends to 0 for the three cases we studied, in the limit of ideal gas. Then a fast rise of adsorption is observed, in which the curves for different number of particles are concave and cross each other (see inset of Fig.  8). The concave region ends up when ρ center vanishes. At high densities, the adsorption is higher for greater number of particles with a monotonic behavior in all the cases. This goes in line with with the behavior of the sur-  face tension for N = 3, 4 (see Fig. 6), in which a smaller surface tension gives rise to a higher adsorption in the walls of the cavity. Figure 9 presents the (Helmholtz) excess Free energy F ex = F − F ideal per particle for the three cases considered in this work. The exact analytical results are given for N = 2, 3 and the simultation-assisted result is provided for N = 4. It is worth noticing that this curve is very difficult to obtain directly by simulations for such a wide range of number densities. Within the simulation approach, this would mean performing a numerical thermodynamic integration for each point of each curve, with typically ten to twenty runs per free energy value [46]. The combined simulation-assisted approach provides this thermal quantity in all the relevant range of number densities. The excess free energy per particle for a hard sphere system is, of course, F ex = −T S ex . Therefore, it gives an idea of how the entropy is reduced for the different number of particles, when the number density increases in the spherical cavity, and the available number of possible states of the system reduces. The free energy curves cross each other, showing the influence of the confinement in the free energy of a few particles at high densities.

VI. CONCLUSIONS
In this work we studied the properties of 2, 3 and 4 hard spheres confined in a spherical cavity at constant temperature. We used analytical calculations and event driven molecular dynamics simulations to study the system for a wide range of number densities. We described the system with variables and concepts widely used in thermodynamics and macroscopic theory of liquids, and provided a comprehensive characterization, covering very different physical regimes from ideal gas to highly-packed. The results we obtained are interesting in the context of theory of liquids, but also for new developments in material science and synthesis of hollow nano-particles and rattles. In this minimal system, only the excluded volume is present, and changes in free energy can only be due to variations of translational entropy. This range of phenomena is common to any experimental system and will have a role worth-understanding, independently of the characteristics of the confining cavity or the nature of the confined colloids or nano-particles. We characterized the density and pressure tensor components profiles as a function of number density and analyzed the distinctive behavior of an estimation of the surface tension for different number of particles. The adsorption was also characterized. In many quantities we found a rich and non-trivial behavior, which was ultimately rationalized with the geometrical differences of these few body systems. We compared carefully analytical and simulation results whenever possible, to validate both methodologies. A simulation-assisted approach was also used to obtain an analytical quasi-exact expression for the free energy and pressure on the wall for the case N = 4, where a full analytical calculation was not found. Using this, we calculated and compared the free energy of the system for all the studied number of particles.
This study will be extended by adding attraction among particles, which is another important feature present in real systems. In this context of minimal models, this will be achieved by using square-well potentials, which in addition to the hard wall repulsion, have a short range attraction. Also the properties of the wall can be modified to model more realistic systems, for example adding a wall attraction. This could be an interesting direction to understand these highly confined systems.

ACKNOWLEDGMENTS
We thank Iván Paganini for fruitful discussions. The financial support through grants PICT-2011-known terms for τ 4 (see Eqs. (3), (4) and Table I). We approximate τ 4 as its truncated form to order R −2 eff by assuming  5) and (8) is obtained the analytic expression to fit βP w R eff . To ensure ∆τ 4 = 0 we should restrict the fit to values R eff > 1.5. However, we obtain a better description analyzing the low density behavior up to R eff = 1.4. Performing this careful fit, our best estimation is c 4,J = −1.5162(500). Once c 4,J is fixed, we analyze fitting functions with two parameters. By fixing c 4,−1 = 0, we found our best fitting coefficients c 4,K = −1.51495(600) and c 4,−2 = 0.4181(80) with a regression coefficient r 2 = 0.999999997. The four HS system in a spherical cavity reaches the highest number density at a tetrahedral packing configuration. Focusing in the limiting behavior of the wall pressure we fit where d 4i stands for the four fitting parameters. The best estimate for d 40 is 8.99145(80) which we replace by the natural number d 40 = 9. Once d 40 is fixed we obtain our best estimate for d 41 given by −7.1745(600) which we represent by d 41 = −43/6. The other two fitted coefficients were d 42 = 88.7675 and d 43 = 414.303 with uncertainties of 1.0 and 15.0, respectively. Using the obtained low density branch of βP w we evaluated its order four series expansion at R eff = 1.4. A fitting method was implemented that reproduce the local behavior of βP w at R eff = 1.4, its asymptotic properties at R eff → R m , and use three extra parameters to fit βP w for R ef f ∈ (0.6124, 1.4). In this way we obtained an expression for βP w (R eff ) for all the range of possible values of R eff ∈ (R m , +∞). The obtained function has two analytic branches that smoothly join at R eff = 1.4 with a discontinuity at the fifth order derivative. The expression of Z 4 for the high density branch, that has an extra non-analytic point (a discontinuity) at R eff = 1/ √ 2 (seeIV), is found from Eqs. (5) and (7). The expression for R eff < 1.4 is with pf = 1 if R eff ≥ 1/ √ 2, pf = 1/2 if R eff < 1/ √ 2 and q 1 = −43/6. The remaining coefficients are given in Table II