Extreme Poisson’s Ratios of Honeycomb, Re-Entrant, and Zig-Zag Crystals of Binary Hard Discs

: Two-dimensional (2D) crystalline structures based on a honeycomb geometry are analyzed by computer simulations using the Monte Carlo method in the isobaric-isothermal ensemble. The considered crystals are formed by hard discs (HD) of two different diameters which are very close to each other. In contrast to equidiameter HD, which crystallize into a homogeneous solid which is elastically isotropic due to its six-fold symmetry axis, the systems studied in this work contain artiﬁcial patterns and can be either isotropic or anisotropic. It turns out that the symmetry of the patterns obtained by the appropriate arrangement of two types of discs strongly inﬂuences their elastic properties. The Poisson’s ratio (PR) of each of the considered structures was studied in two aspects: (a) its dependence on the external isotropic pressure and (b) in the function of the direction angle, in which the deformation of the system takes place, since some of the structures are anisotropic. In order to accomplish the latter, the general analytic formula for the orientational dependence of PR in 2D systems was used. The PR analysis at extremely high pressures has shown that for the vast majority of the considered structures it is approximately direction independent (isotropic) and tends to the upper limit for isotropic 2D systems, which is equal to + 1. This is in contrast to systems of equidiameter discs for which it tends to 0.13, i.e., a value almost eight times smaller.


Introduction
Poisson's ratio (PR) is one of a few elastic moduli that can be used to describe deformation of elastic media. It is defined as a negative ratio of the transverse strain to the longitudinal strain, that occur as a result of an infinitely small change of the longitudinal stress [1]. In the case of isotropic three-dimensional systems, due to thermodynamic stability, PR can take values in the range ν ∈ [−1, 0.5]. This is different from the case of two-dimensional (2D) systems for which ν ∈ [−1, 1]. Both of the mentioned ranges result directly from the equation [2,3]: where: B-bulk modulus, µ-shear modulus and D-dimensionality of the system. Thermodynamic stability requires that both of the above moduli are positive (B > 0, µ > 0), leading to: One can easily see that substitution of D = 3 or 2 to the above expression results in the aforementioned ranges. It is worth noting, however, that the extreme positive value of PR in isotropic 2D systems (i.e., ν = +1) is not only a theoretical consequence of the above analytical formula, but it also has been presented on simple crystalline structures [4].
In this work we continue study described in the previous work [4], where we raised the issue of extreme values of PR of isotropic planar crystals (approaching the upper limit of Equation (2) for 2D systems, i.e., ν ≈ +1) consisting of binary hard discs (HD). Such systems are made of HD of two sizes which are very close to each other. Previously, we introduced a few simple isotropic structures modeled using mentioned particles. Now, we introduce two models based on the honeycomb lattice, analogous to a few of studied previously, but with much increased sizes of their unit cells, what reflects in much more computing power required for their examination. Furthermore, we modify these isotropic (honeycomb) models, gradually transforming them into atomic re-entrant structures [5][6][7], what in the case of structures made of continuous ribs, could lead to an auxetic effect (negative value of PR) [8][9][10][11][12][13][14][15][16][17][18]. The negative sign of PR in materials means that the elastic response of such structures to, e.g., tensile uniaxial stress, is an increase of their transverse size which is clearly counterintuitive behavior. It is followed by many interesting consequences [19][20][21][22][23], which are the reason why they can be applied in many fields, e.g., nanocomposite foams used in chemical sensors [24], electromagnetic shielding [25], damping applications [26] or auxetic foams, structures and textiles used in protective sport equipment [27]. Studies of simple models reveal properties exhibited by the structures investigated and can supply information on mechanisms leading to them. Transition between honeycomb and re-entrant structures and accompanying changes in the elastic properties of such systems have not been investigated before for the interatomic interactions considered in this work.
Due to the introduction of mentioned modifications, studied models are no longer isotropic, and their elastic properties can, in general, depend on the direction. The main goal of the present work is the analysis of changes in PR characteristics caused by the transition from regular honeycomb to re-entrant atomic structure, modeled using binary HD. It is worth noting the role of symmetry of the binary patterns in the studied models. The work includes both: (i) patterns with six-fold symmetry (honeycomb geometry) that maintain the isotropy of elastic properties [1], as well as (ii) patterns with two-fold symmetry (re-entrant and zig-zag geometries) which are anisotropic. Moreover, the considered systems can be classified as elastic metamaterials. First, even if the disc size difference is very small, all models containing isotropic patterns achieve very high PR, unattainable by equidiameter HD. Secondly, some systems containing anisotropic patterns behave like isotropic at very high pressures. Thus, like in typical metamaterials [6,28], it is obvious that the structure of the studied models is very important.
The structure of the article is as follows. In the next Section 2 the models studied and the methods of solving them are presented. In the Section 3 results obtained are discussed. The last Section 4 contains the summary and conclusions.

Models and Methods
In this section, divided into three subsections, we briefly describe Section 2.1 the models studied, Section 2.2 the simulation method used and Section 2.3 the simulation details.

The Model
All the structures discussed in this work consist of HD with two slightly different sizes (determined by their diameters σ i ). The hard interparticle interaction potential u ij for such binary mixtures can be described by [4]: is mean value of diameters of discs i and j. HD constitute the 2D counterpart of the three-dimensional (3D) hard spheres. Both the systems are fundamental models of condensed matter phases (crystals [29][30][31] and fluids [32]) and melting in 2D and 3D, respectively. They can also be used to model colloids [33] and granulates [34,35].
The simulations were conducted for two different honeycomb-based structures, S61 and S91 and their "inversions", S61i and S91i (see Figures 1 and 2). The honeycomb model was chosen due to its obvious isotropy (for small deformations) resulting from six-fold symmetry [1]. The second reason for choosing this geometry among the others presented in [4] is related to the possibility of a relatively easy and "smooth" transition to the reentrant geometry, as shown later in Figure 3. As for the "inverted" structures, the authors also wanted to verify how the mutual replacement of "larger" and "smaller" atoms forming the binary patterns affects the elastic properties of the considered structures.

S61
S61i Figure 1. Images of S61 (left) and S61i (right) binary disc systems with an isotropic pattern of six-fold symmetry. The value '61' comes from the number of "white" (on the left) or "black" (on the right) atoms inside a "frame" made of "black" or "white" atoms, respectively. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively. Figure 2. Images of S91 (left) and S91i (right) binary disc systems with an isotropic pattern of six-fold symmetry. The value '91' comes from the number of "white" (on the left) or "black" (on the right) atoms inside a "frame" made of "black" or "white" atoms, respectively. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively.

S91 S91i
Sizes of typical samples representing the studied systems can be seen in Figures 1 and 2, and Figures A1-A4. The numbers after 'S' in their names come from the quantity of discs with one diameter that are contained in a hexagonal "frame" formed by discs with the second diameter (e.g., open circles surrounded by a "frame" formed by black dots, on the left in Figure 1). The "inverted" structures (on the right in Figures 1 and 2) were created by replacing black dots by open circles and vice versa, and the letter 'i' was added to their names. As it was just mentioned, black dots and open circles refer to HD with slightly different size. The first ones (black dots) stand for "larger" atoms with diameters equal to (1 + δ)σ, while the second ones (open circles) refer to "smaller" atoms, with diameters equal to (1 − δ)σ, where σ is the unit of length of the simulated systems and δ = 0.0005. Of course, such a small diversity of sizes of binary particles is impossible to see with the naked eye in the presented figures.
The honeycomb-based S61 and S91 structures were only a prelude to our main analysis. Furthermore, we modified the atomic structure of these "base" configurations, gradually transforming them into re-entrant structures. The transition example is shown in Figure 3, starting from the structure S61. As one can see, a single side of a "black" hexagon in S61 structure consists of six atoms. By focusing only on the horizontal sides and moving to the next structures (the direction is indicated by arrows) it can be seen that, depending on the structure, the length of these sides changes: 8 (S61-z1), 10 (S61-z2) and so on. The numbers k = 1, 2, 3 after 'z' (S61-zk) and 'r' (S61-rk) indicate how many atoms have been added to each side of any horizontal "line" forming the frame of the base (S61) structure. As the length of the horizontal sides of the frame changes during the transition, the slanted sides are "pushed" inward, gradually modifying the geometry from honeycomb to re-entrant. Using the same procedure, it is possible to transform the S61i structure, and after a small modification (due to different sizes of the unit cells), also the two remaining structures, S91 and S91i. Due to the large number of additional structures obtained using the procedure described in the Figure 3, we do not present them here, but all of them are included in the Appendix A. . Example transition from honeycomb-based S61 structure (isotropic, with six-fold symmetry) to re-entrant (r) structures (anisotropic, with two-fold symmetry). Transition structures are called "zigzag" (z) structures, due to the shapes of their vertical sides.

The Method
The simulations were carried out using Monte Carlo (MC) method in the isobaricisothermal ensemble (N pT) with variable shape of the periodic box. The method was described in detail elsewhere [36][37][38] and here we discuss only a few of its key aspects from the point of view of the conducted research.
All the simulated particles were contained in a 2D parallelogram box, on the edges of which periodic boundary conditions were applied. The box of periodicity can be described by a square matrix, the columns of which are the vectors of its sides: The non-diagonal elements, which were kept equal to each other making the matrix in Equation (4) symmetric, a y = b x were very close to zero throughout the simulation, thus the shape of the periodic box was almost rectangular (with only a small deviation). Marking the box matrix of the reference (average) state as H, representing the mean matrix h during the MC simulation (i.e., H ij = h ij ), accumulated from the moment in which the system reaches the thermodynamic equilibrium state, the deformation tensor ε can be calculated from the expression [38]: where I is the identity matrix and the dots represent the matrix products. Then, the elastic compliance tensor S can be obtained by analyzing the fluctuations of the strain tensor elements: where: T-temperature, k-Boltzmann constant, V p = | det(H)|-average (2D) volume of the system at a pressure p, and . . . means averaging in the N pT ensemble.
Only the honeycomb-based structures (S61, S91 and their "inverse") introduced in the previous subsection are isotropic and their PR has the same value regardless of the direction. However, the modified structures ('z' and 'r') are anisotropic and their PR can, in general, depend on the direction. Knowing all the S ijkl elements, it is possible to calculate the value of PR in any direction on the plane in which the deformation occurs, determined by the unit vectorn. In the 2D spacen is fully defined by a single angle φ:n = (cos φ, sin φ). The reaction of the system takes place in the lateral directionm, so:m = (− sin φ, cos φ). After substituting these versors into the definition of PR, one can obtain [39,40]: where the Voigt notation was used (S ij = S ji ), i.e.,: Due to the periodicity of ν(φ) in Equation (7), with the period equal to π, in the further discussion we assume φ ∈ [0, π].
The elastic properties of the systems were analyzed for a range of reduced pressures p * , determined by the expression: p * = 10 1+i/3 , for integer i ∈ [1,12] or even i ∈ [1,15] in some cases. This means that the reduced pressure was always within the range: 21.544 < p * ≤ 10 6 . The smallest value of p * was much above the melting point of equidiameter HD system. Therefore, in each case the systems formed crystalline structures.
As already mentioned in Section 2.1, the parameter δ, responsible for the binary nature of the studied systems was equal to 0.0005 throughout this work. It was not our goal to investigate the effect of this parameter on systems properties, so it remains an open question and could be the subject of future research. In the present work, we focus on a single δ value and study the effect of the transition between honeycomb and re-entrant atomic structure of binary HD (see Figure 3).
For each structure, 10 independent simulations have been conducted, whose typical length was equal to 10 7 cycles (trial steps per particle), preceded by an equilibration stage of 10% of this length. Although such simulation lengths were sufficient in most cases, some of the structures at very high p * needed even longer runs (see Appendix B for more detailed information).

Results and Discussion
As mentioned at the end of previous Section, for some of studied structures at very high p * , the typical length of a simulation run (10 7 cycles) was not enough. At the beginning of this Section, it is worth to clarify this issue. In Figure 4, analysis of degree of equilibration of the system throughout the simulation, for structure S91 at p * = 10 5 is shown. The average components of the simulation box of the system, h ij , were calculated using small "datasets", consisting of 100k cycles only. The larger number of a "dataset" means that it was collected at a later stage of the simulation run. Performing such an analysis, when using the MC method, allows one to determine how long simulation runs are required for studied systems. It can be seen that the calculated quantities of S91 at p * = 10 5 , even after such a run of 3.3 × 10 7 cycles, still did not yet reached the equilibrium state. It should be noted, however, that in the typical run, equilibration took much less than 10 6 cycles, which were discarded prior to calculations of physical quantities.  pressure p * = 10 5 . The averages were calculated basing on short excerpts from the MC simulation (each consisting of 100k cycles) and plotted as a function of the number of such a single "dataset". Increasing number of a "dataset" represents the extending length of the simulation.
In Figures 5-9, the angular dependencies of PR are shown. The results of ν(φ) of all the structures studied, could roughly be grouped into 10 sets. Dependencies of structures in a given group showed slight differences, but were qualitatively similar, which was the case for their classification. In order to simplify the discussion, here we present the charts only for a single representative from each group (it is always the one first mentioned (bold), in the following paragraphs). The charts for all the studied structures are included in the Appendix B.
In Figure 5a results for group consisting of: S61, S61i, S91 and S91i are shown. As it was expected from the honeycomb-based structures, which were isotropic, the PRs were independent of φ, having the same value for any angle value. It should be noted that the isotropy of the studied systems was not artificially forced in any way, e.g., by enforcing expected relations between the coefficients in Equation (8), such as: S 11 = S 22 , S 13 = S 23 = 0 or S 33 = 2(S 11 − S 12 ). In each case, all the elements of the elastic compliance tensor S from Equation (6) were calculated independently, based on the simulation results. In Figure 5b results for group consisting of: S61-z1, S61-z2, S91-z1 and S91-z2 are shown. A slight anisotropy can be seen, apparently for the intermediate p * studied, i.e., for i ∈ {6, 7, 8}). At low (i < 6) and high (i > 8) pressures, PRs were almost isotropic, like in the first group. When discussing the later groups of structures, we will focus on intermediate-to-high pressures only, skipping the low ones (i 5), because, as one may notice, at low pressures all the structures behaved approximately the same.
In Figure 6a,b result for groups consisting of: (a) S61i-z1, (b) S61-z3 and S91i-z1 are shown. PRs of structures in these groups were visibly anisotropic in the intermediate p * region, having consecutive minima and maxima at φ = n π 2 and n π 2 + π 4 , respectively, with n being an integer. These two groups were very similar, but we decided to distinguish them from each other, due to the evident difference in peak-to-peak values. At high pressures, PRs began to manifest an isotropic nature, which was a common feature of both the considered groups. In Figure 7a,b results for groups consisting of: (a) S61i-z2, (b) S61i-r1, S61i-r2, S61i-r3, S91i-r1, S91i-r2 and S91i-r3 are shown. The dependencies in these groups were quite similar to the previous two, with the difference, however, that at high pressures they still remained anisotropic. As the pressure grew, the central minima became narrower and sharper, but their peaks seemed to tend towards certain values. The difference in these limit values was the reason why it was decided to distinguish the discussed groups, and they were approximately equal to: (a) 0.296 for S61i-z2 and (b) 0.538 for S61i-r1. It should be noted that the nature of the curves in Figures 6 and 7 was qualitatively different, as in the case of the former one, the central peak displacement with the increasing pressure was not inhibited before reaching the value of 1. In Figure 8a results for group consisting of: S61i-z3 and S91i-z2 are shown. The anisotropy of the systems at high pressures was evident and the minima/maxima distribution was the same as in the previous two groups. However, as it can be readily seen, their peak-to-peak values are much smaller and the maxima do not tend to 1.
In Figure 8b results for group consisting of: (S61-r3, S91-r3) are shown. As it can be seen from the Figure 3, structure S61-r3 (and S91-r3, with the similar transition procedure presented in the figure, but for S91-based systems) comes directly after S61-z3 (S91-z2), during the transition process. Both, S61-z3 (see Figure 6b) and S91-z2 (see Figure 5b) are almost isotropic at high pressures, wherein the former showed a more pronounced anisotropy at intermediate pressures. Transition from S91-z2 to S91-r3 led to a subtle change in the nature of the PR dependence at intermediate pressures (i ∈ {6, 7, 8}), where the peak-to-peak values increase a little, although at high pressures it is still isotropic (see Figure A29). In the analogous case of transition from the last "zigzag" (z3) structure to the first "re-entrant" (r3) with S61-based structures, the greater overall degree of anisotropy ended in a more complex characteristic (see Figure 8b). As can be seen, the maxima at φ = nπ (n being an integer) at pressures corresponding to (i ∈ {9, 10}) exceeded ν = 1, which was the upper limit for isotropic 2D systems. Presumably, at even higher pressures, the characteristics would tend to the isotropic one, corresponding to this limit, as in the S91-r3 case. The last two groups of the "re-entrant" structures are shown in Figure 9a (S61-r1, S91-r1) and Figure 9b (S61-r2, S91-r2). The PR dependencies of both of these groups contained minima and maxima at intermediate pressures, typical for the structures presented previously in Figures 6 and 7. Furthermore, the structures with narrower unit cells showed higher peak-to-peak values than those with the wider ones. As pressure increased to very high values (i > 8), the PR dependencies of r1 and r2 structures began to differ in an even more interesting way. The structures with narrower unit cells behaved in a similar way to those presented in Figure 7a,b, but their central minima tended to much higher PR values of approximately 0.95 (S61-r1) and 0.85 (S91-r1). While the narrower structures had a slight recess around φ = π 2 at high pressures, the wider ones, on the other hand, had a small hill. As it was mentioned in the beginning of this Section, the sufficient length of a simulation run grew immeasurably as pressure increased to very high values, thus, the question of whether these hills would flatten out at even higher pressures than applied in this work, remained open. Nevertheless, based on the results available, at least one interesting feature of the studied systems can be seen. Namely, they exceeded the value +1, which was the upper limit of PR for isotropic systems. It is possible, since the systems under consideration were anisotropic. The mean values of PR as functions of the inverse of the reduced pressure, (p * ) −1 , are plotted in Figure 10a for S61-based, and in Figure 10b for S91-based structures. As it can be seen, in both cases only some of the "inversed" structures (marked in red in the figures) did not tend to ν = 1. Two of them, S61i-z3 and S91i-z2, tended for a certain value close to 0.4, for which the "mustache"-shaped characteristics in Figure 8a were responsible. Another seven structures, S61i-z2 and all the "inversed re-entrant" (r1 to r3, based on both, S61i and S91i), achieved high values of PR, but they were hindered in reaching the exact value of ν = 1 by the central minima visible in Figure 7a,b. As was mentioned above, as the p * increased, the minima became narrower, but they never completely disappeared as in the case of S61i-z1, S61-z3 or S91i-z1 (see Figure 6a,b). S61i-z2 S61i-z3 S61i-r1 S61i-r2 S61i-r3  Figure 10. Mean values of PR ν as a function of (p * ) −1 for: (a) S61-based and (b) S91-based binary discs structures. The images of "zigzag" (z) and "re-entrant" (r) structures are included in the Appendix A. Black symbols and lines refer to the "regular" structures (i.e., in the left columns in Figures 1 and 2, Figures A1-A4), and the red ones refer to the "inverted" structures (in the right columns in the mentioned figures).

Summary and Conclusions
Elastic properties of two-dimensional (2D) crystalline structures formed by hard discs (HD) of two different diameters were studied by computer simulations using the Monte Carlo method in isobaric-isothermal ensemble. The angular dependencies of Poisson's ratio (PR) were presented together with their mean values, expressed as functions of the inverse of the reduced pressure. 26 of the studied structures, based on a honeycomb geometry (S61 and S91) were classified into 10 groups whose characteristics show qualitative differences (see . While these differences are evident at intermediate pressure values, most often they disappear at relatively low and at the high ones, as most of the structures strive for isotropy. Although none of the studied structures is either axetic or partially auxetic, the obtained results reveal an interesting and surprising behaviour which may be exploited when constructing new metamaterials at nano-level. The PRs of the vast majority of the studied structures at very high pressures tend to +1, which is the upper limit of 2D isotropic systems. Only two of the structures certainly do not reach this value. It is worth to add that, for the pressure range studied, seven other structures attain PRs around 0.995, due to the nature of their angular dependencies showing a non-vanishing deep at the external stress oriented vertically. The tips of these peaks tend to intermediate values that differ between the three different groups in which the phenomenon occurs. Nevertheless, a very interesting result of the work is that many structures with very different geometries behave almost identically under high pressure conditions. One should also stress that the behaviour of PR observed for the structures studied in this paper is qualitatively different from the dependence of PR observed in systems of identical discs. Namely, in the latter systems, PR decreases monotonically with increasing pressure reaching the value close to 0.13 [41] at the infinite pressure limit (close packing), what is almost eight times less than in the systems studied here.
The last but not less important result is related to the structures with the re-entrant geometry. Namely, the structures with wider unit cells (marked as r2 and r3), at some (high) pressures show the PR values exceeding its upper limit for 2D isotropic systems. This was possible due to the anisotropic nature of the systems under certain conditions. The authors speculate that this anisotropic nature may tend to an isotropic one (with PR equal to +1) at even higher pressures, but the computational cost of the simulations that could demonstrate this behavior is so high that this question will be answered in future.

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

Abbreviations
The following abbreviations are used in this manuscript: isobaric-isothermal ensemble r re-entrant z zigzag

Appendix A. Images of Modified S61 and S91 Structures, Studied in the Work
In this appendix the structures which were not defined in the main part of this paper are presented. The structure names are explained in figure captions.

S61-z3
S61i-z3 Figure A1. Images of S61-zk (left) and S61i-zk (right) "zigzag" binary disc systems with an anisotropic pattern of two-fold symmetry. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively. The numbers k = 1, 2, 3 after 'z' indicate how many atoms have been added to each side of any horizontal "line" forming the frame of the base (S61 or S61i) structure.

S61-r1
S61i-r1 S61i-r3 Figure A2. Images of S61-rk (left) and S61i-rk (right) "re-entrant" binary disc systems with an anisotropic pattern of two-fold symmetry. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively. The numbers k = 1, 2, 3 after 'r' indicate how many atoms have been added to each side of any horizontal "line" forming the frame of the base (S61 or S61i) structure.

S91-z2
S91i-z2 Figure A3. Images of S91-zk (left) and S91i-zk (right) "zigzag" binary disc systems with an anisotropic pattern of two-fold symmetry. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively. The numbers k = 1, 2 after 'z' indicate how many atoms have been added to each side of any horizontal "line" forming the frame of the base (S91 or S91i) structure.

S91-r1
S91i-r1 S91i-r3 Figure A4. Images of S91-rk (left) and S91i-rk (right) "re-entrant" binary disc systems with an anisotropic pattern of two-fold symmetry. Black dots and open circles represent the HD of various sizes: "large" (with a diameter equal to (1 + δ)σ) and "small" (with a diameter equal to (1 − δ)σ), respectively. The numbers k = 1, 2, 3 after 'r' indicate how many atoms have been added to each side of any horizontal "line" forming the frame of the base (S91 or S91i) structure.