Full phase diagram of active Brownian disks: from melting to motility-induced phase separation

We establish the complete phase diagram of self-propelled hard disks in two spatial dimensions from the analysis of the equation of state and the statistics of local order parameters. The equilibrium melting scenario is maintained at small activities, with coexistence between active liquid and hexatic order, followed by a proper hexatic phase and a further transition to an active solid. As activity increases, the emergence of hexatic and solid order is shifted towards higher densities. Above a critical activity and for a certain range of packing fractions, the system undergoes MIPS and demixes into low and high density phases; the latter can be either disordered (liquid) or ordered (hexatic or solid) depending on activity.

We establish the complete phase diagram of self-propelled hard disks in two spatial dimensions from the analysis of the equation of state and the statistics of local order parameters. The equilibrium melting scenario is maintained at small activities, with coexistence between active liquid and hexatic order, followed by a proper hexatic phase and a further transition to an active solid. As activity increases, the emergence of hexatic and solid order is shifted towards higher densities. Above a critical activity and for a certain range of packing fractions, the system undergoes MIPS and demixes into low and high density phases; the latter can be either disordered (liquid) or ordered (hexatic or solid) depending on activity.
Active materials are out-of-equilibrium systems in which the dynamics of their elements break detailed balance [1]. Examples can be found in living systems, e.g. the collective motion of large animal groups [2,3], bacteria swarming [4], and the formation of traveling fronts of actin filaments [5], as well as in synthetic ones, like self-propelled grains [6] or self-catalytic colloidal suspensions [7]. Despite such diversity, the emergence of activity-induced collective behavior is captured by minimal models that yield accurate descriptions and shed light on their universal character. A key example is the Active Brownian Particles (ABP) model which considers spherical self-propelled particles with only excluded volume interactions [5,8,9,[11][12][13]. A hallmark of active particle systems is that at high enough density and activity, self-propulsion triggers a motility-induced phase separation (MIPS) into a low-density gas in coexistence with a high-density drop [5,[11][12][13][14][15][16][17], resembling the equilibrium liquid-gas transition but in the absence of cohesive forces and without a thermodynamic support [18,19].
Although active particles can in principle move in 3D, in most experimental set-ups they are confined to 2D. Most studies of 2D ABP focused on MIPS, and have therefore been largely restricted to intermediate densities [5,[11][12][13][14][15][16][17][18][19]. In contrast, their solidification, or melting, has received little attention [20,21], and the connection between the high Pe behavior and the equilibrium physics as Pe → 0 has been, surprisingly, disregarded. In particular, the fate of 2D melting (with its intermediate hexatic phase) under active forces, has only been investigated for dumbbell systems [22], where MIPS is continuously connected to the passive liquid-hexatic coexistence. This result sheded new light on the very nature of MIPS and showed the importance of exploring the full phase diagram at high densities. In this Letter, we address this issue in the paradigmatic ABP model.
Melting in 2D is a fundamental problem that has remained elusive despite decades of intensive research [23,24]. The transition was initially claimed to be first order [25] and later argued to follow a different scenario, with an intermediate hexatic phase, separated by continuous transitions mediated by the unbinding of defects [26][27][28]. More recently, numerical simulations [1,3,29] followed by experiments on colloidal monolayers [32], clarified the picture. They indicate that melting of passive hard-disks takes place in two steps: as the packing fraction is increased, a first-order transition between the liquid and hexatic phases occurs, followed by a continuous Berezinskii- Kosterlitz-Thouless (BKT) transition between the hexatic and the solid. The hexatic phase exhibits quasi-long-range orientational order and short-range positional one, while the solid phase has quasi-long-range positional and long-range orientational order. Liquid and hexatic phases coexist close to the liquid phase, within a narrow interval of packing fractions.
Here we examine how activity affects the phase behavior of 2D systems of isotropic particles, from the dilute regime to close packing (φ cp ≈ 0.91). We establish the complete phase diagram of 2D ABP spanning a broad range of activities, see Fig. 1. We show that the twostep melting scenario at Pe = 0 is maintained at finite but small activity, with a coexistence region between active liquid and hexatic phases (black area). Above, an active hexatic phase exists for all the explored activities (blue sector). Strikingly, active disks arrange in a hexatic phase in a larger density range than passive ones. At higher densities, orientational long-range and positional quasi-long-range order emerge for any activity, signal- ing the presence of an active solid phase (orange region). The liquid-hexatic and hexatic-solid transitions shift towards higher densities with Pe, meaning that activity destabilizes the ordered phases. At high enough activity (Pe 35), we identify the boundaries of MIPS using both pressure measurements and density distributions (black and white symbols). The MIPS region broadens as activity increases and eventually crosses the hexatic and solid transition lines. Such results show that (i) MIPS prevails over the hexatic and solid phases and (ii) MIPS generates a phase separation between a dilute and a high-density phase, which can either be liquid, hexatic or solid, as activity is increased. We consider N overdamped ABP, in a square box with volume V = L 2 and periodic boundary conditions. They self-propel under a constant modulus force F act along n i = (cos θ i (t), sin θ i (t)) and obey with r i the position of the center of the ith particle, r ij = |r i − r j | the inter-particle distance, and a short-ranged repulsive potential, U (r) = 4ε[(σ/r) 64 − (σ/r) 32 ] + ε if r < σ d = 2 1/32 σ and 0 otherwise. The terms ξ and η are zero-mean Gaussian noises that verify The units of length, time and energy are given by σ d , τ = D −1 θ and ε, respectively. We fix D θ = 3γk B T /σ 2 d and vary the packing fraction φ = πσ 2 d N/(4V ) and Péclet number Pe = F act σ d /(k B T ) by tuning L and F act at fixed γ = 10 and k B T = 0.05. The integration of Eqs. (1) used the velocity Verlet algorithm implemented in LAMMPS [33,46] . Simulations ran with N = 256 2 particles, scanning the parameter space φ ∈ [0 : 0.9] and Pe ∈ [0 : 200]. With less (N = 128 2 ) and more (N = 512 2 ) particles we explored finite size effects.
The equation of state. Our first estimate of the phase boundaries is given by the φ dependence of the mechanical pressure [18,34] with ∆P = P − P G and P G = N k B T /V the ideal gas pressure. The first term, P act , quantifies the effect of F act , the so-called active or swim pressure [35,36]. The second one, P int , is the standard virial term due to particle interactions. The definition in Eq. (2) is a state function for isotropic ABP such that P (φ) defines an equation of state [36]. (This does not hold generically in active systems for which the pressure can, for instance, depend on the details of the interaction between the particles and the confining walls [37].) In the dilute limit we recover the ideal gas law P V = N k B T eff = N k B T (1 + Pe 2 /6), at an effective temperature that is compatible with the one that stems from the fluctuation-dissipation relation in the late diffusive regime [38][39][40]. The equation of state for zero and weak Pe is shown in Fig. 2 (a). P (φ) is roughly flat in a narrow φ interval for Pe 3. A zoom over this area in the Pe = 1 case evidences a double loop structure characteristic of phase coexistence, see Fig. 2 (b). Although the equalarea Maxwell construction, that allows to directly extract the binodals, cannot be readily applied for Pe > 0 [19,36], we use it by extension of the passive disks analysis [1], as a first identification of the coexistence region (black dots in Fig. 1). Beyond Pe = 3, we do not find evidence for coexistence until the high-Pe regime where MIPS is attained. For Pe 35 the P (φ) curves become flat in between two densities. Representative curves at 10 ≤ Pe ≤ 50 are displayed in Fig. 2 (c). As it has been recently reported [18,34], at very high Pe, the pressure drops abruptly at the vicinity of MIPS (see Fig. 2 (d)), as a consequence of the existence of a metastability region with a very large nucleation barrier [18]. We obtain the limits of MIPS with an extrapolation of the flat part of P (φ) across the pressure jump (or spinodal), as illustrated in Fig. 2 (d) for Pe = 100. Previous numerical studies used the local density probability distribution functions (PDF) to locate the MIPS region, see e.g. [11,22]. For the sake of completeness, we searched for the limits of a double peak structure of these PDFs, finding the open symbols in Fig. 1  . The vertical and horizontal boxes represent the location of the parameters in the phase diagram (see Fig. 1). An exponential fit to the liquid data and an algebraic decay with power η = 1/4 are shown in (a).
Orientational order and the hexatic phase. We put the orientational order to the test using the hexatic order parameter ψ 6 (r j ) = N −1 j Nj k=1 e iθ jk , where θ jk is the angle formed by the segment that connects the center of the jth disk and the one of its kth (out of N j ) nearest neighbor (found with a Voronoi tessellation algorithm) and the x axis. We studied its correlation function g 6 (r = |r j − r k |) = ψ 6 (r j )ψ 6 (r k ) / ψ 2 6 (r j ) and kurtosis or Binder parameter U 4 = 1 − ψ 4 6 (r j ) /(3 ψ 2 6 (r j ) 2 ), see Figs. 3 and 4, respectively.
We use the change of behavior of g 6 (r), from exponential (active liquid, in black) to algebraic r −η (active hexatic, in blue), as a criterion to locate the hexatic transition (blue symbols in Fig. 1). In the hexatic (blue) region the power law decay is maintained, with exponent η taking a value close to the BKT η = 1/4 at the transition but varying with φ and Pe. These data are compatible with the behavior of the Binder cumulant, U 4 , that in the scale of the main panel in Fig. 4 has a common intersection point, proving the transition. The zoom in the insert shows a weak remanent N -dependence that would be compatible with a first order phase transition [43,44]; however, the accuracy of our data is not enough to draw such a conclusion and, moreover, a second order transition is consistent with the absence of phase coexistence found above Pe ≈ 3. As illustrated in Fig. 3 (b), activity shifts the emergence of orientational quasi-long-range order to higher densities. Orientational order and coexistence. The maps of the local hexatic order parameter and the PDFs of its modulus, shown in Fig. 5, provide clues to understand the difference between the two sectors with phase separation at low and high Pe. Close to Pe = 0 (a) the PDF is bimodal, with two peaks of roughly the same height for this choice of parameters. The map in the insert proves the existence of a ramified but large (of the order of the system size) region with the same local hexatic order. Under the dynamics this region changes form but the portion of surface that it occupies remains stable. These results are in perfect correspondence with the data for the local densities (see the SM in [42]). In the MIPS region, instead, the map shows many different colors associated to diverse local orientational ordering that do not extend over a long distance, even at long times. Under the dynamics the color pattern changes considerably, with breaking and recombination of blocks. Differences in the maps are translated into differences in the PDFs. The secondary peak close to |ψ 6,j | = 0.9 in Fig. 5 (b) is due to the interfaces between areas with almost perfect orientational order. Additional maps in other sectors of the phase diagram, PDFs of |ψ 6,j |, correlation functions and global hexatic order parameter Ψ = N −1 | j ψ 6,j | measurements are given in the SM.
Positional order and the solid phase. Since it is hard to assert whether g 6 acquires long-range order or does not decay at the length-scales of our finite-size box, we looked for solid quasi-long range positional order, that should be evidenced by an algebraic decay of at the wave vector q 0 at the maximum of the first diffraction peak of the structure factor S(q) = N −1 i,j e iq·(ri−rj ) . The change in the C q0 decay, from exponential (hexatic) to algebraic (solid) for several Pe and φ, see e.g. Fig. 6, yields the orange points in the phase diagram above which lies the solid. Activity introduces non-equilibrium fluctuations that destabilize order and melt the solid. In (a) we show an exponential and an algebraic decay with power η = 1/3 (corresponding to the exponent predicted by the KTHNY theory [26][27][28]).
Summarizing, we established the full phase diagram of Active Brownian Hard-Disks, with active liquid, hexatic and solid phases, as well as coexistence and MIPS.
First, we proved that the overall scenario of 2D melting of passive disks is maintained for small-enough Pe. Weak activity acts as a perturbation on the passive case that destabilizes passive order, similarly to what was found in [45] for a system of softer disks (no coexistence in the passive limit) evolved with Monte Carlo dynamics. This is shown by the fact that increasing Pe both the liquidhexatic and hexatic-solid transitions shift to higher densities, and the liquid-hexatic coexistence region shrinks and eventually disappears. Such behavior can be due to the effective softness introduced by activity (quantified by the ratio between the active and potential forces Γ = ε/(σ d F act )), since, in equilibrium, particle softness reduces the liquid-hexatic coexistence region and eventually distroys it, rendering the hexatic-liquid transition continuous [3].
At high Pe, the MIPS region opens up on top of the hexatic and solid transition lines (differently from what was shown in [45]) and prevails the emergence of hexatic and solid order. In most of the MIPS region, many finitesize patches with different hexatic order coexist at any moment, but the large activity makes them regularly rearrange via breaking and recombination, very differently from what happens at low Pe. Above the point at which the hexatic transition line crosses the MIPS binodal, activity triggers phase separation between a low-density gas and a high-density hexatic, or solid, at higher Pe.
The discontinuity between the coexistence regions for Active Brownian Disks is distinct from what was found for active dumbbells, for which the large Pe phase separation was continuously connected to the zero Pe one. This difference could be due to the fact that dumbbells have a non-convex geometry that eases jamming and the formation of local orientational order. It would be interesting to study systems made of elements that interpolate between the disk and dumbbell geometries, and see how the topology of the phase diagram transforms from the one in Fig. 1 to the one in [22].
To conclude, our results provide a firm basis to rationalize the phase behavior of dense active matter and understand how self-propulsion affects the liquid and solid phases of matter on general grounds. The scenario we established here could be experimentally tested in, for instance, monolayers of self-propelled Janus colloids. This work was possible thanks to the access to the MareNostrum Supercomputer at  Supplementary Material for "Full phase diagram of active Brownian disks: from melting to motility-induced phase separation" In this Supplemental Material (SM) we display further evidence for the various phases and transitions explained in the main text. We organize the SM in two sections in which we expand the analysis of four observables: the local density in Sec. S1, and the local hexatic order parameter, hexatic correlation functions and global hexatic order parameter in Sec. S2. We discuss their behavior and implications.

S1. THE LOCAL SURFACE FRACTION
We sampled the local surface fraction φ i in the following way. We first divided the system in square boxes (or cells) of linear size σ. We then calculated a coarse-grained local density φ i associated to each cell i by computing the mean density over a circle of radius R centered at each cell. From the statistics of these local φ i values we constructed a probability distribution function (PDF). The choice of the specific value of R depends on the point in the Pe-φ plane under consideration; indeed, we adapted the coarse-graining to the heterogeneities of the system. Except for the simulations at Pe = 1, for which we used R = 20, the choice R = 5 was used in all other cases. Representative local density PDFs are shown in Figs. S1 and S2, where data for Pe = 1, Pe = 10, Pe = 50 and Pe = 200 are plotted. The first case corresponds to the low activity limit and the global packing fractions are chosen so that they lie within the black coexistence region in the phase diagram in Fig. 1 of the main text. The second case corresponds to the intermediate Pe region where no coexistence is observed, and the packing fractions are chosen, in particular, in the vicinity of the liquid-hexatic transition. The latter two cases are, instead, beyond the critical point towards MIPS.
The local density PDFs provide further evidence for phase coexistence both in the low-Pe (hexatic-liquid coexistence) and high-Pe (MIPS) regimes (see the phase diagram in Fig. 1 in the main text). In the low Pe limit, Fig. S1 (a), we see that the three curves displayed, corresponding to global densities varying over the very narrow interval [0.715 : 0.72], are bimodal with the weight under the two peaks slowly transferring from the one at low φ i to the one at high φ i for increasing φ. For higher Pe values we plot curves for a larger range of variation of φ. No two peak structure is found beyond Pe 3, that is to say, beyond the ending point on the black region of the phase diagram at low Pe values, and before the critical point for MIPS, the second black region in the phase diagram, is reached. Figure S1   The open circular symbols in the phase diagram in Fig. 1 were obtained from the location of the global packing fractions that limit the region with a double peaked PDF of local densities, in a very satisfactory agreement with the data stemming from the pressure measurements shown with filled black dots.
We note that from local density measurements we cannot investigate whether the dense phase has orientational or positional order. We addressed this point with local hexatic order parameter computations and position correlation functions, respectively.

S2. THE HEXATIC ORDER
As discussed in the main text, we calculated the local hexatic order parameter ψ 6,j = ψ 6 (r j ) = N −1 j Nj k=1 e iθ jk (S1) for each disk in the system, where N j are the first neighbours of the disk j. In order to do so, we built up the nearest-neighbors network by means of a Voronoi tessellation. We studied the local orientational order in the system constructing maps of this order parameter and its PDF. With this quantity we also computed correlation functions and we defined a global order parameter. We discuss all these measurements in this Section.

A. Maps of local hexatic order
In the following figures we visualize the local hexatic order using the method proposed in [1]: first, we project the complex local values ψ 6,j onto the direction of the mean orientation N −1 N i=1 ψ 6,i , where the sum runs over all particles in the sample. We then associate a color code to each bead according to this normalized projection. Regions with orientational order have uniform color. The dominating orientational ordering is painted in dark red, and the hierarchy follows the scale shown at the extreme right of the panels in Fig. S3, S4, S5 and S6. Due to the six-fold symmetry of the ordered state, blue regions are hexatically ordered along a lattice which is rotated by π/2 with respect to the one of the dark red regions. Green spots, which correspond to zero local hexatic order parameter in the color code, represent regions rotated by π/4 from the red ones. The information stemming from the PDFs of local densities in Figs. S1 and S2 is complemented by Figs. S3, S4, S5 and S6, where we show some representative snapshots of the system illustrating the nature of the different regimes reported in the main text. The snapshots encode the maps of local hexatic order, as explained in the previous paragraph. All the pictures correspond to the state of the system after letting it relax for about 10 6 MDs from a fully random initial condition.
In the low Pe case displayed in Fig. S3, where the three panels span the coexistence region of the phase diagram, we see large (non compact) regions with red color that correspond to the same local hexatic order that are surrounded in a rather disordered way by regions with no global hexatic order. These features are very similar to the ones seen in the active dumbbell system studied in [2] and in the passive disk models with sufficiently hard repulsive potential studied by Krauth and collaborators [1,3,4]. In Fig. S4 we display maps at Pe = 10, an activity for which we do not see coexistence, and the transition pattern is simpler with one line separating active liquid and active hexatic phases and another one separating the active hexatic from the active solid phase. To start with, the color map in panel (a) looks different from the ones shown in the three panels in Fig. S3, with no large red zone in this configuration, that corresponds to a low global packing fraction (φ = 0.5) and that we interpret as an active liquid one. The intermediate panel (b), obtained for φ = 0.820, displays an almost uniform, relatively light, red pattern and it lies in the active hexatic phase. Finally, a much darker uniform red state is shown in the last panel (c), in which φ = 0.860, and the system is an active solid, as confirmed by the analysis of the positional correlation functions shown in the main text.
The case Pe = 50 is above the critical Pe value for which MIPS occurs. In the phase diagram in Fig. 1 in the main text, one can see a reentrance of the active liquid above the MIPS region before entering the hexatic phase at an even higher global packing fraction. Typical color maps of the local hexatic order parameter are shown in Fig. S5 for global densities that lie well within the MIPS region (a), in the reentrant active liquid (b) and in the active hexatic (c). It is hard to assess from the snapshot in panel (b) what is the nature of the system for these parameters, although it is clear that there is no predominant hexatic order in the sample. The conclusion about reentrance was drawn from the analysis of other observables, notably, the pressure and correlation functions. At Pe = 200, the case considered in Fig. S6, MIPS generates the de-mixing of the system into a very low-density (φ low ≈ 0.07) and a rather high-density (φ high ≈ 0.9) phase with, on top, local hexatic order of different kinds. This is illustrated in the snapshots in panels (a) and (b) by the dense regions of particles sharing the same color, in coexistence with very low density regions. As largely reported in the MIPS literature [5], the dense phase induced by activity is subjected to anomalously large density fluctuations: it continuously breaks and reforms, giving rise to the observed distribution of orientationally ordered patches of finite size (instead of a single uniform one). We confirmed this claim by following the time evolution of these states (not shown). Panel (a) is also shown as an inset in Fig. 5 in the main text. Next to it, another inset explains the pattern of domains with different orientational order: the map of |ψ 6,i | is shown with a red scale and interfaces between domains with different ψ 6,i are highlighted in green. For high enough global packing fraction, φ = 0.860, the phase separation is between an active gas and an active solid, see panel (c) with a clear hole in the upper right corner.

B. PDFs of the modulus of the local hexatic order parameter
We now make the analysis of the local hexatic order quantitative by tracing the PDFs of the modulus of the corresponding local order parameter.
At low Pe, Fig. S7, the PDF has a bimodal form, in full correspondence with the one for the local density shown in Fig. S1. The peak at a high value of |ψ 6,j | corresponds to regions with local hexatic order while the one at low value of |ψ 6,j | represents the disordered regions. (We note that a disordered liquid phase does not have a vanishing modulus of the local hexatic order parameter but a unimodal distribution with average and typical values that increase with the packing fraction.) Increasing Pe beyond the end of the coexistence region connected to the passive limit (end of the black region in the phase diagram connected to Pe = 0), the PDF of |ψ 6,j | dramatically changes form, and has only one peak that displaces with increasing packing fraction from a low to a high value of |ψ 6,j | but never develops a double peak structure. These facts are shown in Fig. S7. Qualitatively, the form of the PDF is the same as the one for the local density shown in Fig. S1 (b). Therefore, from this measurement we cannot locate the phase transition towards the phase with hexatic order. This has to be done from the analysis of the orientational correlation functions and the Binder parameter, as explained in the main text. Finally, inside the MIPS region, the PDF of the modulus of the local hexatic order parameter acquires, again, a multi peak structure, see Fig. S8. The one at low values of |ψ 6,j | corresponds to the disordered and in some cases almost empty regions. Note that, basically, the position of the maximum of this peak does not depend on φ within the scale of this figure and it coincides with the one of the homogeneous phase at very low φ. The peak at |ψ 6,j | close to one indicates that there are regions in the system with almost perfect local hexatic order. Below this peak appears a second one, of lesser height, that is associated to the disks that are close to the interfaces of the perfectly ordered domains (of finite size). This fact is proven by the map in the second insert in panel (b) of Fig. 5 in the main text. This figure is to be compared to Fig. S2 where the PDFs of local density for the same parameters are shown. See the text for a discussion of the form of these curves, especially the secondary peak appearing below the one close to |ψ6,j| = 1.

C. Hexatic correlation function
Data for the correlation function of the local hexatic order parameter measured a distant points in space, g 6 (r) = ψ * 6,j ψ 6,k |rj −r k |=r / |ψ 6,j | 2 , at Pe = 10 and various global packing fractions φ, and various Pe and φ = 0.8, were shown in Fig. 3 in the main text. The change from exponential to algebraic was used as a criterium to locate the transition from the active liquid to the active hexatic phases. Here we exhibit the behavior of these correlation functions at a higher value of Pe where the change operates within the MIPS region and informs us about the nature of the dense phase that has either short, quasi or proper long range hexatic order.
In Fig. S9 we plot the correlations of the hexatic order parameter at Pe = 50 in (a) and Pe = 200 in (b), and we follow its form for varying φ. Differently from what was shown for the coexistence region at low Pe, for most densities the correlations decay very fast, as exponentials, where MIPS is between a dilute and a dense phase but, clearly, the latter has no long-range nor quasi long-range orientational order. Two curves in both panels do not decay so fast: the black and green sets of data. Concerning the black curve in (a), it was obtained for a global density that falls, according to our measurements, in the reentrant active liquid phase and it is still exponential. The black curve in (b) is within the MIPS area but close to the line beyond which the separation involves a dense phase with hexatic order. The green curves instead correspond to a global density that is very close to the one of close packing and the system behaves as a solid in these length scales.

D. Global hexatic order
We used the local hexatic parameter to compute the global hexatic order parameter that we show in Fig. S10, as a function of φ for several Pe values and system sizes. In panel (a), where we used a system with N = 256 2 particles, for each Pe, the global hexatic order parameter displays a sharp increase around a value of φ that increases with increasing activity. This shows that activity shifts the emergence of hexatic order toward higher densities. The position of the critical curve (Pe, φ) separating the active liquid (Ψ = 0) from the other phases with various criteria yield values that are in agreement with the determination of the hexatic transition line from the orientational correlations and the Binder cumulant (see the main text). The curves in (a) retain roughly the same form for all Pe values until, say, Pe = 40 where we see that they detach from zero at small densities developing a kind of plateau, the height of which increases with increasing activity. The way in which the highest value Ψ = 1 is reached in the ramping part of the curves Ψ(φ) is also modified at high Pe. In panel (b) we have checked how do the curves for fixed Pe = 50 depend on the number of particles in the system. From the main panel one sees that the form of the curves gets closer to the one of the weaker activity systems for increasing system size. We have also indicated with vertical dotted lines with different color the important density values: the limits of MIPS in black, the transition to the hexatic phase (as obtained from the correlation functions of ψ 6,j ) in blue and the one to the solid (as obtained from the positional correlation functions) in orange. In the inset in panel (b) we study the infinite N limit of the plateau value. We make four choices of the packing fraction φ = 0.4, 0.5, 0.6, 0.7, at which we trace the value of Ψ against 1/N . Although we only have three system sizes to work with, and it is not possible to conclude on the finite N dependence and N → ∞ limit of the data beyond any doubt, a naïve comparison to the law 1/N 1/2 is quite acceptable. This suggests that, in the infinite size limit, the global order parameter vanishes in the MIPS region below the (blue dotted) line in the phase diagram that indicates the entrance into coexistence with hexatic order.