Enhanced two-dimensional nematic order in slit-like pores

The effect of out-of-plane positional freedom is examined on the stability of two-dimensional (2D) nematic order of hard non-spherical particles using the second virial density-functional theory. The particles are allowed to move and rotate freely in the plane of confining walls and can move between the two parallel walls. The wall-to-wall distance (H) is varied between the strictly 2D and the two-layer forming cases, i.e. σ < H < 2σ, where σ is the particle’s shortest length. As expected, we observe that more and more particles are required for the formation of 2D nematics with increasing H when the rod-like particles are hard ellipsoids. Surprisingly, we found that the opposite tendency is observed in the case of hard cylinders, i.e. fewer and fewer particles are needed to stabilize the nematic order with increasing H. This paradox can be understood by projecting the three-dimensional system into a 2D mixture of particles having position-dependent aspect ratios and molecular areas. However, the complex phase behaviour found for plate-like cylindrical particles with increasing H cannot be explained in terms of the same simple geometrical arguments. Video Abstract: Enhanced two-dimensional nematic order in slit-like pores


Introduction
The lack of long-range positional and orientational order makes the two-dimensional (2D) systems peculiarly interesting from the academic point of view [1]. The observed short and quasi-long range orders are due to the fact that the long-wavelength fluctuations destroy the long-range orientational and positional order with the dimensional restriction. Despite the presence of these strong fluctuations, 2D systems can exhibit a very rich phase behaviour, structural changes, Kosterlitz-Thouless (KT) continuous transitions, and first-order phase transitions [2][3][4][5]. For example, the freezing of hard disks is accompanied by a partially ordered hexatic phase [6], that of hard squares by a tetratic phase [7], and a solid-solid transition occurs between orientationally disordered and ordered solid structures [8]. Even two intermediate mesophases can be sandwiched between the isotropic fluid and solid phases of bulk superdisks, whose shape is between the disk and the square [9]. However, the richness of phases can be enhanced by introducing out-of-plane positional fluctuations by making the system quasi-two-dimensional (q2D) [10][11][12][13][14][15]. The orientational and positional alignment is particularly important in nanodevice applications such as the nanowire arrays and nanotransistors, where the aim is to fabricate very dense and ordered monolayers with nanotubes and nanowires [16][17][18][19][20].
The main features of the 2D nematic phase (nematics) of thermalised systems are the short range positional order, the lack of long-range orientational order, and the algebraically decaying orientational correlations [15]. This situation is also common in the nematic phase of active matter [21,22], even though there are some indications of long-range nematic order in the suspensions of swimming filamentous bacteria in very thin layers [23]. Note that the exclusion of true long-range orientational order in 2D systems can be proved rigorously only for special pair potentials, where the orientational and positional freedoms can be decoupled from each other [24]. Unfortunately, this decoupling is not possible in the case of rod-like and plate-like particles such as cylinders and ellipsoids. To elucidate this issue, several simulation Hard cylinders with lengths L and diameters D located between two parallel flat planes. The particles are allowed to move freely between the walls, but they are restricted orientationally as they can rotate freely only in the x-y plane (planar order). The width of the pore is chosen to maintain the q2D character of the pore, i.e. D H 2D. studies were performed with 2D rod-like particles confirming the lack of true long-range orientational order [10,11,25,26].
The isotropic-nematic (IN) transition is proved to be continuous through a KT disclination unbinding mechanism for several shapes [10,25,26]. The widely used density functional theories predict a second-order phase transition and long-range nematic order [15]. In addition to this, these mean-field theories underestimate the IN transition density and overestimate the extent of orientational order in the nematic phase. Contrasting with these failures, the mean-field theories reproduce quite accurately the Monte Carlo simulation results obtained for the equation of state in both isotropic and nematic phases [27]. Moreover, they correctly describe the shape dependence of IN transition density in nonspherical hard body fluids [15].
The experimental setup of strictly 2D systems is still a big challenge. The possible ways to restrict the orientational and positional freedoms of the particles in two dimensions are the followings: (1) confinement into nano-size slit-like pores, (2) adsorption at smooth surfaces, and (3) confinement at air-liquid interfaces [28][29][30][31][32][33][34][35][36][37][38]. With these techniques, it is possible to create q2D systems where the out-of-plane positional and orientational freedoms are not completely vanished. By using non-spherical colloidal particles, 2D nematic phases are found in several experimental studies, and the order of IN phase transitions is continuous or first-order [28,36,37]. To get closer to the experimental setup and to understand the differences between the experimental and theoretical results, it is worth studying q2D systems where the effects of out-of-plane positional and orientational fluctuations are considered. Simulations show that the two-step melting scenario of the hard-sphere monolayer is not affected by the out-of-plane positional fluctuations when the perpendicular particle motion is within half the sphere's diameter [39]. They also show that the 2D IN transition of hard spherocylinders preserve the KT type nature of strictly 2D hard body systems when slightly increasing the pore size [40]. It was also observed that the 2D transition density increases with the thickness of the pore for both flexible and rigid rods, i.e. widening the pore destabilizes the nematic phase [40,41]. Indeed, there are several simulations and theoretical works devoted to studying the capillary nematisation, the smectic order, the layering transitions, the surface ordering transitions, and the surface-induced biaxiality for different shapes and anisotropies in slit-like geometries [42][43][44][45][46][47][48][49][50][51][52][53][54][55]. Nonetheless, a stability analysis of the 2D nematic phase in the presence of out-of-plane positional and orientational fluctuations is still missing for several molecular shapes and sizes.
This work aims to fill this gap for cylinders and ellipsoids. It is generally accepted that the number of particles needed for stabilizing the nematic phase must increase with the widening of the pore as the available room for the particles increases [40,41]. Here we show that this expectation is valid for hard ellipsoids, where the out-of-plane positional freedom weakens the orientational coupling between the freely moving particles. Surprisingly, the rod-like (prolate) hard cylinders exhibit the opposite behaviour as the orientational coupling becomes stronger with a gradual increment of the out-of-plane positional freedom. Thus, it turns out that fewer and fewer particles are needed for the stabilization of nematic order. This unexpected result is also found in the case of plate-like (oblate) cylinders. We produce these results by placing rod-like particles between two parallel flat walls, where the extremely small wall-to-wall distance enforces the particles to align parallel to the surface of the walls (planar order). In the case of plate-like cylinders, the same type of order, which is often called edge-on, can be achieved between two parallel walls when an extra external electric field is applied in the direction of the surface's normal to overcome the entropy-driven face-on order effect of the walls. With this technique, it was possible to set up a q2D fluid of plate-like particles and study the 2D isotropic-nematic transition of hard rectangles, since the interacting cross-sections of the plates have this shape [56]. We use the second virial density functional theory to locate the IN transition and to measure the extent of the nematic order as a function of the pore width for a given 2D density. We note that our present study is closely related to a previous one on confined hard spherocylinders, where both the orientational and positional ordering are examined using replica-exchange Monte Carlo simulations [57]. An interesting result of this work is that the number of spherocylinders needed to stabilize the nematic order does not depend on the pore width in the q2D regime. Therefore, in this regard, hard spherocylinders behave between cylinders and ellipsoids.

Models and theory
We examine the stability of 2D nematic order of monodisperse hard cylinders and that of hard ellipsoids in very narrow slit-like pores, where the distance between the two parallel flat walls (H) is increased gradually up to the limit at which the system abruptly loses its q2D character. To satisfy the q2D fluid condition, the range of the pore width is restricted to D H 2D in the case of hard cylinders with lengths L and diameters D (see figure 1), while the confined hard ellipsoids with σ major and σ ⊥ minor lengths fulfill σ ⊥ H 2σ ⊥ . At the lower limit (H = D for cylinders and H = σ ⊥ for ellipsoids) the confined system turns into the strict 2D case of hard rectangles and hard ellipses, respectively. The out-of-plane translational and orientational freedom can be switched on continuously by increasing H. By doing so, we increase the available room for all particles while the number of particles is kept fixed. Therefore, one can naively think that the particles can rotate more freely while the system becomes more disordered with the widening of the pore.
In the vicinity of the hard walls, the out-of-plane orientational freedom is suppressed, which results in a drying-out effect of the walls at low densities. It can be shown that the density linearly increases from the wall to the middle of the pore in the low-density limit, which is the case of non-interacting particles. In the presence of particle-particle interactions, the density profile is determined as the result of the competition between orientational and packing entropy. With increasing density, the excluded volume interactions between the particles become more and more dominant, which yields their alignment and pushes them towards the walls to maximize the available room. Therefore, at some point, the drying-out effect is replaced by wetting of the walls [57]. These results highlight the fact that the suppression of out-of-plane orientational freedom cannot be applied at low densities but only at relatively high densities when the particles form two interacting fluid layers near the walls. These two layers are coupled, i.e. they cannot move independently from each other for D H 2D (or σ ⊥ H 2σ ⊥ ). For H = 2D, the cylinders of the lower (upper) layer interact as hard needles with the cylinders of the upper (lower) layer, while in the corresponding case, H = 2σ ⊥ , ellipsoids interact as point particles with those in the other layer. We find it convenient to use the 2D number density as a control parameter, which is defined as ρ 2D = N/A, where N is the total number of particles and A is the surface area of the walls. We consider the following two possible scenarios: (1) the case of strong adsorbing walls (the adsorption energy is infinite) and (2) the case of non-adsorbing ones (the adsorption energy is zero). As we find qualitatively similar behaviors in these limiting cases we believe that cases with finite adsorption energies at any temperature should also behave similarly. In the first case, the particles are attached to the surfaces of the walls in such a way that 50% percent of them are at each side. In this case, particles can move and rotate freely in the x-y plane, while the out-of-plane orientational freedom is frozen. In the second case, we keep the assumption that the particles prefer the planar arrangement for their orientations, but the particles are allowed to move freely between the two confining walls. This case allows us to examine the effect of out-of-plane positional freedom, while the effect of out-of-plane orientational freedom is neglected. This approximation is based on the above-mentioned suppression of these degrees of freedom as a result of the wetting effect due to the interaction between the particles in a dense system.
From the possible density functional theories, we have chosen the second virial one, which is simple and capable to correctly capture the shape and pore width dependence of the IN transition properties of confined nonspherical particles [57]. Note that advanced density functionals such as the fundamental measure ones can provide more accurate transition properties. Nonetheless, they are also mean-field and exact only at low densities as the one we are implementing [58][59][60]. In the second virial formalism, the grand potential (Ω) is a functional of the local density (ρ) as follows (1)) , (1) where β = 1/k B T is the inverse temperature, 1 = (x 1 , y 1 , z 1 , ϕ 1 , θ 1 ) is an abbreviation for the position and orientation of a particle (and analogously for 2), where x, y, and z are the Cartesian coordinates, ϕ is the azimuthal angle, and θ is the polar angle. Furthermore, ρ(1) is the local density, f M (1, 2) is −1 (0) for two overlapping (non-overlapping) particles having locations 1 and 2, respectively, μ is the chemical potential and V ext is the external potential due to confinement. To determine the equilibrium grand potential and the local density, the functional derivative of equation (1) with respect to ρ(1) must be zero, i.e. δβΩ/δρ(1) = 0. We can show that the equilibrium density profile satisfies the following equation, The solution of equation (2) must be substituted into equation (1) to get the equilibrium grand potential of the confined system. The strong adsorption condition through the external potential prescribes that z = ±(H − σ)/2 and θ = π/2. The first condition forces the particles to be only on the lower and upper surfaces of the slit-like pore, while the second one ensures that the orientations of the particles are in the x-y plane. Note that σ corresponds to D and σ ⊥ for cylinders and ellipsoids, respectively. We assume that the pore is infinite in the x-y plane and that there is no in-plane positional order, i.e. the local density is independent of the x and y coordinates as we are interested only in the stability of 2D nematic order. Using these conditions we get from equation (2) that where i = 1 corresponds to z = −(H − σ)/2 and i = 2 to z = (H − σ)/2. It is also used in the above derivation that the excluded area between two particles being in layers i and j can be obtained from A ij ex (ϕ 1 − ϕ 2 ) = − dx 12 dy 12 f M (1, 2), where x 12 and y 12 are the coordinates of the center-to-center vector between the two particles (see the appendix for details). Equation (3) can be simplified by using the assumption that the density distribution of the lower and the upper layers are the same for symmetry reasons, i.e. ρ 1 (ϕ) = ρ 2 (ϕ) = ρ(ϕ)/2. It should be stressed that the calculations for the non-adsorbing cases do not use this symmetry assumption and never produce a thermodynamically stable symmetry-breaking solution. Thus, we believe that this assumption does not restrict generality. It is also worth using the normalization condition of the local density, d1ρ(1) = N. In the case of strong confinement this condition simplifies as follows: It is also useful to introduce the orientational distribution function, which is normalized as follows: f(ϕ) := ρ(ϕ)/ρ 2D and 2π 0 dϕ f (ϕ) = 1. After a straightforward calculation, equation (3) can be rewritten in terms of f(ϕ) in the following form The solution of this self-consistent equation provides the equilibrium orientational distribution function. We solve numerically the integrals of this equation using the trapezoidal quadrature and iteration to get f(ϕ). To measure the extent of the nematic order it is useful to introduce the nematic order parameter as follows: S = 2π 0 dϕ f (ϕ) cos (2ϕ), which is zero (non-zero) in the isotropic (nematic) phase. To locate the density of the IN phase transition, we perform a bifurcation analysis with cos(2ϕ) function in equation (4), which is described elsewhere [61]. For this purpose, we employ the Fourier expansion of the excluded areas, which can be written as where the Fourier coefficients are given by Here, δ is the Kronecker delta function. Equation (4) can be rewritten as a self-consistent equation for the order parameter by integrating both sides with cos(2ϕ). We get where S k = 2π 0 dϕ f (ϕ) cos (2 k ϕ) is a generalised order parameter. Note that S 0 = 1 and S 1 corresponds to the standard nematic order parameter (S). In the vicinity of the IN bifurcation point, the order parameters S k with k 1 are close to zero, therefore the approximation exp(−x) ≈ 1 − x can be applied to equation (6). Using this first-order Taylor expansion result, we get S = ρ 2D 4 A 11,2 ex + A 12,2 ex S, which is valid for S → 0. After rearranging this equation, we end up with the following expression for the bifurcation density of the IN transition: Note that ρ IN corresponds to the IN transition density in the case of second-order phase transitions, while it provides only an estimation for first-order phase transitions. We present the results of equations (4) and (7) in the Result section for ellipsoids and cylinders, respectively. For the non-adsorbing case, the out-of-plane positional freedom is turned on. Here, the summation for the two positions is replaced by a definite integral for the z variable in equation (3), i.e.
where the local density (ρ(z 1 , ϕ 1 )) is now position and orientational dependent and A ex |z 1 − z 2 | , ϕ 1 − ϕ 2 is the excluded area between two particles in that case, when they are |z 1 − z 2 | apart along the z-axis and they have orientations ϕ 1 and ϕ 2 . In order to use ρ 2D = N/A as an input of equation (8), we resort to the normalization condition of the local density, After integrating out both sides of equation (8) with z 1 and ϕ 1 , we get an expression for e βμ , which helps us to express the local density as a function of ρ 2D as follows We solve equation (9) numerically using the trapezoidal quadrature for both variables (z and ϕ) and iteration to determine the equilibrium density distribution at a given ρ 2D . To solve equation (9) we also need the expression of the excluded area, which is presented in the appendix for ellipsoids and cylinders. In addition, we calculate the global nematic order parameter (S), the local density (ρ(z)), and the nematic order parameter (S(z)) profile from the equilibrium ρ(z, ϕ). These quantities are given by and To locate the isotropic-nematic bifurcation density (ρ IN ) we start the iterative solution of equation (9) with a nematic trial function for ρ(z, ϕ) at such a high ρ 2D to guarantee the phase to be nematic (S > 0). Then, we gradually decrease ρ 2D and solve iteratively equation (9), while searching for the first ρ 2D value Figure 2. Isotropic-nematic phase transition of prolate (κ > 1) and oblate (κ < 1) particles in the ρ * -H * plane for strong planar adsorption. Phases are nematic above the curves and isotropic below them. The results are obtained from equation (7) for cylindrical and ellipsoidal shapes. The curves for oblate and prolate ellipsoids coincide with each other.
having S = 0. The collection of these densities as a function of H are ρ IN (H) for our confined system. We present our results for ρ IN , ρ(z), S and S(z) in the next section, where the unit of the length is σ, i.e. H * = H/σ and z * = z/σ. We make the 2D number density dimensionless by ρ * = ρ 2D min(L 2 , D 2 ) for cylinders and ρ * = ρ 2D min(σ 2 , σ 2 ⊥ ) for ellipsoids. As we consider oblate and prolate shapes, this definition is useful since it forces both shapes to have the same ρ * IN in the limit H * = 1.

Strongly adsorbing case
We start this section by showing the phase behaviour found for strong adsorption on the walls, where half of the particles are adsorbed tightly to one surface and the other half to the other. Thus, in this case, the particles form two fluid layers without out-of-plane orientational freedom, and no particles are placed in the interstitial region. The distance between the two flat surfaces is limited to maintain the q2D feature of the pore, i.e. the layers interact with each other. Therefore, we can see the two-layer fluid as a 2D non-additive binary mixture of anisotropic particles, where the in-layer interactions would be like ones, whereas the inter-layer interactions would be unlike ones. In other words, in equation (4) the contribution of like interactions would be given by A 11 ex , while that of unlike ones would correspond to A 12 ex . This analogy highlights the fact that we can envisage a 2D-like behaviour even when the system has two well-defined layers. Figure 2 presents our results for prolate and oblate shapes, where the H dependence of ρ IN is shown for hard ellipsoids and hard cylinders. In the H * = 1 limit, we get the strictly 2D system of hard ellipses and that of rectangles, since particles can interact only through their central cross-sections. In the H * = 2 limit, the two-layer fluid of ellipsoids corresponds to a 2D equimolar mixture of ellipses of equal size and shape, where the unlike interactions take place between point particles. In addition, the two-layer fluid of cylinders could be seen as an equimolar mixture of rectangles of equal size and shape with needle-needle unlike interactions. We observe that the IN phase transition is always second-order on the level of the theory, which should be a KT-type continuous phase transition according to simulation studies [11,25,40]. As we mentioned earlier, this feature of 2D phase transitions cannot be captured by the present mean-field theory [15]. Starting with the ellipsoids, we can see that the widening pore destabilizes the nematic order by shifting the IN transition in the direction of larger densities (see the dashed curves of figure 2). This is the expected behaviour, since it indicates that more ellipsoids are required to induce the nematic order as the pore becomes more spacious with increasing H. In the limit of H * = 2, we observe that one must double the number of particles (and ρ 2D ) to get the nematic phase since the fluid layers of the two walls completely decouples from each other. This additive behaviour can be easily understood in terms of the 2D mixture analog, with the emerging of point-point unlike interactions explaining their decoupling. As the number of the ellipsoidal particles is N/2 at both walls, the number of ellipsoids must be duplicated to regain a stable nematic phase. We can also see in figure 2 that the IN curves of prolate and oblate ellipsoids are identical provided that a very strong external edge-on aligning field is applied for oblate ellipsoids to maintain the same out-of-plane orientation. This perfect agreement is due to a couple of facts. On the one hand, there is an oblate-prolate symmetry for the excluded area (or volume), which results in the same IN transition packing fractions for oblate and prolate ellipsoids even for the three-dimensional bulk case [62][63][64][65]. On the other hand, we have always used the particles' shortest length as the unit of distance. It is also visible that the stronger shape anisotropy, which corresponds to increasing κ for prolate and decreasing κ for oblate shapes, stabilizes the nematic order. In the κ → ∞ (thin prolate) and κ → 0 (needle oblate) limits, ρ IN vanishes in accordance with MC simulation results of 2D hard ellipses [11].
The orientational order of hard cylinders is strikingly different from that of ellipsoids as ρ IN decreases monotonically with the widening of the pore for prolate shapes (κ > 1). This is also true for oblate cylinders (κ < 1), when H is below a certain κ dependent value of H. However, for oblates the nematic phase is always destabilized in the vicinity of H * = 2. The decreasing ρ IN means that fewer particles are needed in the pore for the stabilization of the nematic, despite that more room is available with increasing H. In the case of κ = 5, figure 2 shows that the required amount of particles to form the nematic order is about 20% less for the widest pore (H * = 2) with respect to the narrowest one (H * = 1). When comparing this percentage with the case of κ = 10 (lower panel of figure 2), it turns out that it goes down. Indeed, the decrease in the number of cylinders falls to 10% for κ = 10. This suggests that the stabilization effect of increasing H on the nematic phase weakens with increasing κ. The different behaviour of ellipsoids and cylinders can be clearly seen for the nematic order parameter curves as a function of H at a fixed surface coverage (ρ 2D ) (see figure 3). Recall that the two interacting layers move away from each other with increasing H. In the case of ellipsoids, increasing H results in a lowering nematic order parameter, while the opposite is found in the case of prolate cylinders, which is a strikingly different behaviour. As can be seen in figure 3, while S is a monotonic increasing function which is almost doubled for the cylinders, it is a monotonic decreasing function and reaches zero for ellipsoids as H * goes from 1 to 2. The oblate cylinders depict a curious behaviour as their S follows the trend of prolate cylinders for low H * (H * 1.5) and that of ellipsoids for high values of H (1.5 H * 2). As mentioned, the monotonic decreasing trend of S with H for hard ellipsoids at a fixed number of particles is the expected behaviour, while the increasing S(H) trend of cylinders seems illogical at a first glance.
A possible explanation for the unusual ordering behaviour with increasing H of hard cylinders can be given by using the following simple arguments. It is a well-known fact that the increasing shape anisotropy enhances the nematic order. In the seminal paper of Lars Onsager, it is shown that the packing fraction of the IN transition of 3D objects in bulk vanishes for κ → ∞ [66]. The same tendency is confirmed for the nematisation of hard ellipses and hard discorectangles (stadiums) using MC simulations [11,25]. In the case of hard needles of length L, which can be seen as the κ → ∞ limit of hard rectangles, it is observed that ρ IN L 2 ≈ 7 [10]. This suggests that the IN packing fraction of rectangles behaves as η IN = ρ IN a ≈ 7/κ, where a = LD is the area of the rectangle. From this equation, we can see that ρ IN is proportional to κ −1 a −1 , i.e. increasing the shape anisotropy (κ) and particle's area (a) lowers ρ IN . In the slit-like pore, there are competing in-layer and inter-layer interaction areas and shape anisotropies, which results in a very sensitive ordering behaviour. To explain the observed ρ IN trend with H, we make use of an empirical extension of ρ IN ∼ κ −1 a −1 , given that our system is analogous to an equimolar non-additive binary mixture, where the like interactions correspond to that of in-layer particles and the unlike interactions to that of inter-layer ones. For this purpose, we simply replace κ and a with their corresponding effective values, which are defined as κ eff := (κ in−layer + κ inter−layer )/2 and a eff := (a in−layer + a inter−layer )/2. Here, we consider the in-layer and inter-layer cross-sections of the bodies in parallel orientations. Using the average values of the interacting cross-sections we employ ρ IN,eff ∼ κ −1 eff a −1 eff to explain the shape and pore width dependence of ρ IN .
One can easily show that ρ IN,eff does not depend on H for a rectangular prism with lateral faces parallel and perpendicular to the walls. In this case, the in-layer and the inter-layer cross-sections are always the same rectangles. Therefore, the adsorbed prisms in the pore behave exactly as the 2D fluid of hard rectangles, and figures 2 and 3 would depict horizontal lines. The case of parallel prisms is in high contrast with the cases of ellipsoids and cylinders. Figure 4 depicts the H dependence of ρ IN,eff for ellipsoids and cylinders (parallel prisms would also show a horizontal line here). We can see that ρ IN,eff increases for ellipsoids because κ eff is independent of H as κ in−layer = κ inter−layer = σ /σ ⊥ , while a eff decreases with H as a in−layer = (π/4)σ σ ⊥ and a inter−layer = (π/4)σ σ ⊥ [1 − (H − σ ⊥ ) 2 /σ 2 ⊥ ]. This shows that the ellipsoids participate with decreasing interaction area (a eff ) in the formation of the nematic phase, as a consequence of weakening layer-layer coupling (decreasing a inter−layer ).
The situation is quite different in the case of prolate cylinders, because the inter-layer aspect ratio diverges at H * = 2, while the inter-layer interaction area goes to zero as follows: As the in-layer aspect ratio and interaction area are given by κ in−layer = κ = L/D and a in−layer = LD, respectively, ρ IN,eff decreases with increasing H * . Therefore, the stronger anisotropic interactions (increasing κ inter−layer ) between the neighboring layers are responsible for the enhancement of the nematic order. In the H * = 2 limit, the needle-needle interactions take place between the two adsorbed fluid layers, which renders the IN transition density to be zero. Note that triangular prisms also show the same effect when they are adsorbed to the wall with their sides. In this case, the effective width of the particles linearly goes to zero (the needle limit), and therefore, the enhancement of the nematic order should be even more pronounced than for cylinders.
Conversely to the prolate cylinders, the IN density behaviour of oblate cylinders cannot be explained with the help of the above arguments. Namely, the fact that the nematic order is strengthened with the widening of the pore up to a certain value of H (H * ≈ 1.4-1.5), and then is replaced by a strong destabilization effect up to H * = 2 (see κ = 1/5 and 1/10 cases in figures 2 and 3) cannot be explained by the simplified argument based on κ eff and a eff . The effective density is now given by ρ IN,eff ≈ κ eff /a eff , Figure 5. Isotropic-nematic phase transition of prolate (κ > 1) and oblate (κ < 1) particles in the ρ * -H * plane when the adsorption energy is zero. The phase is nematic above the curves and isotropic below them. Results are obtained from equation (9) for cylindrical and ellipsoidal particles. Oblate and prolate ellipsoids yield the same curves. Figure 6. The global nematic order parameter as a function of the pore width for confined ellipsoids and cylinders when the adsorption energy is zero. The 2D density is fixed for all shape anisotropies as follows: ρ * = 0.3 for κ = 5 and κ = 1/5, while ρ * = 0.06 for κ = 10 and κ = 1/10. The curves are the results of equation (12). Oblate and prolate ellipsoids yield the same curves.
because κ < 1 for oblate shapes. The cross-section areas of the particles decrease monotonically with increasing H since the slices of the cylinders move away from the largest central rectangular slice. As a eff decreases with H, the effective aspect ratio must compensate the destabilization effect of a eff to stabilize the nematic order. This can be achieved with decreasing κ eff , which would correspond to a more anisotropic rectangular shape. However, we get that κ eff increases by widening the pore up to H = D + √ D 2 − L 2 , where the shape of the cylinder's slice changes from oblate to prolate (κ inter−layer = 1). This shows that the stabilization of the nematic phase with increasing H for oblate cylinders cannot be explained in terms of the shape anisotropy and the cross-section area of effective 2D objects, as we always get that ρ IN,eff increases with H (see figure 4). The fact that the nematic phase is less stable with increasing H for D + √ D 2 − L 2 < H < 2D does also contradict the above arguments because the inter-layer cross-sections of the oblate cylinders become more needle-like as H → 2D, which should stabilize the nematic phase. Here, our simple explanation fails since it does not take into account that the contact points of two cylinders attached to different walls do not lie on a plane at a given H. Therefore, it is not possible to approximate the bilayer of oblate cylinders with a non-additive equimolar mixture of hard rectangles. It is curious how the 3D character of the group of points arising from the contact of two cylinders attached to different walls plays a key role for oblate shapes but not for prolate ones.

Non-adsorbing case
For obvious reasons, all the results for the strongly adsorbing and non-adsorbing cases coincide when H * = 1 (compare figures 2 and 5, and 3 and 6). Furthermore, we observe similar trends in the orientational ordering properties in the two cases. The nematic phase is destabilized for ellipsoidal shapes, while it is stabilized for the cylindrical ones by increasing the pore width. Figure 5 shows that the out-of-plane positional freedom weakens the destabilization (stabilization) effect of the walls for ellipsoids (cylinders). For oblate cylinders, the drop of ρ IN goes below 10% with κ = 5 when doubling H * from 1 to 2, which should be compared with the drop of more than 20% occurring for the strong adsorption. In the case of ellipsoids with κ = 5, the increment of ρ IN is about 23% at H * = 2, which was 100% in the strongly adsorbing case. These weakening stabilization and destabilization tendencies are due to the fact that the particles are distributed continuously between the two walls, and thus the average out-of-plane distance between the particles decreases, resulting in a continuous distribution of interacting areas and shape anisotropies. For ellipsoids, the average interacting area increases, while the average aspect ratio of the interacting areas keeps constant. Therefore, the increasing interaction areas causes the weakening destabilization of the nematic order. In the case of cylinders, the aspect ratio is responsible for the weakening stabilization of the nematic phase because the cross-sections of the interaction areas become less anisotropic. The same weakening tendencies happen in the case of oblate cylinders as the average out-of-plane distances between the particles decreases and the interactions between cylinders are closer to the interactions of 2D system of hard rectangles. It is also clear that the out-of-plane effects get even weaker with increasing the shape anisotropy of the particles because the effective shape anisotropy can change less with the spreading of the particles between the walls.
These findings manifest in the behaviour of the nematic order parameter (equation (12)), which is shown in figure 6. We can see that the change of S is smoother here than in the case of strong adsorption. The nematic phase of ellipsoids survives in the whole range of H, the growth of the nematic order weakens with H for prolate cylinders, and the destabilization of the nematic order occurs at higher H for oblate cylinders. The accompanying local density, equation (10) and local order parameter, equation (11) show some interesting features for prolate particles. Figure 7 shows that both ρ(z) curves (ellipsoids and cylinders) depict U-shaped profiles peaking at the walls and having a minimum at the pore center. These density profiles indicate that the particles accumulate at the walls despite the positional freedom in the pore. This entropic adsorption is more intense for ellipsoids than for cylinders. However, the orientational order of ellipsoids and cylinders exhibit qualitative differences. Cylinders are more ordered at the walls than in the middle of the pore, contrasting with the behaviour of ellipsoids. This different behaviour can be understood in terms of the interacting areas and anisotropies. The ellipsoids at the middle of the pore must interact with those at both sides having the same aspect ratio and relatively large interaction area. This contrasts with the absorbed ellipsoids since they have small interaction area with the ellipsoids adsorbed at the other wall. Therefore, ellipsoids in the middle of the pore feel a denser environment than those adsorbed, explaining the S peak at the middle. The case of cylinders is different because the interaction area decreases while the interaction anisotropy increases with distance. For example, for H * = 2, two cylinders staying at the opposite walls interact as two hard needles. Therefore, the stronger orientational correlation is responsible for the higher nematic order at the walls. Although the cylinders in the middle plane of the pore feel a denser environment, the interaction areas are less anisotropic and the extent of the nematic order is smaller.
Our results show that the phase behaviour of strictly 2D anisotropic systems cannot be reproduced accurately in real experiments when the magnitude of the out-of-plane positional fluctuations are at least in the order of the short length of the particles and the particles are weakly or moderately anisotropic. Note that out-of-plane orientational and positional fluctuations are always present in monolayers for different experimental confinements, such as slit-like pores, adsorbing solid surfaces, and liquid interfaces. Therefore, the experimental results for the orientational order and for the order of IN phase transitions of confined mesogenic particles should be considered with some reservations as strict 2D results. It may happen that the first-order nature of IN transitions observed for some experimental 2D setups [28,36] is partially due to the presence of out-of-plane positional fluctuations. To have negligible out-of-plane fluctuation effects on the nematic order, the particles must be very elongated or very thin in q2D experiments. Note that as we go beyond the H = 2σ limit, the coupling between the fluid layers of the walls is lost and the number of particles must be increased to induce nematic order even for the cylindrical shape. This means that the nematic phase stabilization effect of the out-of-plane positional fluctuations is present only for H 2σ, and for particles having a cross-section anisotropy that increases as moving out from their geometrical center.

Conclusions
We have shown that the general belief that the 2D nematic order is always destabilized by adding out-of-plane positional degrees of freedom to the anisotropic particles is simply not true for both rod-like and plate-like cylinders (see figure 8). By placing rod-like hard cylinders into a slit-like pore, we observed enhanced nematic order by increasing the wall-to-wall distance for a fixed 2D number density. The strongest ordering effect with increasing H is found for the strong adsorption case, where the particles are attached to the surface of the walls. Surprisingly, the nematic order parameter can be doubled by doubling H at a given 2D density. If the out-of-plane positional freedom is turned on while keeping the planar alignment, the above-mentioned effect persists but weakens. We expect that the incorporation of out-of-plane orientational freedom would have additional weakening effects on the enhancement of the nematic order of hard cylinders. Nonetheless, we do not expect a qualitative change of the trend. The reason is twofold. On the one hand, as simulation results show [57], the packing entropy of the cylindrical particles lead to a density profile that has a maximum very near the walls, and the density is smaller in the middle of the pore when H is smaller than twice the width of the particles. Consequently, the walls reduce the out-of-plane orientational fluctuations of the vast majority of the particles, that is, one rotational degree of freedom of the particles is effectively suppressed. On the other hand, the contribution of accessible polar angles is reduced by increasing the shape anisotropy of rod-like particles, and the elongated particle shape is necessary to obtain the nematic order. Moreover, we highlight here that other solid shapes such as triangular prisms may exhibit even more intensively this phenomenon.
The outcomes of the present density functional theory were compared elsewhere with Monte Carlo simulation results for confined freely rotating hard spherocylinders [57]. They agree qualitatively and predict that the stability of the nematic order is practically independent of H in the presence of both out-of-plane positional and orientational degrees of freedom in the q2D regime. Therefore, hard spherocylinders and square prisms can be considered as intermediate cases between the nematic stabilizing and destabilizing shapes. We found that ellipsoids are a typical destabilizing shape for the nematic order with increasing H. In brief, anisotropic 3D objects can be classified as stabilizers, neutral and destabilizers of the nematic order with increasing the pore width.
Our geometrical rationalization of the phenomenon by using the interacting slices of the particles works well for prolate ellipsoid and cylindrical shapes. However, it fails for prolate spherocylinders because it predicts the stabilization of the nematic order with H as the interaction slices of the particles become more anisotropic with increasing the inter-layer distance. Therefore, our empirical formula ρ IN,eff ∼ κ −1 eff a −1 eff cannot be applied safely even for prolate shapes. Oblate shapes are more complicated, as both stabilization and destabilization trends are observed in very narrow pores for a given shape. In this case, the majority of the contact points tend to lie farther away from a z = const. plane, and the effective shape of the cylinder cannot be well approximated by a simple rectangle. Instead, the sides of the rectangles depend on the relative positions of the cylinders and this effect turns more important as κ decreases (see the appendix).
To detect experimentally the out-of-plane positional fluctuation induced enhanced nematic ordering, we believe that the granular system of vertically vibrated granular cylinders can be an ideal playground. Along this line, the vertically vibrated granular cylinders in circular cavity exhibits orientationally ordered phases with defect structures [67]. To avoid the destabilization effect of the boundary, we think that the cylinders should be placed into elongated rectangular cavity, where, at least near the middle of the longer side of the cavity, the symmetry of the nematic phase matches with that of the boundary and the wall induced defects has a negligible effect.
We believe that our functional theory predictions for oblate shapes can be checked experimentally by inserting plate-like particles into a narrow slit-like pore, while applying an external field to force the plates to have an edge-on orientation. In this direction, a q2D rectangular system has been already set up experimentally by ordering plate-like cylinders into the edge-on direction by applying an external electric field [56], where the particles settle in the bottom of the cell and form isotropic and nematic phases.
An interesting collateral result of our work is the excluded area between two identical hard ellipsoids or cylinders. This quantity is required for computing the excluded volume (and second virial coefficient) through definite integration in the z coordinate. In turn, the excluded volume is a key parameter for the location of the IN transition density of 3D systems [15]. It is also useful for testing the overlap routines of hard body systems, which must be implemented in Monte Carlo simulations [68]. As the excluded volume for both ellipsoids and cylinders has been already reported using different geometrical methods [62,66], we could check the correctness of our method, which turned out to exactly reproduce the published results [62,66]. The overlap check is time-consuming for both hard cylinders and ellipsoids, and so new overlap detection algorithms are constantly being developed [69][70][71][72][73][74]. To ensure the reliability of these overlap routines, it is useful to calculate the excluded volume with the combined use of overlap routines and MC integration methods as it was done for hard cylinders elsewhere [68]. Our excluded areas make it possible to perform deeper tests for the overlap detection algorithms.
Finally, we would like to provide some outlook into the applications of the orientationally ordered monolayers. To construct complex nanostructures, it is very important to increase the surface density, the size of the nematic mono-domain, and the orientational order on a substrate [20]. For the case of single-walled carbon nanotubes, one possible solution is to shorten the nanorods, because the carbon nanotubes are very long and wavy [75]. Another option is the patterning of the substrate surface with stripes, which works very efficiently for increasing the density and the nematic order of gold nanowires [76]. In this regard, our work suggests another possible way to enhance the nematic order, i.e. by allowing the particles to slightly move apart from the plane of the substrate, or to slightly increase the width of the confining interface. To observe this enhancement of the nematic behaviour the shape anisotropy of the cross-section area of the particles must increase as moving away from its geometrical center, as it is the case for cylinders.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Excluded area of hard cylinders and ellipsoids
Because of the central role of the excluded volume of hard bodies in the second virial density functional theory there are many result in the literature related to the excluded volume in 3D [77] and the excluded area in 2D [78]. Even in his seminal paper [66], Onsager presented an exact computation for the mutual excluded volume of two hard cylinders. However, in our case, the system is not invariant under the translation along the z-axis, and therefore, we need the excluded area of the objects to apply the Onsager's second virial theory. This fact is explained in detail in the description of the theory given in section 2. The excluded area, A ex , is the area of the xy-plane figure enclosed by the curve which is drawn by the center of a particle while it moves around another fixed one in contact with it, for a parallel configuration, at a given z, and with a fixed orientation (for cylinders see the blue curve of figure A1). The contact vector σ joins the center of both particles in contact and points from the fixed one towards the other. In addition, φ is the azimuthal angle of σ. Thus, σ(φ) is a natural parametrization of the curve which encloses the excluded area. We find that certainly depends on the difference between the z coordinates of the particles, z 12 = z 2 − z 1 , and on the relative orientation of the particles, ϕ 12 = ϕ 2 − ϕ 1 . As it can be seen, only σ xy , the x-y component of σ, contributes to the integral, i.e. σ can be replaced by σ xy in the above formula. In the following we assume that z 2 z 1 , in the opposite case one can change the label of the particles.

A.1. Cylinders
In practice, it is difficult to handle σ xy as a function of φ. While particle 2 moves around 1 (the fixed one), and φ runs from 0 to 2π, the contact point describes a 3D trajectory (see the red curve in figure A1). We denote the z coordinate of this contact point by ζ, which is a function of φ. In our computation, we use the fact that σ xy can be more easily expressed in terms of ζ than of φ. Therefore, we use ζ for the parametrization of the σ xy curve instead of φ, and we change the integration variable of equation (A.1) accordingly. With this trick, we can avoid some mathematical complications and yield manageable expressions for A ex . The only problem is that the ζ(φ) function is not invertible everywhere, but this issue can be overcome.
In the following, we use the radius of the cylinders (R = D/2) besides its length (L). By cutting both cylinders by the z = ζ plain, we obtain two rectangles in contact (let us avoid the degenerated case in which one of them is a line segment). These are the yellow rectangles shown in figure A1. Furthermore, in figure A1(b) their projections into the z = z 12 plane are shown in blue color. It is easy to see that these rectangles have a side of length L another side with length d i (ζ) = 2 R 2 − (ζ − z i ) 2 , where i = 1, 2 is the particle index. Certainly, the second rectangle is rotated ϕ 12 with respect to the first one. Figure A1. A cylinder moves around a fixed one while they are in contact. The difference between their z coordinates, z 12 = z 2 − z 1 and between their orientations, ϕ 12 = ϕ 2 − ϕ 1 are constants (for simplicity, we have set z 1 = 0 in the figure). The yellow rectangles denote the section of the cylinders defined by the z = ζ plane containing the contact point. The projections along the z-axis of these rectangles into the z = z 12 plain gives the blue rectangles of the panel (b), their centers are connected by the x-y projection of the contact vector, σ xy , indicated by the green arrow. The blue curve denotes the border of the excluded area which is a plane figure. Despite this, the red curve, which is the path of the contact point drawn as the particle moves around the other one, is clearly a 3D curve. The figure related to the R < z 12 < 2R case, the 0 < z 12 < R case is slightly different (R is the radius of the cylinder). A multimedia version of this figure, which shows the whole path of the moving cylinder, is accessible online in the HTML document.
For cylinders, the φ ∈ [0, 2π] interval can be split conveniently into several parts. There are subintervals where ζ is a constant while φ (and thus σ) changes. These are the cases when the edge of the base (which is a circle) of the moving cylinder slips on the cylindrical surface of the fixed one. In these cases, the tangent line to the base is contained by the tangent plane of the cylindrical surface passing through the contact point. At a given z 12  Furthermore, there are also subintervals where ζ is a piecewise strictly monotonic (therefore invertible) function of φ (see the arc parts of the red curve of figure A1). In these cases, the vertices of the two rectangles touch each other. Therefore, the x-y projection of the contact vector, σ xy can be easily expressed in terms of ϕ 12 and the sides of the rectangles, L and d i (ζ). By changing the integration variable to ζ, the integral given by equation (A.1) can be performed. A lengthy but straightforward calculation yields the following final result, A ex (z 12 , ϕ 12 ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ L L sin ϕ 12 + 4 R 2 − ζ 2 0 + 4 R 2 − (z 12 − ζ 0 ) 2 cos ϕ 12 + 4 sin ϕ 12 if z 12 R L L sin ϕ 12 + 4 R 2 − ζ 2 0 + 4 R 2 − (z 12 − ζ 0 ) 2 cos ϕ 12 where 0 ζ 0 R is the solution of equation (A.2). Note that the total length of the subintervals where ζ(φ) = ζ 0 is a constant decreases with decreasing κ. Therefore, the contact points, on average, tend to lie farther away from the z = ζ 0 plane, i.e. the red curves shown in figure A1 get farther away from a plane figure. That is the reason why the approximation proposed in section 3.1 fails for the oblate case.
By integrating equation (A.3) respect to z from −2R to 2R one finally obtains the excluded volume. The result coincides exactly with the Onsager's result, which confirms the correctness of our procedure. Indeed, since our method is independent of his derivation, we can also claim that our outcomes are proof of the correctness of Onsager's result (his results are debated in [68] based on MC integration methods). We think our results are also useful for testing overlap detection algorithms.

A.2. Ellipsoids
Here, we present the results for the general case of triaxial ellipsoids with principal semi-axes a, b, and c. Note that in the main text we study only the spheroids, where a = σ /2 and b = c = σ ⊥ /2.
The method is similar to that described for cylinders, nevertheless, there are some differences. Ellipsoids are convex bodies with smooth surfaces, therefore at the contact point, they always present a well-defined common tangent plane. In other words, the normal vectors of the surfaces of the contacting ellipsoids are antiparallel with respect to the other. This leads to an equation useful for determining ζ, the z coordinate of the contact point via the quartic equation By cutting the two ellipsoids with the z = ζ plain we get two contacting ellipses of different sizes (for the degenerated cases they are points), where the second one is rotated by ϕ 12 respect to the reference one. The x-y projection of the contact vector of the ellipsoids, σ xy , points from the center of the reference ellipse to the center of the other. Unfortunately, the determination of this vector is not as easy as in the case of rectangles, but it can be straightforwardly solved [79,80]. In this case, there is no gain on changing the variable integral of equation (A.1). The final result can be written as follows, A ex (z 12 , ϕ 12 ) =