Ordering, clustering, and wetting of hard rods in extreme confinement

The effect of out-of-plane orientational and positional fluctuations is examined in the phase behavior of hard spherocylinders confined between two parallel walls. The stability of isotropic, nematic, and solid phases is studied for aspect ratios (κ = 1 + L/σ , where σ and L are the diameter and length of the cylinder) of 8, 10, and 16, while the width of the slitlike pore, H , is set in the quasi-two-dimensional regime, σ < H 2σ . Using replica exchange Monte Carlo (REMC) simulations and Onsager theory we provide evidence for the stabilization of the nematic order with increasing pore width. The minimum surface density necessary to exhibit nematic order remains nearly the same with increasing H , while the upper bound of nematic order is postponed to higher densities in detriment of the single-layer-solid phase. We prove that a drying-out effect takes place in the pore as the density increases linearly from the wall to the middle of the pore in the ideal gas limit. This behavior is kept for the isotropic phase and partially in the nematic one, while a wetting behavior is observed at the walls mostly for the solid phases. For κ = 8, we have also observed a compressibility signal of a transition between an isotropic fluid and a cluster-like fluid. This cluster phase is destabilized by increasing H . Finally, we build a master diagram, where the drying-wetting, solid-solid, and close packing density curves are independent from κ .


I. INTRODUCTION
The rich phase behavior of hard body fluids depends very strongly on the shape of the particles [1] and the confining dimension [2,3]. It is now well understood that the anisotropic shape of particles is responsible for the formation of mesophases such as the the nematic, smectic, and columnar. Similarly, confinement gives rise to very rich phase behavior. For example, depending on the interaction between the confining walls and the particles, several phenomena may arise such the wetting, evaporation, layering, and different types of anchoring [4][5][6][7][8]. Probably, the most important consequences of strong confinement are related to the fact that the nature of orientational ordering of anisotropic particles is quite different in two and three dimensions. The three-dimensional (3D) nematic phase is accompanied by long range orientational order, while the two-dimensional (2D) one is quasi-long ranged and the orientational correlation decays algebraically. In addition to this, in case of particles with sufficiently large anisotropy, the 3D isotropic-nematic (I-N) phase transition is usually first order, while the 2D nematic transforms continuously to the isotropic phase through the Kosterlitz-Thouless disclination unbinding mechanism with decreasing density [9,10]. The strange nature of 2D nematic phase was explored with several simulation studies [11][12][13][14][15][16][17] and even with experiments [18,19]. Note that similar differences are present in the properties of 2D and 3D solid phases, too [20][21][22][23]. * godriozo@azc.uam.mx The simplest way to investigate the crossover between 2D and 3D phase behavior of hard body fluids is to confine the particles between two parallel hard plates, where setting the wall-to-wall distance (pore width) from infinity to the shortest dimension of the particle allows to interpolate between 3D and 2D bulk properties. That is, by tuning the pore width it is possible to examine the expected intermediate behavior between 2D and 3D orientational and positional ordering of rod-like particles. It is also possible to do experiments with hard ellipsoids (or spherocylinders), where the confined hard ellipsoids corresponds to the 2D bulk system of hard ellipses in the case of most extreme confinement [24]. Comparing the strictly 2D and 3D I-N phase transitions, it is observed that the minimum shape anisotropy required to stabilize the nematic ordering is very different. In the case of hard ellipsoids, confinement stabilizes slightly the nematic order, which is neither 2D nor 3D nematic, because the aspect ratio (the ratio of the lengths of major and minor axes) must exceed 2.8 in 3D bulk [25,26], while the minimal aspect ratio is only 2.4 in 2D bulk [27]. This slight stabilization effect of hard flat walls on the nematic ordering of the ellipsoids was confirmed by Monte Carlo simulations [28,29]. The opposite trend happens in the case of hard spherocylinders, i. e. confinement destabilizes the nematic ordering, since the minimal aspect ratio κ (the total length divided by the diameter) is 4.8 in 3D [30,31], while it is around 8 in 2D. A possible explanation for this behavior is that the cylindrical part of the particle behaves as a hard rectangle in 2D confinement, which is known to favor the formation of clusters and, depending on the shape of the particle tip, even the formation of a tetratic phase is fea-sible when the anisotropy is moderate [32][33][34]. The effect of central and terminal shape parts of rod-like particles was examined in a quasi-2D granular system, where the monolayer of hard rods was vibrated vertically [35]. It was found that the cylindrical shape of the central part of the particles is responsible for the tetratic correlation, while rounded or tapered tips favor nematic order [35]. Monte Carlo simulation studies also confirm these findings, given that the tetratic phase was found in the system of weakly anisotropic hard rectangles [36,37], while only short ranged tetratic correlations were observed in the case of hard discorectangles [32]. Therefore, it is logical to consider the tetratic correlation in strictly confined hard spherocylinders (hard discorectangles) as a key factor for the destabilisation of the nematic phase.
It must be also taken into account that the properties of 2D and 3D nematic phases are qualitatively different. This can be seen in the system of hard spherocylinders placed in a slit-like pore, where the hard flat surfaces induce a planar adsorption and a surface phase transition between orientationally disordered and ordered surface layers (2D I-N ordering) with increasing density in the pore. This surface induced 2D I-N transition is not related to the capillary I-N transition of 3D nematisation, because the 3D I-N transition terminates at a critical density larger than the 2D I-N transition density [38][39][40]. Monte Carlo simulations found that the 3D nematisation can be strengthen with narrowing pores in such an extent that the capillary isotropic phase vanishes at pore widths smaller than about the double of the rod length so that only the surface induced nematic phase survives in very narrow pores [41]. The 2D I-N transition was studied for the fluid of very long hard spherocylinders, which are placed between two planar hard walls [41][42][43]. It was found that the effect of adding out-of-plane positional and orientational freedom destabilizes the 2D nematic ordering, which means that more and more particles are needed to form the 2D nematic phase with increasing wall-to-wall distance in the pore [44]. This is due to the fact that the out-of-plane orientation freedom decreases the effective in-plane shape anisotropy of the particles, while the positional out-of-plane ordering increases the available room for particles. From a 2D point of view, the confined rods can be considered as a polydisperse mixture of 2D hard discorectangles, where the polydispersity takes place both in length, due to the increase of orientational freedom, and diameter, due to the increase of out-of-plane positional freedom. The effect of polydispersity can be stronger for shorter rods, where the nature of 2D I-N transition changes from continuous to first order [45].
In this work, we study the phase behavior of hard spherocylinders which are placed between two parallel walls using replica exchange Monte Carlo simulations and Onsager theory. The length, L, and the diameter, σ, of the cylindrical part of the particle is chosen such that the aspect ratio, k = 1 + L/σ, is fixed to 8, 10, and 16. The wall-to-wall distance, H, is restricted to stay between σ and 2σ, σ < H < 2σ, which makes the system be in the quasi-2D region, where two independent fluid layers cannot evolve in the slit. This means that layers adsorbed on opposite walls directly interact with each other. We pay special attention to the adsorption properties, the cluster formation, and the stability of the ordered phases. We show that almost the same number of particles per unit area, ρ, is enough for the formation of nematic order with H = σ and H = 2σ, even though the available volume for the particles is double for H = 2σ. In addition to this, wetting of the walls occurs at lower ρ values with increasing the pore width H, whereas the opposite tendency can be seen for the single-layer-solid (SLS) formation curve. All these results are summarized in a phase diagram, which share similar features for all studied κ values. Indeed, the wetting lines, the maximal packing density lines, and the two-layer-solid (TLS) formation lines approximately collapse to define master curves.

II. SIMULATION DETAILS
We have performed replica exchange Monte Carlo (REMC) simulations of confined hard spherocylinders. This technique was developed to enhance sampling from equilibrium, specially in the cases where the free energy landscape is uneven [46,47]. It is based on the definition of an extended ensemble, Q ext = i Q i , Q i being a partition function of ensemble i. In our case, we make use of an expansion of the isobaric ensemble, Q ext = i Q(N, P i , T ), where N , T , and P i are the number of particles, temperature, and the three dimensional pressure of ensemble i, respectively [48]. Here, N and T are fixed for all ensembles defining Q ext . This expansion is useful to deal with systems of hard particles [49][50][51].
As mentioned, spherocylinders are hard particles made up of two hemispheres of diameter σ connected by cylinders of length L. For these particles, the aspect ratio is then given by κ = L/σ + 1. The overlapping algorithm for spherocylinders consists on a simple method, where the distance between the nearest points of two line-segments with length L must be larger than σ [52]. In addition, the overlapping of spherocylinders with walls is simply achieved by comparing the distances between the center of the spheres at each end of the particle and the corresponding wall, which cannot be less than the hemisphere radius, σ/2. Our system is then built by placing N spherocylinders between two parallel hard walls separated a distance H and perpendicular to the z axis. We build two different starting arrangements. On the one hand, we consider N r different loose initial configurations, where particles positions and orientations are taken at random while avoiding overlaps. On the other hand, we also build N r identical compact initial configurations, where particles are packed parallel to the walls in order to yield the densest packing, see the Supplemental Material (SM). That is, by stacking side-by-side particles alternating the z-position to get in contact with each wall while forming columns, which are then assembled by making the hemi-spheres of the one column enter into the holes of the next one. In the first case, the system evolves by increasing density, whereas in the second case the system of replicas follows an expansion.
We have applied periodic boundary conditions only in the x and y axes, given that the z axis is bounded by the walls. We have also employed Verlet lists to gain efficiency. The values of βP i follow a geometric progression, with a minimum pressure of βP min = 0.02σ −3 and a maxi- Each of the N r replicas performs the sampling of the isobaric ensemble at a fixed P i value, where conventional displacement, rotation, and volume change trials are carried out. This is implemented in Graphic Processing Units (GPUs), being each replica handled by a single core, allowing for the definition of a large number of replicas. Next, information is gathered by the CPU where replica swaps are tried with acceptance probabil- and V j are the pressures and volumes of replicas i and j, respectively. We allow swap trials of replicas sampling from ensembles i and j with |i − j| ≥ 1, in such a way that we obtain an average acceptance around 30%. This is done as explained in detail elsewhere [53].
In all cases we run the simulations long enough to ensure the system or replicas reaches a stationary state, presumably equilibrium. After achieving this condition, we perform several averages over the system of replicas, including probability density functions and structural properties. For most runs we have defined a relatively large number of replicas, N r = 230, and N = 128, 160 and 144 for κ = 8, 10 and 16, respectively. Note that this technique takes advantage of density fluctuations that increase with decreasing the system size. Consequently, it does not provide an efficient way to study long-range correlations among particles, but yields a detailed equation of state, βP (ρ), and isothermal compressibility, χ(ρ). However, in a few cases we have performed runs with N r = 80 and N 1400 to observe quasi-long range correlations.

III. RESULTS
In the case of H → σ, the confined system of spherocylinders by parallel plates should produce the same results as discorectangles. This last case is studied in reference [32] for several κ values and in particular for κ = 16. We consider this case as a reference to check the correctness of our implementation. For this purpose, we perform simulations with κ = 16 and H = σ, by compressing loose random configurations and expanding tightly packed structures. The comparison of our outcomes with those given in reference [32] is shown in panel a) of Fig. 1, where the 2D pressure, βP H, is plotted against the 2D number density, ρ = N H/V . Note that each of our curves is constituted by 230 points, and this is why we employ lines to represent them. Compression and expansion runs yield the same stationary curves, which in addition coincide with the data reported by Bates and Frenkel [32]. Hence, we can safely conclude that our implementation is correct. In addition, the isothermal compressibility, as obtained from density fluctuations through χ = N ( ρ 2 − ρ 2 )/ ρ 2 , also produce identical curves for compression and expansion runs. The same occurs for the nematic order parameter where φ i is the angle between an arbitrary fixed axis on the xy plane and the projection of the axis of the i th particle on this plane. Notwithstanding, we have observed marked differences of this parameter when comparing runs starting from loose random configurations and compact ordered structures during a long equilibration time, despite the fact that βP (ρ) and χ(ρ) were nearly the same at much shorter times. The same was observed for the tetratic order parameter, S T , which is defined by replacing 2φ by 4φ in Eq. (1) (not shown). For 0.038 < ρ < 0.045, the former starting condition produces some configurations with lower S N and higher S T values than the second one, which yields arrangements where all particles are clearly aligned with a well defined director. The fact of having different configurations during long simulation times signals that equilibrium is hard to achieve. Indeed, the compression curve for S N and S T were observed to evolve very slowly with REMC cycles, approaching the expansion one. From here on, we only show expansion curves, which are observed to produce equilibrium with a smaller computing effort. The isothermal compressibility curve, χ(ρ), shown in the panel b) of Fig. 1 shows two peaks. The first one is linked to a sudden increase of the nematic order parameter, S N , which is a signature of the I-N transition. In addition, S N > S T for all ρ values. Snapshots corresponding to the isotropic and nematic fluids are depicted at the left of Fig. 1. Furthermore, this first χ peak, located at ρ IN = 0.021 ± 0.001, is closely preceded by the inflection point of S N (ρ). We observe that the density distribution is always unimodal for all pressure values around this point signaling a high order transition (not shown), which is consistent with a Kosterlitz-Thouless disclination unbinding mechanism [9] and agrees with strictly 2D results [32]. Orientational correlations in the 2D nematic phase decays algebraically, and so S N also diminishes with increasing the system size. Nonetheless, we expect the shape of curve S N (ρ) to hold with increasing N and so, the location of its inflexion point should remain invariant for finite systems. The second peak of χ corresponds to the nematic-solid transition reported in the phase diagram given elsewhere [32] and occurs at ρ NS = 0.041 ± 0.001. As can be seen, it also coincides with a smaller increase of S N (ρ), where the particles further align to form the solid-like structure. Here, the appearance of columns and rows is observed. Recall that in two dimensions there are no true crystals, given that thermal fluctuations destroys long range order [54]. According to the unimodal density distributions we have obtained, although stronger than the I-N transition, the nematic-solid transition is not first order.
For the same κ = 16, fixing H to 1.8σ produces the results shown in Fig. 2. It is observed a somewhat uneven but smooth increasing βP H as a function of ρ, which translates into a χ(ρ) showing the three maxima signaled by the vertical dashed lines. With increasing pressure, the first one corresponds to the I-N transition, the second one to the nematic-solid transition, and the third one, a very wide peak occurring at much larger ρ values, is related to a layering process. Note that the isotropic-nematic location, ρ IN , is close to the inflexion point of S N (ρ), which is practically the same as for the case H ≈ σ. In addition, ρ NS is increased a bit and it continues to match a small increase of S N (ρ). In Fig. 2 c) we show, together with S N , the layering order parameter which is given by increasing rate of S L (ρ) is small, which yields a not so well defined inflexion point. This point seems to coincide with the χ(ρ) maxima and we take it as the transition between a SLS and a TLS. The three snapshots inserted in Fig. 2 corresponds to the dry SLS (SLSd), wet SLS (SLSw), and wet TLS (TLSw), from top to bottom. This last structure shows neighboring particles belonging to the same column alternating positions in the z direction, which leads to a much larger packing density (see the SM for a comparison between the snapshots of 1LSw and the 2LSw for N =1400).
As mentioned above, a simple criteria for discriminating wetting from drying could be given in terms of S L . Its definition relies on the density profile of the particles centres, ρ(z), or equivalently, on the positional distribution function defined as f (z) = ρ(z)/ρ, which fulfills −H/2 f (z)dz = 1. This positional distribution function is shown in Fig. 3 for κ = 16 and H = 2σ. In this plot we employ different colors to distinguish among the curves belonging to different phases. We also insert as a black line the f (z) corresponding to the ideal gas limit. In this case, the exact calculation yields for H < κσ. This expression is derived in the first section of the SM. We can observe that this limiting formula agrees perfectly with our simulation results for ρ < 0.02. At these densities a drying-out effect takes place because the particles move away from the walls and accumulate in the middle of the pore. This drying-out effect is due to the fact that the density is low enough and the particles can maximize the number of microstates with the gain coming from the orientational freedom when they are located in the middle of the pore. A significant departure from the ideal gas limit is observed for the nematic phase, where the midplane maximum diminishes its height and small peaks emerge close to the walls. For a large enough pressure, the profile looks like the highlighted blue curve, which corresponds to S L ∼ 0 and has the peaks at the walls with a similar height than that at the midplane. Thus, this curve can be considered as a dividing line between the drying and wetting phenomena of the walls for H = 1.8σ. Conversely, the particles wet the walls in the solid phase. For a low pressure solid, wetting is not strong enough to produce a couple of well defined layers at the walls because several particles are still located close to the midplane. For a large enough pressure particles are totally excluded from the midplane and peaks at the walls are fully developed. In between these two cases, peaks at the walls continuously grow with increasing pressure, while the depletion at the middle of the pore becomes stronger. This yields a smooth grow of S L (ρ), which signals a broad solid-solid transition. This is due to the fact that the SLS-TLS transition imply not only the formation of two layers but also their arrangement to yield a structure where particles alternate their z positions.
Our results deviate substantially from previous simulation results, which were performed in much wider pores. It was observed that the particles wet the walls and a thin nematic film forms in the vicinity of the walls at such densities, where the middle of the pore is still isotropic [41,42]. The change of the equation of state, the compressibility, and the order parameter curves with increasing H can be seen in Fig. 4. At low ρ values, all curves collapse and produce the first peak of χ(ρ) with practically the same ρ IN . As mentioned, this peak is correlated to a sudden increase of S N . It should be noted that for the same 2D number density ρ and for increasing H, S N slightly decreases. This is a result of the gain in possible polar angle values, which decreases the 2D shape anisotropy (the effective kappa in the xy plane). The nematic-solid transition, the one observed at ρ ∼ 0.04, shifts to the right with increasing H. The main consequence of these results is that the out-of-plane orientational and positional freedom can support the formation of 2D nematic ordering. Thus, the density range for the nematic phase widens. It is surprising that practically the same amount of particles is enough to induce nematic ordering when we double the pore width and the available room for the particles. At even larger densities, the curves for low and high H behave very differently. For H < 1.6σ the shape of βP (ρ) is nearly kept, strongly contrasting with cases having H > 1.6σ, where the systems yield much larger densities. This increase in density is accompanied with a wide peak of χ(ρ), as discussed for H = 1.8σ, which in turn is linked to an increase of S L , indicating that the particles start locating close to the walls. Note that the density at which wetting commences, ρ W , shifts to the left with increasing H. Nonetheless, except for H = 2σ, wetting always starts once the solid phase is formed. Note also that for H > 1.8σ the location of the χ(ρ) peak signaling the solid formation is not clear. Hence, in this case we have performed simulations with N 1400 and N r = 80 to determine, from the decayment of the radial distribution peaks, the transition density. We have followed the same procedure for all κ values to determine all fluid-solid transitions occurring above the wetting line.
For H σ, decreasing κ leads to the disappearance of the nematic fluid phase. According to Bates and Frenkel [32] this happens close to κ = 8. In this case, we observe S N curves reaching values close to 0.5 and some snapshots where a clear director can be defined, as shown by the top right snapshot of Fig. 5. Nonetheless, at the inflexion point of S N (ρ), neither the χ(ρ) curve nor its derivative,χ(ρ), show any sign of an I-N transition. In addition, in this case we have also performed simulations with N ≈ 1400 showing much lower S N values and a fast decay of orientational correlations, confirming the absence of a nematic fluid. So, κ = 8 seems to be just below the κ limit for obtaining a nematic fluid phase. Note also that the nematic phase is frustrated by the appearance of clusters forming a mosaic-like arrangement, which makes S N (ρ) to decrease while S T (ρ) increases up to the point of yielding S T (ρ) > S N (ρ). Recall that these clusters appear even by starting from tight compact parallel configurations and so, one can safely expect them to belong to equilibrium. Furthermore, the appearance of these clusters is accompanied by a tiny peak of χ(ρ), probably more noticeable by the zero value ofχ(ρ), suggesting the existence of an isotropic to tetratic transition, as the one formed by hard rectangles [37], although in this last case S N (ρ) remains close to zero for all ρ. We should mention that Bates and Frenkel described this phase as globally isotropic, with strong local positional and orientational correlations between particles [32]. They have indeed ruled out the existence of a tetratic phase for κ = 6 and ρ = 0.125 by studying the decay of the orientational order with distance for systems with N = 1000. We have confirmed this with our N ≈ 1400 simulations for κ = 8 and H σ. So, in view of these results, it is difficult to name this phase otherwise than globally isotropic with strong tetratic correlations. Nonetheless, given the existence of this tiny χ peak, which persists with increasing H up to H = 1.7σ as shown in Fig. 6 b), let us call it cluster phase. It should also be mentioned that tetratic correlations are also found for similar systems [35,55], where particle tips are shown to have an important impact on the presence of the tetratic order. We speculate that there exist a system with κ 8 which forms a nematic phase at H = 2σ and produces clusters at H = σ.
Results for κ = 8 and σ ≤ H ≤ 2σ are shown in Fig. 6. As for all κ, the curves collapse for low ρ values, as the system approaches the ideal gas limit, and they shift to the right with increasing H for large ρ values. In addition, the χ peak signaling the isotropic-cluster transition shifts to the right and decreases its height, vanishing at H = 1.8σ. This is accompanied by a smoothing of the S N (ρ) drop as well as the disappearance of the well of the difference of the nematic and tetratic curves, S N (ρ)−S T (ρ). General phase diagram of hard spherocylinders confined by two walls. ρ * = ρ/ρcp, where ρcp is the 2D close packing density for H = σ. Inaccessible densities are shown by the hatched area. Blue, green, and red lines correspond to κ = 8, 10, and 16, respectively. Solid squares, circles, and diamonds correspond to SLS-TLS transition, nematic-solid transition, and I-N transition, respectively. Open diamonds correspond to a probable clustering transition for κ = 8. Triangles correspond to wetting line, i. e. the density location at which SL = 0. This line divides the diagram into wetting and drying regions, which are colored with bluish and reddish tones, respectively. Labels I, N, SLS, and TLS correspond to isotropic, nematic, single-layer-solid, and two-layer-solid, respectively. The final letter in the labels, w and d, make reference to wetting and drying, respectively. Blue labels Iw and Id correspond only to the κ = 8 case.
shows an I-N transition while producing some cluster-like configurations.
The outcomes of the Onsager theory are compared to data from simulations in Fig. 8 for κ = 16 and 10. The implementation of the Onsager theory is given in in section II of the SM. It is observed that the general trends of the equation of states, βP (ρ)H, and the nematic parameter curves, S N (ρ), are qualitatively captured by the theory. In addition, for increasing H it captures the shift of βP (ρ)H toward smaller values, and for increasing κ the shift of S N (ρ) to reach lower densities. Also, the theory predicts the independence of the I-N transition from H, which is practically confirmed by simulations. It seems that the theory improves the prediction of ρ IN for increasing κ.
Finally, we build a general phase diagram that depicts ρ * = ρ/ρ cp as a function of H, ρ cp being the close packing density for H = σ. The original phase diagrams as a function of ρ and η are given in the SM. This general phase diagram is given in Fig. 8. We hatch the region of inaccessible densities and include the I-N transitions, the nematic-solid transitions, and the SLS-TLS transitions. Note that the maximal packing density curves obtained for the different κ values collapse to define a master curve, and this fact allows us to define a common inaccessible region. The same happens with the obtained wetting lines, which collapse to define a single curve. This let us define the same wetting and drying areas for all κ values. These are colored with pale red for drying and light cyan for wetting. A similar collapse is also observed for the two layer solid formation curves. As a consequence of the wetting and drying regions, we add an extra letter to the phases labeled in Fig. 8. For instance, the Nd and Nw labels for κ = 10 and 16 stand for drying and wetting nematic. No transition is observed between Nd and Nw as only a simple change in S L (ρ) occurs with increasing density. Similarly, we can identify an isotropic drying (Id) and an isotropic wetting (Iw) phase for κ = 8. The SLS phase can also be drying or wetting. In the case of low H, the drying is preferred because the particles are located mainly in the center of pore to increase the orientational entropy. Note that the wetting region vanishes completely as H approaches σ. The SLSw appears at lower densities than the TLS, which is always wetting. Hence, we observe the following phase sequence with increasing density for H = 1.6σ and κ = 10 or 16: 1) isotropic, 2) nematic, 3) SLSd, 4) SLSw and 5) TLSw. It is worth noting that the SLS density range decreases its size with increasing H, vanishing at H = 2σ. At this point, the nematic phase changes directly into a TLS.
Conversely to the close packing, wetting, and solid-solid transition curves, the I-N and the nematic-solid transitions shift toward lower densities with increasing κ. In other words, increasing κ produces a relative easiness to form ordered phases. We can see in the master diagram that the density range of the nematic widens at the expenses of the SLS with increasing pore width, while the density range of the isotropic phase is nearly constant, in good agreement with theoretical predictions (see Eq. (18) of SM). This means that the same number of particles per unit of area is needed to yield a nematic phase, irrespective of H. This is a quite unexpected result given the obvious increase of accessible volume with H. Finally, the kink in the compressibilty we have obtained for κ = 8, associated a possible isotropic-cluster transition, disappears for H > 1.7σ. We have observed that increasing H disfavors the tetratic correlations.

IV. CONCLUSIONS
We have built the phase diagram of strongly confined spherocylinders with anisotropy 8, 10, and 16 in slit-like pores. We have constrained the particles to the pore width range σ < H < 2σ. For this purpose, we have performed replica exchange Monte Carlo simulations, exact calculations in the ideal gas limit, and bifurcation analysis of the Onsager-theory. The stability regions of the isotropic, nematic, and solid phases have been determined to see the effect of varying the shape anisotropy and the pore width. In addition, the upper bound of the solid phase is bordered with the maximimal density of the close packed solid structure.
The three different phase diagrams corresponding to the different particle's shape anisotropies, κ, share very similar features (see the last section of the SM). This allows to build a general phase diagram where some curves seems to be independent of κ. In particular, the maximal accessible packing density, the lines splitting the drying and wetting regions, and the transitions toward a TLS, are very close to define a master curve. Conversely, the I-N and the nematic-SLS transitions are clearly κ dependent.
Our results also show that the isotropic-nematic curves do not depend on H, i. e. the same number of particles is needed to induce nematic ordering at H = σ and H = 2σ. This is against the naive linear expectation that doubling the pore width should be accompanied with doubling the number of particles to get the a similar behavior. Moreover, the out-of-plane orientational and positional fluctuations are believed to work against the orientational ordering as the 2D projected area of the hard spherocylinder is less anisotropic. Conversely to the I-N density dependence with H, the nematic-solid transition is postponed to higher densities as positional and orientational fluctuations are favored. Therefore, the stability region of the nematic ordering is enhanced by widening the pore. This behavior contrasts with the Monte Carlo results reported for L = 320σ, where the I-N transition density increases with the pore width for H > 2σ (note that the σ < H < 2σ region was not examined [44]).
Another important aspect of our work focuses on the drying or wetting behavior of the particles at the walls. In the ideal gas limit, we show that the local density increases linearly from the wall to the middle of the pore. This behavior results in a drying-out effect at the walls in the isotropic phase, whereas wetting of the walls is observed mostly for the solid phase. The nematic phase also exhibits a drying-out effect except close to H = 2σ, where wetting emerges at the walls.
Finally, we have observed for κ = 8 a compressibility signal, compatible with a weak higher order phase transition between an isotropic fluid and a cluster-like fluid, neither of them showing quasi-long range orientational order of any type. This cluster phase is destabilized by increasing H. We speculate the same destabilization may appear when considering particles with different degree of sharpening of their tips.
We would like to emphasize that our results show that the experimental detection of 2D nematic ordering can be easier in wider pores because the system can be less dense and the cluster formation is not entropically favored. It is also feasible that the same stabilization of the nematic phase by widening the pore could occur for particles without rounded ends, such as cylinders.