Onset of criticality in hyper-auxetic polymer networks

Against common sense, auxetic materials expand or contract perpendicularly when stretched or compressed, respectively, by uniaxial strain, being characterized by a negative Poisson’s ratio ν. The amount of deformation in response to the applied force can be at most equal to the imposed one, so that ν = − 1 is the lowest bound for the mechanical stability of solids, a condition here defined as “hyper-auxeticity”. In this work, we numerically show that ultra-low-crosslinked polymer networks under tension display hyper-auxetic behavior at a finite crosslinker concentration. At this point, the nearby mechanical instability triggers the onset of a critical-like transition between two states of different densities. This phenomenon displays similar features as well as important differences with respect to gas-liquid phase separation. Since our model is able to faithfully describe real-world hydrogels, the present results can be readily tested in laboratory experiments, paving the way to explore this unconventional phase behavior. Auxetic materials are characterized by a negative Poisson’s ratio - they become thicker when stretched. Here, authors reach the limit of auxetic behavior, known as hyper-auxetic, in polymer networks, hitting an unconventional mechanical critical point.

T he mechanical response of a material subjected to uniaxial strain in the direction orthogonal to the deformation is quantified via the Poisson's ratio ν, defined as the negative ratio between transverse and longitudinal deformation. For the most common three-dimensional materials ν is positive, so that these expand (contract) in response to a compressive (extensional) strain. This situation is schematically illustrated in Fig. 1(a). On the contrary, auxetic materials are characterized by negative values of ν, meaning that they become thicker perpendicularly to the deformation axis, as shown in Fig. 1(b). Auxetic behavior has been so far reported in a large variety of systems, including foams, polymers, fibers, tendons, and crystals [1][2][3][4][5][6][7] . Recently, a strong research interest has been devoted toward auxetic metamaterials in which the elastic properties can be tailored by geometrical design [8][9][10] or by pruning methods 11 . Besides geometrical reasons, a negative ν can also be obtained by exploiting critical behavior and phase transitions, as in the case of ferroelastic materials in the vicinity of the Curie point 12,13 .
Within linear elasticity theory 14 , the appearance of a negative ν can be related to a decrease of the bulk modulus K with respect to the shear modulus G, namely to an isotropic softening of the material. A vanishing K echoes the divergence of the isothermal compressibility occurring at a gas-liquid critical point. However, the presence of a finite shear modulus, as found in polymer networks, such as hydrogels, may induce a negative ν. Pioneering evidence of a negative Poisson's ratio has been reported for these systems close to the so-called Volume Phase Transition [15][16][17] : in this case, a variation in temperature changes the affinity of the polymer to the solvent, favoring monomer-monomer aggregation, in full analogy with the gas-liquid critical point, but with the additional constraint of network connectivity. Another thermodynamic parameter that influences the network properties without affecting monomeric interactions is pressure, or tension. Indeed, theoretical works have addressed the occurrence of auxeticity in two-dimensional models of networks under tension 18,19 .
It is important to notice that these studies have flourished about twenty years ago, but the interest in hydrogel and microgel networks has increased again in the last few years, thanks to advances in chemical and in silico synthesis. In particular, it became recently possible to tune the amount of branching points (crosslinkers) to yield ultra-low-crosslinked networks [20][21][22][23] . In parallel, numerical efforts have been able to realize fully-connected, disordered networks with arbitrary density and crosslinker concentrations 24,25 .
In this article, we numerically investigate the elastic properties of polymer networks under tension and demonstrate that auxeticity naturally emerges in the ultralow-crosslinked limit. Combining stress-strain and equilibrium simulations, we show that low-density polymeric networks exhibit a nonmonotonic behavior of K, as well as ν, the latter becoming increasingly negative with reducing crosslinker concentration c. This phenomenology is found for both ordered diamond-like and disordered hydrogel realizations, indicating that there is no need of a specific topology to observe auxeticity in polymer networks. Remarkably, we do not find that this behavior continuously evolves down to c → 0. Rather, it hits a mechanical critical point that we name "hyperauxetic" point, at a finite c = c *~0 .35% where ν = − 1. Owing to the fact that K and G would become negative, this value denotes the lowest limit of mechanical stability, even though this condition is not yet fully understood 7,26,27 . At this point, we detect the onset of a coexistence between two different states: a low density and a high density one. These results call for an analogy with the well-known gas-liquid phase transition in attractive fluids, but with two important differences: (i) the lack of attraction between the monomers due to the good solvent conditions of the polymer networks and (ii) the presence of critical-like density fluctuations, which do not seem to obey Ising-like statistics within the present numerical resolution.

Results
Auxetic behavior of ordered and disordered hydrogels. We start by calculating the elastic properties of diamond networks (Diam) for different values of the crosslinker concentration c for negative pressures starting from P = 0. Although the diamond network is a simplified model 28 , it can still describe some relevant phenomena in experiments at a qualitative level 29,30 . For each studied state point, we independently evaluate three moduli: the bulk modulus K is obtained from equilibrium NPT runs, whereas the Young modulus Y and the Poisson's ratio ν are estimated from strainstress simulations, as described in the Methods section. We report ν in Fig. 1(c) for hydrogels with different c, while the corresponding K and Y are shown in Fig. S1 of the Supplementary Information. All moduli display a similar behavior: they initially decrease, then reach a minimum at intermediate values of pressures, and finally, increase again for very negative P. The minima become more and more pronounced for lower and lower c, remaining visible at all c for ν and K, while disappearing for Y when c ≳ 3%. Remarkably, we find that the different networks display a positive value of ν both for P = 0 and for very large negative pressures, while auxetic behavior is observed for c ≲ 5% in a finite range of tensions, that become smaller and closer to zero pressure as c decreases. We denote the minimum value reached by ν as ν min and its corresponding pressure as P min (see inset of Fig. 1).
So far we exclusively discussed ordered networks. However, in real-world realizations, hydrogels are intrinsically disordered and frequently made of chains whose length is exponentially distributed 31,32 . Being able to prepare hydrogels with these features, as described in Methods, we find that disordered networks (Dis) show the same phenomenology as ordered ones   Fig. 1 Auxetic behavior of ordered and disordered hydrogels. Illustration of (a) standard (ν > 0) versus (b) auxetic (ν < 0) behavior. Here, the red arrows indicate the uniaxial deformation (stretching) that is applied on the initial state of the system, represented by the dark cube that is identical in both cases. Following external strain, the system deforms perpendicularly in the directions indicated by the black arrows, leading to two different final states, represented by the light polyhedra; (c) calculated Poisson's ratio ν as a function of negative pressure P for different values of crosslinker concentration c = 1, 3, 5, 7.5% both for ordered (Diam, full symbols) and disordered (Dis, empty symbols) networks. Inset: zoom of c = 1%, where the minimum values of pressure and Poisson's ratio, P min and ν min , respectively, are indicated for the Diam-1% case. Pressure is given in units of k B T/σ 3 as described in Methods. Source data are provided as a Source Data file.
when subjected to tensions, with auxeticity also emerging for c ≲ 5% (see Fig. 1). Interestingly, we observe lower values of ν min for disordered hydrogels with respect to ordered ones at the same crosslinker concentration, as shown in the inset of Fig. 1. In particular, we find ν min ≈ − 0.8 for the Dis-1% network. We ascribe this effect to the higher structural heterogeneity characterizing disordered systems, independently of the specific network topology. Indeed, the same qualitative phenomenology is observed for all examined independent realizations (see Fig. S2). Since the preparation procedure is based on the self-assembly of patchy particles 24,25 , we cannot easily obtain disordered networks with smaller values of c, due to the nearby occurrence of phase separation. Thus, for the lower degree of crosslinking, we focus only on diamond networks, having shown in Fig. 1 that there is no major qualitative effect of the underlying topology on auxeticity, in line with previous results for ordered or partially ordered topologies [8][9][10][11] .
Hyperauxeticity and mechanical instability. By further decreasing c for Diam-N we still observe a progressive decrease of ν min . To have a better understanding of the mechanical behavior of the system, we report in Fig. 2(a) the negative transverse (λ ⊥ ) vs the longitudinal strain λ k À Á (see Methods for definition) for c = 0.35% at different values of pressure. We clearly see that, for P = 0, the slope, that is precisely our numerical estimate of ν, is positive. Then, by progressively reducing P, ν becomes negative, down to ν min ≃ − 1 within numerical uncertainty. At this point, the system has reached the limit of mechanical stability, so that state points with ν < − 1 are not allowed. We therefore consider c * ≃ 0.35% as the critical fraction of crosslinkers at which a mechanical critical point is encountered and define state points with ν = − 1 as hyper-auxetic. Interestingly, we find that, upon further increasing tension, the slope starts to increase again, as shown in Fig. 2(a). The behavior of ν vs P is reported for three ultra-low values of c respectively above, at, and below c * in Fig. 2(b). These findings indicate that the system reaches a hyperauxetic condition also for c < c * . This phenomenology thus occurs for ultra-low-crosslinked networks in general, with the system reaching ν min = − 1 at a small, finite negative pressure.
Critical-like nature of the transition. The behavior discussed so far close to hyperauxeticity may echo what happens close to a thermodynamic second-order phase transition, such as the gasliquid one, where a homogeneous system becomes thermodynamically unstable due to the divergence of the isothermal compressibility. To avoid this, the system thus separates into two phases characterized by a different density. It is now interesting to investigate by which mechanism ultralow-crosslinked hydrogels deal with the presence of the mechanical instability and which similarities or differences with respect to the gas-liquid scenario occur. To this aim, we report the behavior of the density fluctuations with respect to time for the Diam -0.35% system at P = P min in Fig. 3(a), detecting the onset of critical-like fluctuations. It is evident that, close to the mechanical instability, the system fluctuates between two states, an expanded and a compressed one, respectively illustrated in the corresponding snapshots of Fig. 3(b, c). We thus separately calculate the elastic moduli of the two states and find that the bulk modulus is much smaller in the expanded case (K EXPANDED~1 .4 × 10 −6 k B T/σ 3 ) as compared to the compressed one (K COMPRESSED~5 .0 × 10 −6 k B T/σ 3 ). A similar behavior was also detected for the Young modulus with the compressed state having it significantly larger than the expanded one. However, the Poisson's ratio does not change so much, being ν ≃ − 1 for the expanded state and − 0.84 for the compressed one, as shown in Fig. S3. These data also confirm the absence of relevant anisotropic effects or preferential directions within our system and suggest that a hyper-auxetic behavior is found in both states. Hence, the mechanism of a density jump is the one allowing the system to avoid the mechanical instability in full analogy with gas-liquid phase separation. Such an analogy is evident when looking at the density fluctuations as a function of time for different pressures around P min for the Diam -0.35% network that are shown in Fig. 4(a).
We now examine what happens as a function of c and plot the behavior of the (negative) pressure against density in Fig. 4(b) for all studied diamond networks. It is important to note that, in the present model, the crosslinker concentration plays the role of temperature in phase-separating fluids: as c becomes smaller, P becomes progressively flatter and a critical-like point is observed for these putative equation of states, indicated by a red diamond in the figure, which depend on a geometric rather than a thermodynamic control variable. For c ≤ c * , a clear discontinuity in density is observed, signaling a first-order-like transition between the two states. To get better microscopic insights of this behavior, we report the average bond length l bond as a function of P in Fig. 4(c). From a geometrical perspective, lowering the pressure towards more negative values has the effect of stretching the chains, as demonstrated by the growth of l bond . However, it is important to note that this behavior arises only when P ≈ P min , i.e. when ν starts to increase again beyond its minimum value. Since it is well-known that entropy plays a fundamental role for phase behavior of polymer systems 33 , we also quantify the effect of entropy by calculating the average single-chain entropy s l within a Langevin-approximation (see Methods and Ref. 25 ). Relying on this controlled approximation, we detect a decrease of s l by roughly one order of magnitude going from the compressed to the expanded state, as shown in Fig. 4(d). The entropy further shows critical-like fluctuations close to P min , as shown in Fig. S4.
On the other hand, when we monitor the average total potential energy e t , we find no remarkable change with pressure and, importantly, no critical fluctuations, as shown in Fig. S5(a). We thus focus on the non-bonded potential energy e nb , pertinent only to non-bonded particles, shown in Fig. S5(b), which instead manifests critical-like fluctuations. This correspondence is confirmed by the scatter plots, reported in Fig. S5(c) and (d), of total energy and nonbonded energy with density, respectively. Clearly, correlation is completely absent between e t and ρ, while e nb is correlated with density, similarly to the single-chain entropy as discussed in the SI.
Density distributions and comparison with Ising statistics. We now examine in more detail the nature of the density fluctuations and report the distribution of the density PðρÞ in Fig. 5(a) for different negative pressures close to the onset of the mechanical instability for the Diam -0.35% network. We notice that, close to the transition, the system displays a bimodal distribution, that is highly asymmetric. This is even more evident from the fact that data are reported on a log-log scale to improve visualization. In particular, we find the presence of a broad high-density peak and a narrow low-density one. Aiming to build a correspondence with thermodynamic phase separation, we next calculate the order parameter M, equivalent to the one used to describe the gas-liquid transition, that is composed by density and energy fluctuations 34,35 . Since we found that for the present system the total potential energy is not correlated with density, we define M = ρ + se nb , where we only consider the non-bonded potential energy and s is a mixing parameter 34 . In Fig. 5(b), we plot the distribution of the order parameter PðMÞ with zero average and unit variance for P = P min and different values of s. We find that the presence of the mixing term is able to reduce the asymmetry of the original PðρÞ, but still not completely. In particular, the heights of the two peaks become comparable for s = 0.9, but the difference in their variance is retained at all studied s.
In order to compare with the expected 3D Ising universal distribution, we apply the single histogram reweighting technique, as discussed in the Methods. This is shown in Fig. 5(c) for state points close to the mechanical instability at two different values of c. The resulting PðMÞ for ultra-low-crosslinked hydrogels are characterized by a slightly asymmetric shape, not perfectly centered in zero with rather broad peaks. We thus find a qualitative disagreement with the Ising behavior, independently of c, that might be due either to an intrinsic difference of the network system with respect to associating molecules or to insufficient sampling. Assuming true the first case, we may speculate that the presence of the network connectivity may bias the way in which the density fluctuates as compared to unbound particles. Alternatively, the deviation may be attributed to the fact that, in the present simulations, c cannot be varied in a continuous way, as normally done with temperature close to the gas-liquid critical point, thus hindering a proper exploration of the critical properties of the transition. However, we note that we found deviations from the Ising expectations for all examined c values, moving either below or above c * , finding no systematic improvement. Future work will be needed to properly address this issue.

Discussion
The present results, obtained by means of extensive simulations to calculate the elastic properties of ultra-low crosslinked polymer networks at negative pressures, report the emergence of auxetic behavior for c ≲ 5% independently of the specific geometry of the network. It is important to note that, when we started our investigations at not-too-small values of c and detected the onset of auxeticity, we expected to find a continuous behavior of the system until approaching the limit of stability (ν = − 1) at P = 0. This was indeed suggested by looking at the c-dependence of P min and ν min , that is reported in Fig. S6: while ν min decreases logarithmically, P min is found to obey a quadratic power-law behavior for all examined networks which spans a range of more than three decades in c. Instead, surprisingly, we found that the limit of the mechanical stability of the system, and the associated hyperauxetic behavior (ν = − 1), occurs at a finite crosslinker concentration, c * ≃ 0.35%, even for a regular network such as the diamond one. Hence, for even smaller values of c, the network, being unstable, undergoes a transition between two states, a compressed and an expanded one. Such a transition is also accompanied by large density fluctuations reminiscent of critical ones in gas-liquid phase separation, as shown in Fig. 4(a).
We then examined in more detail the nature of this phase transition, that is clearly distinct from the widely investigated Volume Phase Transition (VPT) of thermoresponsive hydrogels. Indeed, the latter occurs due to the change of the underlying polymer-solvent interactions at a characteristic temperatures. Previous studies on the VPT of hydrogels have already focused on the associated critical properties, that were tentatively attributed to the Ising university class [36][37][38] . Instead, the present work focuses on networks in good solvent conditions, where the monomer affinity to the solvent does not change and the underlying interactions are always dominated by excluded volume. Thus, the phase transition observed in the present work is the consequence of changing the network connectivity down to very low c, which generates a non-trivial interplay with steric interactions under a small tension. To this aim, it is instructive to focus on the values of the pressure at which ν reaches its minimum for c * , i.e., P min~1 0 −4 k B T/σ 3 , while the system volume V fluctuates around~10 6 − 10 7 σ 3 . This implies that the product PV is roughly comparable to the scale of the non-bonded particle energy~0.02Nk B T and about 3 − 4 orders of magnitude smaller than the total energy scale of the system~20Nk B T (see Fig. S5), where N ≈ 10 4 is the total number of monomers in our simulations. These considerations confirm that the present phase transition is dominated by fluctuations of non-bonded energy and of entropy, with a negligible contribution of the total energy. The important role of (infinite) connectivity and the different interactions with respect to a standard attractive system (e.g. a Lennard-Jones fluid) may thus be the reasons why the critical-like fluctuations of the present ultra-low-crosslinked hydrogels are not found to obey Ising universality class. Further numerical and theoretical work on this issue will be needed in the future. While the former should aim, in particular, to probe the critical fluctuations in a more extensive time and length window as well as to vary c in a more continuous fashion (e.g. by developing appropriate crosslinker insertion/deletion methods), the latter should be devoted to provide an additional description of the mechanical instability, taking into account the connectivity, similarly to what discussed for the Volume Phase Transition 39 .
Finally, it is important to note that these observations are relevant for ultra-low-crosslinked polymer networks, that are nowadays within experimental reach 20,40,41 . Indeed, for example, Poly(N-isopropylacrylamide) (PNIPAm) microgels are synthesized even in the absence of crosslinkers, taking advantage of (rare) self-association of NIPAM monomers, which gives rise to an effective c in the system that is close to zero 20,42 . It should thus not be difficult to realize this also for hydrogels. Given that experimental realizations are necessarily disordered, our numerical predictions suggest that disordered networks should display a slightly larger value of c * , not too far from 1% (see Fig. 1), with respect to the diamond case, which would actually favor the experimental observation of this mechanical critical point and a thorough exploration of its vicinity both from above and from below c * . Notably, our results are based on hydrogel simulations, but it would be interesting to apply our analysis also to ultra-lowcrosslinked microgels. In this respect, the microfluidic approach to microgels synthesis appears to be particularly promising, because it allows to prepare microgels of sizes of the order of 100μm 43,44 . Furthermore, a specific method to calculate their elastic properties, known as capillary micromechanics, was already established, making these ideal model systems to test our numerical predictions. We thus hope that the present results will stimulate novel experimental activity on ultra-low-crosslinked polymeric materials, aiming to verify their peculiar hyper-auxetic behavior and the occurrence of such an unusual phase transition, where mechanical and thermodynamic instabilities appear to be

Methods
Model. We perform Molecular Dynamics simulations of polymer networks made of monomers interacting via the Kremer-Grest potential. Excluded volume for all particles are given by the Weeks-Chandler-Andersen potential 45 : where σ is the monomer diameter, which sets the unit of length, and ϵ controls the energy scale. Defining m as the mass of the particles, the unit time of our simulations is defined as τ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi mσ 2 =ϵ p . Chemical bonds between connected monomers are modeled by a FENE potential V FENE r ð Þ 46 : where k F = 15 is the spring constant and R 0 = 1.5 is the maximum extension of the bond. We consider both ordered (diamond-like) and disordered topologies. In the former case, we prepare systems made up of 8 unit cells, each containing 8 crosslinkers, that are placed on the lattice atom positions and are connected through chains of equal length 47,48 . For each crosslinker concentration c, the network is composed of N = 64/c monomers forming chains of equal length l = (1 − c)/(2c), thus by varying l we change c in a controlled way. To produce disordered networks, we use the method recently developed in Ref. 24,49 , which exploits the self-assembly of binary mixtures of patchy particles with valence f = 2 (monomers) and f = 4 (crosslinkers). We let the system equilibrate at low enough temperature (T = 0.03) through the oxDNA simulation package 50 until 99.9% of the bonds are formed by exploiting a recently devised swap algorithm 51 . Then, we select the largest cluster from which we remove dangling ends and replace patchy interactions with the bead-spring ones (Eqs. (1) and (2)). We obtain systems with final crosslinker concentrations c~1, 3, 5% with deviation from the nominal values smaller than 5%. For c~1% we consider three independent network realizations to assess the dependence of results on the specific topology. For both ordered and disordered networks we perform NPT simulations using LAMMPS simulation package 52 with a Nosé -Hoover thermostat and barostat. Temperature is set to 1.0 throughout the manuscript and is measured in units of energy, i.e. fixing also k B = 1, where k B is the Boltzmann constant. We thus perform simulations at different (negative) pressures employing a timestep δt = 0.003τ.
Calculation of elastic moduli. We perform two kinds of simulations: (i) via equilibrium simulations in which the box is allowed to fluctuate anisotropically we calculate the bulk modulus K from volume fluctuations as K ¼ k B T hVi hV 2 iÀhVi 2 ; (ii) via stress-strain simulations we simultaneously calculate Y and ν. In particular, we first apply a longitudinal extensional strain λ k ¼ À L k À L 0 k Á =L 0 k , where L 0 k and L ∥ are the initial and final box lengths along the deformation axis, respectively. The range of deformation values encompasses λ k 2 0; 0:3 ½ , at which the response is in the linear regime, and we use a fixed strain rate _ λ ¼ 0:01τ À1 . The box is allowed to fluctuate transversally to the deformation in order to guarantee a constant average P. Then, once the system acquires the desired strain, the stress σ ∥ along the deformation axis is calculated from the virial stress tensor and averaged over 10 6 τ, yielding Y = σ ∥ /λ ∥ . The Poisson's ratio is instead extracted from transversal fluctuations through the relation ν = − ∂λ ⊥ /∂λ ∥ , where λ ? λ 2 þ λ 3 À Á =2 and λ 2 , λ 3 are the components of the strain orthogonal to the deformation axis. For each network and each state point, results for Y and ν are averaged over 20 independent deformations in which the same configuration is deformed with different initial velocities taken from the Maxwell-Boltzmann distribution. This procedure is repeated for each configuration by deforming the network over all three directions independently, which are then averaged in order to improve the statistics of the results. Results for K and Y are given in units of k B T/σ 3 .
Langevin approximation for single-chain entropy. In general, the entropy of a chain with n bonds of length b can be written as 53 s n; r ð Þ¼k B ln W n r ð Þ þ A n ; ð3Þ where r ¼ À r x ; r y ; r z Á is the end-to-end vector of the chain, W n r ð Þ is the end-to-end probability density and A n is a temperature-dependent parameter that can be set to zero. In the limit of systems that are submitted to a very large deformation or in a dilute regime, i.e., for r~nb, the end-to-end probability reads as 54 where β = 1/k B T and L À1 r=nb À Á is the inverse Langevin function, defined as L r=nb À Á ¼ coth r=nb À Á À nb=r 55 . Thus, the entropy of a single chain can be expressed as: where, in our case, b is the minimum value of the FENE interaction potential. We use this equation to calculate the single chain entropy of each chain in the network and then average over all chains to obtain the average chain entropy s l that is reported in Fig. 4(d) and in the SI.
Ising Universality class. A universality class groups phenomena arising in different physical systems, that, although being described by diverse microscopic models, exhibit an asymptotic large-scale limit that is characterized by the same invariant critical exponents 56 . In particular, the universality class defined by the three-dimensional Ising model describes the phenomenology of second-order phase transitions as diverse as the ferromagnetic Curie point or the gas-liquid criticality.
We thus attempt to compare the critical behavior of ultra-low-crosslinked networks to the expected behavior of a system belonging to the 3D Ising universality class, relying on the asymptotic expression of the probability distribution of the order parameter P Ising M ð Þ that can be conveniently approximated as 57 : where a = 0.158, c = 0.776, and γ is adjusted to provide unit variance to the distribution.
Histogram reweighting technique. We employ histogram reweighting in order to better identify the putative critical point of our phase transition, since this technique provides a powerful tool to reconstruct the probability distribution of a given observable at a state point P P 0 ; c 0 ð Þfrom equilibrium distributions of close enough state points, as long as thermodynamical control variables vary continuously. This is not the case for our system, for which c assumes discrete values in the ordered system and cannot be finely controlled in the disordered system, as a result of the self-assembly procedure of the network. We are thus forced to perform histogram reweighting only in P and, to this aim, we perform numerous, long-time NPT simulations at each c around P min . During the simulations, we record the behavior of the density ρ and of the non-bonded particle energy e nb as a function of time. Then, our single histogram reweighting method, reported in Fig. 5(c), relies on the following expression: P V; e nb ; P 0 À Á ¼ P V; e nb ; P À Á exp P À P 0 ð Þ V ½ : Thus, for each P 0 , we calculate the histogram reweighting factor exp P À P 0 ð Þ V ½ . We use this expression to obtain the distribution of the order parameter M = ρ + se nb , as discussed in the main text. This is calculated for all the values of ρ and e nb in P V; e nb ; P 0 À Á by varying the mixing parameter s. The set of M values is then rescaled to have a zero mean and a unit variance and compiled into the histogram P M ð Þ. This last is finally scaled onto the Ising curve for each given c by minimizing the mean squared error between the two distribution by varying P 0 .

Data availability
Data supporting the findings of this manuscript are available from the corresponding author upon reasonable request. Source data are provided with this paper.

Code availability
The simulation code is available from the corresponding author upon reasonable request. Most simulations reported in this work have been obtained using a precompiled code properly referenced in the manuscript.
Received: 29 July 2021; Accepted: 21 December 2021; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.