Deformable hard particles particles confined in a disordered porous matrix

With suitably designed Monte Carlo simulations we have investigated the properties of mobile, impenetrable, yet deformable particles that are immersed into a porous matrix, the latter one realized via a frozen configuration of spherical particles. By virtue of a model put forward by Batista and Miller [Phys. Rev. Lett. {\bf 105}, 088305 (2010)] the fluid particles can change under the impact of their surrounding (i.e., either other fluid particles or the matrix) their shape within the class of ellipsoids of revolution; such a change in shape is related to an energy change which is fed into suitably defined selection rules in the deformation"moves"of the Monte Carlo simulations. This concept represents a simple, yet powerful model of realistic, deformable molecules with complex internal structures (such as dendrimers or polymers). For the evaluation of the properties of the system we have used the well-known quenched-annealed protocol (with its characteristic double average prescription) and have analysed the simulation data in terms of static properties (radial distribution function and aspect ratio distribution of the ellipsoids) and dynamic features (notably the mean squared displacement). Our data provide evidence that the degree of deformability of the fluid particles has a distinct impact on the aforementioned properties of the system.


I. INTRODUCTION
Experimental, theoretical, and computer simulation based investigations on fluids confined in disordered, porous confinement have received over meanwhile several decades a steadily increasing share of interest [1][2][3][4]. This is certainly due to the fact that the scenario of confined fluids is ubiquitous in our daily lives, with examples ranging from biophysics over chemical engineering to technological applications. But also from the pure academical point of view such systems have been (and still are) studied extensively in fundamental science which is due to the fact that confined fluids drastically change their properties once they are exposed to the external field of a disordered matrix: this applies not only to the static properties (such as the structure, the thermodynamics, or the phase behaviour) but also to the dynamic properties of the fluid (such as the mean squared displacement of the particles -and hence the diffusivity -, or the dynamic structure factors).
From a theoretical point of view -be it within theoretical frameworks or in computer simulations -such investigations are particularly challenging both from the conceptual as well as from the numerical point of view: it is difficult to properly define theoretical concepts that allow the reliable evaluation of thermodynamic properties of such systems; and it is also challenging to put forward models -both for the liquid as well as for the porous confinement -that faithfully mimic the features of real systems (such as He confined in aerogels [5,6]).
In many of the preceding theoretical investigations the seminal quenched-annealed (QA) model was used: here the confining environment -henceforward termed matrix (introducing the index 'm') -is modeled as an instantaneously frozen configuration of an equilibrated liquid; further, the space left void by this matrix is filled by mobile, fluid particles (with index 'f'). This concept represents a viable compromise for the above mentioned requirements: it is amenable both to computer simulations and theoretical frameworks, while it still captures relevant features of realistic systems at a rather faithful level. The related theoretical, statistical mechanics based framework was pioneered by Madden and Glandt [7,8] and later by Given and Stell [9][10][11]: within these concepts the system is considered as a very peculiar mixture of the matrix and of the fluid particles. Observable quantities are calculated via a double-averaging procedure: in a first step physical properties are obtained via averaging over the degrees of freedom of the fluid for a given matrix configuration; subsequently the obtained values are averaged over all possible matrix configurations that are compatible with 2 the macroscopic system parameters. Within this formalism the correlation functions between the particles are related via an Ornstein-Zernike type equations (the so-called Replica Ornstein-Zernike equations) which can be solved numerically in combination with a suitable closure relation. This approach provides furthermore access to the thermodynamic properties of the system [12][13][14][15][16][17]. And it should be mentioned that Krakoviack [18][19][20] successfully merged this formalism with the concept of mode coupling theory [21], paving thus the way to investigate also dynamic properties of confined fluids. Computer simulations, on the other hand, literally execute the above outlined double averaging procedure: typically five to 20 different, but equivalent matrix configurations are selected to realize the "outer" averaging procedure. All in all, numerous studies have been performed during the past decades on QA systems to investigate their static [15,[22][23][24][25][26], thermodynamic [12-14, 16, 27], or dynamic [28][29][30][31][32][33][34]36] properties.
In most of these investigations both types of particles are assumed to be spherically symmetric and fixed in their shape. This assumption is presumably not too restrictive for the matrix particles as they are fixed in their positions, anyhow. However the situation is different for the mobile, fluid particles: in particular in the realm of soft matter physics these particles are usually viewed as coarse-grained, "effective" models whose shape can in principle vary as a consequence of the internal dynamics of the constituent entities of such molecules. In the vast majority of scientific investigations such "effective" particles are considered as spherical; still, in a few contributions the non-spherical shape of these molecules has been taken into account, for instance, in polymers [37][38][39][40], dendrimers [41,42], or patchy particles [43]. Within the framework of "effective" particles it is therefore not surprising that the shape of such "effective" particles adopt their shape to their immediate surrounding -be it another mobile particle or the rigid matrix.
In an effort to take this feature within the QA scenario properly into account it is more appropriate to introduce to a certain degree flexibility in the shape of the mobile particles.
In this contribution we address this issue for the first time, taking benefit of the availability of a simple, yet realistic and powerful model of particles of variable shape: our investigations are based a model of deformable, impenetrable particles (proposed by Batista and Miller [44,45]) which are allowed to change -at the cost of some energy penalty -their shape in computer simulations: to be more specific, an initially spherical particle can be deformed into an ellipsoid of revolution, characterized by its aspect ratio x. As specified in Refs. [44,45] the energy penalty for deformation is proportional to βκ ln 2 x, where βκ is the deformability (or stiffness) of the particles (β being the inverse temperature). Batista and Miller have shown that this model is readily amenable to Monte Carlo simulations, introducing suitable, additional rotational and particle deformation "moves." We have integrated this model into the framework of QA systems: the mobile, impenetrable particles are now deformable, based on the above outlined model of Batista and Miller. We apply in our Monte Carlo simulation in an alternating order the following MC "moves" to the fluid particles: (i) conventional particle displacement, (ii) rotational moves of the ellipsoidal particles, and (iii) deformation moves. Now that the particles are no longer spherical the verification of the overlap criterium is computationally more cumbersome and considerably more time consuming than in the spherical case: here we have taken benefit of a concept proposed by Vieillard-Baron [46] in his seminal work and have implemented the related, numerically very efficient algorithm proposed by Perram and Wertheim [47].
For the "outer" average of the above mentioned double averaging procedure within the QA framework we have assumed -as a consequence of the high computational cost of the extended steps -five independent matrix configurations. Particle ensembles typically range from 3000 to 4000 particles.
In order to quantify our findings we have focused on the following features: the static radial distribution function and the distribution of the aspect ratio of the fluid particles, and their mean squared displacement. Investigations have been carried out over a significant number of states in the parameter space, spanned by the number density of the matrix particles (Φ m ) and of the fluid particles (Φ f ) and considering a broad range of deformability parameters βκ.
The manuscript is organized as follows: in the subsequent Section we present our model and provide details about the simulation technique; we furthermore define four pathways through parameter space along which we have investigated properties of our system. The relevant data are discussed in Section III, with particular emphasis on the radial distribution function of the deformable, mobile particles, the emerging probability distributions of the aspect ratios of the particles, and on their mean squared displacement. The manuscript is closed with concluding remarks and an outlook to future work.

A. Model
We consider initially an ensemble of N hard spheres of diameter σ, confined in a volume V at a temperature T (with β = 1/(k B T ), k B being the Boltzmann constant); σ is the unit length of the system and is henceforward set to unity. N m of these particles are used to form the matrix and are henceforward kept fixed in their positions. Furthermore, these matrix particles are undeformable, i.e., they keep their spherical shape; their packing fraction, φ m , is defined as The remaining N f = (N − N m ) fluid particles are immersed into the space left void by the matrix configuration; these particles are mobile and deformable: following the concept of Batista and Miller [44,45] they can -subject to an energy penalty -deform into (still impenetrable) prolate or oblate ellipsoids of revolution. Their packing fraction is defined by The framework behind our model is that of a quenched-annealed (QA) system, which has been designed to investigate the properties of fluid particles immersed into a confinement realized via a disordered particle arrangement: the matrix is assumed to be an instantaneously frozen configuration of an equilibrated fluid with its voids being filled by fluid particles [9,11,12,32,34]. The concept of QA systems allows us to calculate thermodynamic averages of observables: working in the canonical ensemble (and following the ideas of previous contributions) the properties of the system (in particular its static and dynamic behaviour) are calculated via a double averaging procedure (see, e.g. [32,34]): assuming a The concept of deformable hard spheres rests on an idea put forward by Batista and Miller [44]: starting from spheres with diameter σ the particles can change their shape and can become ellipses of revolution (spheroids), with semi-axes a = b = c; it must be emphasized that during these shape transformations the volume of the particles is preserved, therefore the particles can be considered as incompressible. For convenience we introduce the aspect ratio x, defined by x = a/c; thus particles are oblate for x < 1, prolate for with orientations Ω 1 and Ω 2 ) show after such a move an overlap or not: in the former case, the particle move is rejected while it is otherwise accepted. For spherical particles the implementation of such a criterion is trivial; however, this is not the case for spheroids: in our program we have implemented a reliable and numerically efficient method for identifying the possible overlap of two ellipsoids, put forward by Perram and Wertheim [47]. This method is numerically more attractive than an implementation of the original criterium put forward by Vieillard-Baron [46]. For the application of this criterium the spherical (matrix) particles are considered as spheroids with x ≡ 1.

B. Monte Carlo simulations
Our canonical Monte Carlo simulations were carried out in a cubic cell, assuming periodic boundary conditions. Three types of "moves" were applied to the particles: (i) conventional particle displacements; the maximum value of displacement is adjusted dynamically during the simulation to guarantee an acceptance rate of the translational moves of 50 %; (ii) rotational moves of the particles; in a trial rotation of a particle its Euler angles are changed by angular increments, whose maximum values are dynamically adjusted to ensure an acceptance rate of the rotational moves of 50 %; (iii) shape deformation "moves" of the particles; following the concept of Batista and Miller 6 [44,45], the potential U (x) which governs the particle deformations is given (in lowest order approximation) by (1) x is the above introduced aspect ratio of the particles and κ is the stiffness (or deformability) parameter; summarizing the details laid out in [44], κ is related to the surface tension γ of the particles via κ = 8 45 σ 2 πγ.
In a shape deformation "move" the ln x-term in Eq. (1) is changed by an additive "displacement" ∆x shape ; this value is dynamically adjusted to ensure an acceptance rate of the shape deformation "move" of 50 %; for details we refer to Refs. [44,45].
Acceptance of all three "moves" is governed by the aforementioned overlap criterion: thus, if a trial move leads to an overlap between two particles, this move is definitely rejected.
For the shape deformation "move" an additional, energy-based criterion is imposed which reads (for details cf. Refs. [44,45]): Thus a shape deformation "move" is accepted -provided that no particle overlap occurs for the deformed particles -by a probability p specified in the above equation.
In our notation a sweep consists of a series of subsequent N f translational, N f rotational, and of N f shape deformation "moves" of randomly selected fluid particles. At this point it should be emphasized that the check for particle overlap is computationally quite expensive (imposing thus limitations on the ensemble size).
In view of the above said, our system is characterized by three parameters: the deformability κ and the number of fluid (N f ) and matrix (N m ) particles. With our choice of κ-values (i.e. βκ = 1, 5, 10 and 100) we cover a broad spectrum ranging from strongly deformable particles (with low βκ-values) to undeformable particles (with βκ = 100); the case of spherical hard spheres is recovered for βκ → ∞.
Simulations have been carried out in ensembles of typically 3000 to 4000. The actual values of N f and N f depended on the location of the investigated state in the parameter space spanned by Φ f and Φ m : in an effort to guarantee a sufficient numerical accuracy of the data we considered in general a number N α = 4000 × Φ α (α standing for 'f' or 'm') of either particle species: thus, for instance, for a system with Φ m = 0.20 and Φ f = 0.30 we considered N m = 800 matrix and N f = 1200 fluid particles. Only for low-density states (i.e., if either of the densities was 0.05) we used N α = 8000 × Φ α (α standing for 'f' or 'm'), instead.
In an effort to scan the relevant regions of the parameter space spanned by Φ m and Φ f we have defined four pathways (denoted by I, II, III, and IV); for this choice we were guided by the diagram of states of spherical fluid particles confined in a matrix of immobile, spherical, particles, as depicted in Fig. 1 of Ref. [32]. These four pathways are shown -along with symbols representing the specific, investigated states -in Fig. 1 and are defined as follows: Note that along each of these pathways the total packing fraction (i.e. Φ tot = Φ f + Φ f ) is less than 0.55, a value which is close to the coexistence density between the liquid and the FCC phase of hard spheres.
In Fig. 1 we also show (as a reminder) the "dynamic" phase diagram of spherical particles, confined in a porous matrix, formed again by spherical particles; this information should only be considered as a help of orientation -for more details cf. Fig. 1 of Ref. [32]. are shown in different colours. The grey lines sketch the "dynamic" phase diagram of a previous study [32] where the dynamic properties of spherical particles, confined in a porous matrix, formed by undeformable, spherical particles, were investigated.
one one side it is the higher numerical costs to identify particle overlap of the spheroids, on the other side it is the fact that we now have to include with βκ an additional system parameter: The protocol of the simulations can be summarized as follows. (i) First, we created the matrix configurations: for a given value of Φ m we performed simulations for N m spherical particles in a simulation box of volume V : starting from a random initial configuration, the system was simulated until equilibration was achieved; such an equilibrated particle arrangement was stored as a possible matrix configuration; continuing the simulation over

III. RESULTS
In the following we discuss the results obtained along the four paths through the (Φ m , Φ f )plane, as specified above and as visualized in Fig. 1. Our discussions are based on the radial distribution function, g(r), the probability distribution of the aspect ratio x of the fluid particles, the so-called aspect ratio distribution (termed ARD), the mean squared displacement (MSD), δr 2 (t), of the fluid particles, and the effective exponent of the MSD, Throughout the time t is "measured" in terms of MC-sweeps.

A. Path I and path II
Since results obtained for our system along path I and II display only small quantitative differences we discuss these data together in one subsection.
Along these pathways the packing fraction of the matrix assumes rather low values, namely Φ m = 0.05 and Φ m = 0.10, respectively. In Fig. 3 we have summarized the results for g(r) obtained for states along path I: now that the centers of two particles can approach each other -depending on their shape -to distances smaller than σ the onset of the main peak of g(r) can assume values that are smaller than σ. As a consequence the typical discontinuous peak in g(r) as observed in a hard sphere system, is now a steady function of the distance r: for obvious reasons this softening is the more pronounced the smaller the value of βκ, i.e., the stronger the particles are deformable (see also the discussion of the ARD below); in contrast, for βκ = 100 the main peak in the g(r) is still very similar to the one found in a pure hard sphere system. Further, also the onset of this peak shows a substantial shift to smaller distances with decreasing βκ which is essentially independent of the fluid packing fraction Φ f : for βκ = 1, for instance, the onset of the main peak occurs already at r 0.7 σ. Analysing the peak position in g(r) we note that for βκ = 100 this position is essentially unaffected by the fluid packing fraction, while for strongly deformable particles (e.g., for βκ = 1) a clear shift in the position of the peak of g(r) towards larger distances with increasing Φ f is observed: the strong increase in the preferred distances between strongly deformable particles indicate that they tend to occupy the available space within the matrix more efficiently; those regions are -for geometric reasons -not accessible for less deformable (or even rigid, spherical) particles, as they are not able to access narrow pores within the matrix. To conclude we point out that for the highest value of Φ f and for the highest value of βκ a shoulder in the side peak of the radial distribution function occurs: this feature indicates the onset of a glassy state, a feature which is agreement with the related investigations of spherical particles confined in a disordered matrix of spherical particles (and widely discussed in Refs. [32,34]). Since this shoulder is not observed for more deformable particles (see the other panels of Fig. 3 for smaller values of βκ) we conclude that a strong deformability of the mobile particles prevents our system from forming a glassy state.
Several of the above features are nicely and consistently complemented by analysing the data of the ARD, which we now discuss along path II; the related data are summarized in Fig. 4. These results indicate in an unambiguous manner how the fluid particles do deform under the influence of the surrounding matrix and of other the fluid particles. Our observations can be summarized as follows: the shape and the width of the ARD as a function of the aspect ratio x is primarily triggered by the stiffness parameter βκ, while on a more subordinated level these features are influenced by the available space, notably by Φ f : stiff particles (for instance with βκ = 100) show an ARD that is strongly peaked around x = 1 with a shape that essentially does not depend on Φ f . As the particles become more deformable (e.g., for βκ = 10 or 5) the ARD is still centered around a value close to x 1, however, the shape of the distribution function broadens considerably; still no sizable shape dependence on Φ f can be observed. Only for strongly deformable particles (i.e., for βκ = 1) the impact of the fluid packing fraction Φ f becomes visible: (i) we find very broad ARDs whose peak positions have shifted to smaller x-values -e.g. for Φ f = 0.05 we find x max 0.6 (corresponding to pronounced oblate particles); (ii) concomitantly, the width of the ARD broadens strongly with decreasing fluid packing fraction: for βκ = 1 and Φ f = 0.05 the "wings" of the ARD still assumes sizable values for x 0.2 and x 2.5: for obvious reasons particles with such extremely small aspect ratios are now able to populate narrow pores inside the matrix which are inaccessible to particles if they were spherical -which brings us consistently back to our above discussion of g(r). On a secondary level, the degree of subdiffusivity is triggered by the fluid density Φ f : for a given βκ-value this dependence is the more pronounced the denser the fluid particles. the mobile particles.

B. Path III
Data for the radial distribution function are shown in Fig. 6. Similar as for the related data collected along the other pathways we observe that the deformability of the particles induces with decreasing βκ a strong softening of the main peak in g(r) and a considerable shift of the onset of this peak towards smaller distances. The related data for the ARD are shown in Fig. 7. Compared to the preceding results we find similar features and dependencies: (i) the lower the value of βκ, the broader the ARD as a function of the aspect ratio x; (ii) only for strongly deformable particles (i.e., for βκ = 5 or even smaller) the packing of the particles (now imposed by the packing fraction of the matrix, Φ m ) starts to have a distinct influence on the shape of the ARD and its peak position shifts with decreasing Φ m to smaller x-values, located at x 0.55.
For the MSDs (shown in Fig. 8) we identify distinctive differences as compared to the data accumulated for states along path III (see Fig. 5). We start again with the longtime diffusivity (typically encountered for 10 4 t). Again the main impact on the MSD originates from the packing fraction of the matrix: as we increase Φ m from 0.05 to 0.20 (i.e., by a factor of four), the distance that the particles are able to cover over a comparable time range differs by more than three orders of magnitude (compare differently coloured curves in each of the MSD-panels of Fig. 8). In contrast, the MSD essentially does not depend on the deformability: particles cover for a given Φ f -value and again after a comparable time span essentially the same distance, irrespective of the deformability βκ (compare equally coloured curves in the different MSD-panels of Fig. 8). Thus we conclude that the long-term diffusivity of the particles within the surrounding matrix is triggered by confinement, while the deformability plays (at least for the value of Φ f that we have chosen) only a negligible role. From the slope of the MSD in the linear regime we can read off the diffusion constant: we observe that for those system parameters where the diffusive behaviour has been regained for the fluid particles of our system, calculated along path IV and considering various combinations of the fluid packing fraction, φ f , and the stiffness parameter, βκ (as labeled).

D. The radial distribution functions across the pathways
To conclude our analysis we have collected the radial distribution functions along all pathways for hardly (i.e., βκ = 100) and for strongly (i.e., βκ = 5) deformable particles and for all density parameters available in Figs. 9 and 10, respectively. The aim of this discussion is to disentangle the impact of the deformability, βκ, from the density parameters, Φ m and Φ f , on the shape of the distribution function.
From the accumulated data it is obvious that βκ is primarily responsible for the shape of g(r), in particular for the main peak: a large value of βκ (Fig. 9) leads to a harshly repulsive peak with an onset close to σ; the position of this peak changes only marginally as either of the densities (i.e., Φ m or Φ f ) are changed. Side shoulders in the second peak of g(r) for systems with high total packing fraction (Φ tot = Φ f + Φ m ) indicate the occurrence of a glassy state for the case of hardly deformable fluid particles. In striking contrast, we find for strongly deformable particles (i.e., βκ = 5, with data shown in Fig. 10) that the main peak of g(r) is now a smooth function with an onset at distances r 0.8σ (or even smaller -cf. Figs. 3 or 6). But also the position of the main peak is no longer restricted to values close to σ, but can reach -depending on the density parameters -values as high as r 1.25σ. We conclude that the strong deformability of the particles which can either lead to center-to-center distances between particles that are smaller than σ, but also to separations that are larger than σ. Possibly the particles can also explore the available space inside the pores of the matrix more efficiently as they can access in their deformed shape spaces that otherwise would be inaccessible for rigid, spherical particles; however, a closer analysis of this conjecture would require a detailed geometric analysis of the voids inside the matrix in terms of a Delaunay decomposition [35], an investigation which would clearly bypass the limitations of this contribution. Of course also the density parameters, Φ f and Φ m , have their impact on the shape of the radial distribution function. Interestingly, if one compares the g(r) of systems with (approximately) the same total packing fraction, Φ tot , and given βκ the respective curves differ only marginally in their shape; thus we conclude that the shape of g(r) primarily depends on the total packing of the system and that the parition of this quantity into matrix and fluid densities plays only a minor role. Possibly similar conclusions were also available in the preceding study of spherical fluid particles in a matrix of hard spheres [32][33][34], but have not been investigated more closely.

IV. CONCLUSIONS
We have investigated in extensive Monte Carlo simulations the properties of impenetrable, deformable, fluid particles that are immersed into a matrix formed by immobile, impenetrable spherical particles. Using the concept of deformable particles, put forward by Batista and Miller [44,45], the mobile constituents of the system can modify their shape within the class of ellipsoids of revolution. The energy change related to this deformation is fed into suitably adapted Monte Carlo selection rules; thus we apply in our simulations -apart from the translational and orientational moves -also shape deformation moves for the fluid particles. Overlap between the (fluid and matrix) particles was detected by using a criterion derived by Vieillard-Baron [46] for elliptic particles; we note that from the computational point of view this criterion is computationally considerably more expensive than the one for simple spherical particles. The above mentioned energy change related to the shape deformation is a function of the aspect ratio x of the particles; its strength is measured in terms of the deformability βκ: for βκ → ∞ we recover undeformable hard spheres, while small βκ values specify easily deformable particles.
Our investigations are based on the quenched-annealed concept [7][8][9][10][11] where physical 21 quantities are evaluated from the particle positions and orientations via a double averaging procedure: (i) in a first step a trace is taken over the degrees of freedom of the fluid particles for a given matrix configuration; (ii) in a second step observables are averaged over a limited number of different, but equivalent matrix configurations; the high numerical cost imposed by the overlap criterion of elliptic particles forced us to limit the second averaging process to five independent matrix configurations per state point investigated. In total, the system is characterized by the packing fractions of the mobile (fluid) and the immobile matrix particles (Φ f and Φ m , respectively) and the deformability parameter βκ.
In an effort to extract as much information as possible we have defined in the parameter space spanned by Φ f and Φ m four pathways which were also used in a previous investigation of spherical fluid particles confined in a matrix of spherical particles [32][33][34]. The pathways are characterized by keeping either of the aforementioned densities constant, while varying the other density over a representative range. For the specification of these ranges we dwelled on the previous investigation which provided some insight into the density range where the system remains in a disordered -i.e., either liquid or at most glassy -state.
Our investigations and our analysis are based on the radial distribution function, g(r), the distribution of the aspect ratio (ARD), the spatial mean squared displacement (MSD), δr 2 (t), and its effective exponent, z(t).
The shape of the radial distribution function is primarily influenced by the deformability of the particles: in systems of strongly deformable, mobile particles the main peak in g(r) is a smooth function of the distance and its onset can be located at distances as 0.7 σ.
In contrast, hardly deformable particles have harsh, hard-particle type main peaks which are located essentially at σ. For a given value of βκ the shape of g(r) is imposed by the total packing of the system. The ARD is a strongly peaked function for high values of deformability and extends over a growing x-range as βκ decreases; only for strongly deformable particles an impact of the packing of the system is observable: for very low densities the ARD can extend over a remarkably broad range of x-values.