Aurélien Perera

Tomaz

Tomaz Urbic

Clustering in complex ionic liquids in two dimensions

HAL is

Introduction

The Coulomb interaction potential in three dimensions (3D) and two dimensions (2D) is the solution of the Poisson equation in their respective dimensions [START_REF] Chaikin | Principles of condensed matter physics[END_REF][START_REF] Caillol | [END_REF]. This interaction differs in each dimensions in an interesting manner, since in 3D it has a 1/r form, while in 2D it has a -ln(r) form [START_REF] Chaikin | Principles of condensed matter physics[END_REF]. Despite this fundamental difference, we have recently shown [3] that both systems have similar structure in the fluid phase, namely the charge ordering property [4]. Charge order is principally a short range feature, which enforces a form of order where charges of opposite sign are disposed in quasi-alternance [5]. In 2D, it leads to a checker-board appearance of the local order [3], due to optimizing of attraction/repulsion of the pair interactions. Even in the low density regime, where clustering dominates the structure of the system, the clusters obey charge order. In the same recent paper [3], we have shown that, replacing the 2D Coulomb interaction by the 3D screened Coulomb form, retained a strikingly similar local charge order, both in the liquid and gas phases. This finding enforces the idea that local order is mainly ruled by the strong pair interaction, despite the very different forms and signs of the full respective Coulomb interactions.

It is with this idea in mind that we wish to study how charge ordering in 2D is affected by the presence of neutral sites attached to the cation. This is motivated by the so-called room temperature ionic liquids (RTIL), such as ethylammonium nitrate, for example, which are liquid at room temperature [6] when ionic system such as NaCl are crystalline. We have previously argued [7] that it is the presence of inert sites attached to one of the charged atoms, which allows the system to prefer the liquid state at low temperatures instead of the solid phase, precisely because the inerts sites hinder the charge order, which otherwise would induce a crystalline state. In particular, this hindered charge order produces local clustering of the charged segments, which, in turn, produces a scattering pre-peak in the structure factors. These effects are well documented for realistic RTIL in 3D [8,9,10,11,12,13,14]. In the present work, we will consider the 3D form of screened Coulomb interaction, but for the strict 2D case, similar to that in our previous study in Ref. [3]. For the 2D case, it would be interesting to know how this screening of the charge order, induced by the presence of neutral sites, is affected, when we expect lesser possibilities of molecular conformations? To this effect, we study here the influence of thermal disorder at various temperatures, but also at various densities, in order to see the influence of clustering. In addition, there are 2 other interesting issues. The first issue concerns the existence of a liquid-gas coexistence at low temperatures, which has been intensely studied for case of the simple monomer system the so-called restricted primitive model (RPM) [15,16,17,18,19,20] and its 2D version [21,22,23,24,25,26,27,28,29,30,31]. The second issue concerns the existence of a Kosterlitz-Thouless (KT) transition [32] in the very low density region, where no free charges exist below the KT horizontal line in the (density, temperature) phase diagram [26,28,29]. We examine here how these 2 properties are affected by the presence of neutral sites.

In the present study, we principally focus on the model illustrated in Fig. 1, namely a single site anion and linear molecular cation with sites tangently attached to each other, with the cation site at on end. Our study clearly indicates that 3D to 2D dimensional reduction hinders the stability of the liquid phase at lower temperatures, such that the molecular ionic system has a stable liquid state for temperatures higher than the simple monomer ionic liquid. This is exactly the opposite behaviour than in 3D. In other words, molecular complexification in 2D restricts the range of the liquid phase. Interestingly, layer-like clustering is favoured at low densities, enforcing a horizontal asymptote in the (density, temperature) phase diagram for the stability of the disordered phase with respect to the cluster "phase". Since this cluster "phase" needs to be specified in relation to the behaviour of the integral equation theory, we can only conjecture about the underlying KT type behaviour in relation to the layer-like association which occurs in this part of the phase diagram.

Models and technical details

In a previous paper [3], we have considered charged system of monomers of equal size, with the 3D form of screened Coulomb interaction. In the present paper, as illustrated in Fig. 1, we consider a monomeric anion together with a dimer made of tangent spheres, one being a cation and the other neutral. All sites are taken to be spheres of diameter σ. The total site-site interaction reads

βv ij (r) = Z i Z j T C T exp(-r/λ) r/σ + 4T 0 T σ r 12 (1)
where β = 1/k B T is the Boltzmann factor, with T the temperature expressed in Kelvin, T C = 55700K corresponds to the temperature in the 3D Coulomb interaction [7] for the choice σ = 3Å, T 0 = 100K is arbitrarily chosen. The screening parameter is chosen to be λ = 2. All these parameters are the same as in Ref. [3]. The valences are Z 1 = -1, Z 2 = +1 and Z 3 = 0. Unlike our previous study [3], where we have used reduced temperatures, we will use here temperatures in Kelvin, which allows to make contact with our previous work in 3D.

It is important to note that, since we use screened Coulomb interactions it would be incorrect to conclude that the influence of charge in the short range order is minor. Indeed, in Ref. [3], we demonstrated that the structure of the system, as witnessed by the correlation functions, both in real and reciprocal space, show the same characteristics as when unscreened Coulomb is used. This remark is even more important since we are comparing the log-Coulomb of the 2D case with a screened version of the 3D case. The reason for this strong influence of the charges comes from the fact that the Coulomb interactions dominate the short range ordering through the factor T C T 0 . Interestingly, this effect is not only important in dense fluid, but also in the gas phase, as we show below in the Results section.

Monte Carlo simulations

All Monte Carlo (MC) simulations are done in the canonical (NVT) ensemble following the same protocol previously outlined in Ref. [3]. With one MC cycle consisting of a tentative move of N particles, 10 6 equilibration moves and 10 7 moves for statistics are performed for each system. Cut off of the potential was half-length of the simulation box. All simulations were performed with N = 100 or N = 200 molecules. Increasing the number of particles had no significant effect on the calculated quantities. The fact that smooth correlation functions with very low noise level are obtained is strong indication that the simulations are well converged.

Integral equation theory

Concerning the integral equation theory (IET) approach, we have used the 2D site-site Ornstein-Zernike (SSOZ) formalism [33,[START_REF] Hansen | Theory of Simple Liquids[END_REF] together with the hypernetted chain (HNC) closure [START_REF] Hansen | Theory of Simple Liquids[END_REF]. The choice for this particular closure in place of others, such as the mean spherical approximation (MSA), self-consistent closures [26,28], or the Hirata-Kovalenko type closure [START_REF] Kovalenko | [END_REF], is that the HNC closure represents the first level of approximation where all correlations higher than rank 2 are neglected [START_REF] Hansen | Theory of Simple Liquids[END_REF]36].

The SSOZ equation consists in the following matrix equation:

SM = I (2)
where the total structure factor matrix S is given by

S = W + ρ 2 H (3)
and

M = W -1 - ρ 2 C (4)
The intramolecular part of the total structure factor is defined through the matrix W as

W =   1 0 0 0 1 J 0 (kσ) 0 J 0 (kσ) 1   (5)
where J 0 (x) is the zeroth-order integer Bessel function. ρ = N/V is the total density of the system. The matrices H and C with respective elements hij (k) and cij (k) are the 2D-Fourier transforms of the pair and direct correlation functions, h ij (r) = g ij (r) -1 and c ij (r), where g ij (r) is the radial distribution function between monomeric sites i and j. The 2D Fourier transform of a function f (r) is defined as

f (k) = 2π ˆ∞ 0 rdr f (r)J 0 (kr) (6)
Since the Coulomb interaction is short ranged, the 2D Fourier transform of the direct correlation functions are well defined at k = 0, and it is not necessary to take the special precautions described in Ref. [3] for unscreened Coulomb interactions.

The HNC closure equation is

g ij (r) = exp [-βv ij (r) + h ij (r) -c ij (r)] (7)
where the exponential term is missing the so-called bridge function b ij (r) , which contains all the high order rank correlations. Setting b ij (r) = 0 is the first level of controlled approximations, which is why this particular closure is interesting.

Other choices, such as the mean spherical approximation (MSA) [START_REF] Hansen | Theory of Simple Liquids[END_REF], represent even higher level of approximation. Further choices, such as self-consistent closures [26,28], or the Hirata-Kovalenko type closure [START_REF] Kovalenko | [END_REF], are uncontrolled approximations, mostly based in empirical methodologies, which could help improve finding numerical solutions when HNC cannot, but such methods cannot help understand nor appreciate the role played by high order correlation through b ij (r).

Both the HNC and the SSOZ equations are approximations. In particular, the SSOZ equation used here, has known deficiencies [37]. It is only through the comparison with simulation that one can assert the range of applicability of these 2 equations. This is the empirical approach that we use here.

These 2 equations (2, 7) are iteratively solved using standard techniques developed for the 2D case [38]. The correlation functions are sampled on a logarithmic grid of 1024 points, and the Fourier transforms are handled through the Talman technique [39,40].

The atom-atom structure factors shown in Section 3 are defined as

S ij (k) = 1 + ρ 2 hij (k) (8)
They are related to the structure factor defined in Eq(3) by removing the intramolecular part. Also we add 1 to the cross terms, instead of the usual δ ij , in order to facilitate the graphical representation.

While the simulations meet no problems even when strong clustering is present, we find that the IET cannot be solved below the no-solution line shown in Fig. 2. As mentioned in Ref. [3], in the case of Coulomb interactions, this behaviour does not appear to be due to the onset of a liquid-gas transition, but to a strong clustering of opposite charges. This clustering increases the first peak of the unlike ions correlations, which is one of the causes for the raise of the corresponding structure factor near k=0, in addition to the usual long range tail of the correlations. These points are discussed in the Results section below.

Results

Phase diagram

Fig. 2 shows the no-solution "phase" diagram as obtained by the HNC approximation. The data from the monomer ionic fluid of Ref. [3] is equally shown in dotted lines. It is seen that the, to the difference of a factor in density, which could match the difference in volume of the two types of system, the shape is nearly the same. The dimer model shows an increase at high density which could correspond to the existence of the solid phase at even higher density. What is more intriguing is the flat asymptote behaviour at very low densities, as shown in the inset with densities in log scale. This behaviour bears some resemblance with that reported by Lomba et al. [28] (inset of Fig. 1), in particular when comparing the behaviour of computer simulations from Ref. [18] with that of IET for the case of charged hard-discs model. This latter model mimics the behaviour of the true Coulomb 2D system, which undergoes a Kosterlitz-Thouless transition in the low density region, where pairs of opposite charges bind and the system becomes electrically neutral for lower temperatures [21,22,26]. In Ref. [3], our results suggested that the theory predicts something very similar, as previously noted [28]. The low density behaviour of the present system is even more striking with the KT-like behaviour. However, snapshots of simulations indicate that, in the vicinity of the no-solution line, despite strong clustering, there are several free charges, which can be explained by the fact that the interaction is screened, and particles far away cannot irreversibly form pairs. At lower temperatures, we observe that all charges form bilayer type clusters and no free charges remain. So, it is tempting to conclude that the approximate theory predicts this binding reminiscent of the KT ion pairing.

The working hypothesis of this report is that the no-solution line delimited by HNC corresponds to a physical line below which the simulations show that clustering occurs. The simulation themselves do not show any sort of singularity except for the existence of marked clustering below this line. This is very similar to what we reported in Ref. [3]. At present we have no specific characterization of what the "phase" below the no-solution line could represent, other than "cluster phase". As can be seen in the energy plot (Fig. 4 and Fig. 5 discussed below), there is no sign of thermodynamic singularity in the vicinity of the no-solution line.

Snapshots

Fig. 3a-d show typical snapshots of the system showing how clusters form below the no-solution line, for 3 different densities ranging from dense liquid ρ = 0.7 to gas phase ρ = 0.01. For each density, 3 temperatures are shown, a high temperature (above and close to the no-solution region), a temperature just below the no-solution region, and the lowest temperature T=500K we simulated. At this temperature we would expect a solid phase in principle. The contrast between the disordered behaviour above the no-solution line and the existence below it of well defined clusters with no or little free particles is obvious and striking. For the dense liquid, we observe that below the no-solution line charge and invert groups show micro-segregation. In 3D it is possible to find solutions with a segregated domain pre-peak in the structure factor, but not in 2D.

The existence of "droplets" for lower densities could suggest that the system is in a 2 phase region, implying a phase separation. However, we could not observe such phase coexistence between a presumed gas and a liquid. The fact that small pieces of bilayers are formed in the dense region indicates that no liquid phase is formed. It is more tempting to suggest the existence of a cluster phase rather.

Thermodynamics

Fig. 4 shows the excess energies per particle (upper panel) and the constant volume heat capacities (lower panel) versus temperature, as obtained from computer simulations. Each curve corresponds to a isochore. In relation to the no-solution diagram in Fig. 2, we have drawn the curves for densities below ρ = 0.1 in dotted lines, and full lines for densities above. The curve for ρ = 0.4 which corresponds to the minimum of the no-solution line is drawn in red. The two highest densities,ρ = 0.76 and ρ = 0.7 are shown in dotted cyan, since they have trends different from the other curves. The yellow dots represent the no-solution line in Fig. 2. For each density, the energy is seen to become more negative as more and more clusters form, which is expected. The heat capacity is seen to become more noisy in the low temperature cluster region, which is expected of the rather small size (N=200) in the simulations. In a way, the appearance of these fluctuations help delimitate the cluster region. However, we do not see any signature or singularity in the vicinity of the no-solution line. Indeed, the smoothness of the curves in the vicinity of each yellow dots clearly indicate that there is no variation of curvature, which would point to a change of energy due to sharper clustering, and would then reflect as a bump in the heat capacity. This is not observed here. It is easy to connect the no-solution points (the yellow dots) in the energy/temperature diagram into a u-shaped curve, but not in the heat capacity diagram, specially for the high density part. This is partly due to the fact that the heat capacities change as the minimum in the phase diagram (shown as a red line in the plots) is crossed.

Fig. 5 shows a comparison of the excess energies as obtained from simulations (in blue) and the theory (magenta), over the temperature range and for 3 typical densities ranging from high (ρ * = 0.7, upper set of curves) to low (ρ * = 0.1, lower set), through an intermediate density (ρ * = 0.4, middle set). We observe a rather good agreement over the range of overlapping temperatures, with a very similar trend of both calculations. The agreement is excellent for medium to low densities. The last low temperature point of the theory corresponds to the no-solution point, and it is clearly seen that there is no trend for any change of curvature, which would suggest a phase change. In agreement with this observation, we do not find any singularity in the heat capacities, neither in the compressibility. This is in contrast with what was observed in the case of an underlying gas-liquid transition in Ref. [28], where a sharp change in curvature of the compressibility and the heat capacity was reported in the vicinity of the no-solution points.

The conclusion we draw in this study is that usual thermodynamic quantities which help signal first or second order transitions, do not show any singularities or marked behaviour as the cluster line, or the no-solution of the IET, are crossed.

Cluster size distribution

Cluster size distributions have been calculated by computer simulations by using a distance criteria, namely the first minimum of the g +-(r) function. Fig. 6 shows these distributions, as the number of cluster N (c S ) of size c S versus the cluster size c S , and for typical 4 densities, ranging from high ρ = 0.76 to medium ρ = 0.4, low ρ = 0.1, to very low ρ = 0.01. For each density, the cluster size distributions for all the temperatures studied from T = 5000K (top curve in red) until T = 500K (lowest curve in green) are displayed. The temperatures corresponding to the no-solution line of Fig. 2 are shown in cyan, and those below this temperature are in grey lines. For all 4 cases, it is clearly seen that there is a gap separating the high temperature distributions from the lower ones. The gray curves show a different behaviour than the higher temperature distributions, with a clear tendency to display higher cluster probabilities for large sizes. For all 4 cases aside the lowest density, it is seen that the gap starts more or less below the cyan curve, which indicates a good agreement with the cluster "signature" given by the theory through the no-solution line of Fig. 2. This signature fails at the very low density, as the no-solution curve (at T = 3500K) is clearly above the gap signature, shown as the gold line at T = 1500K. This discrepancy between theory and simulation can be attributed to the fact that SSOZ equation is known to lead to incorrect predictions for the very low density region [37]. This means, that the cluster region in very low density region, as obtained from the simulations, should saturate at a lower temperature than what is shown from the theory. Otherwise, it would seem that the no-solution prediction of the HNC theory is in rather good agreement with predictions from simulations.

Correlation functions and structure factors

Fig. 7a-e show a comparison between the simulation and the HNC approximation for all the 6 site-site correlation functions (left panel) and corresponding structure factors (right panel), for different densities and temperatures very close above the no-solution line of Fig. 2, which are the most demanding conditions to test the theory. These plots principally illustrate the remarkable agreement between the simulation and IET data. For the dense liquid phase ρ = 0.7, we also show a plot for a high temperature (Fig. 7a), demonstrating that there is not much quantitative difference with the low temperature case in Fig. 7b. Fig. 7e shows a comparison for very low density ρ = 0.01 at T = 3500K, and more particularly an interesting sub-structure which appears for g +-(r) obtained by the IET. It concerns a marked double shoulder feature which appears just at base of the first peak around r ≈ 2σ. This feature is absent from simulation (magenta curve). However, we find a similar feature in the simulation data, but for a lower density ρ = 0.002. This is reported as a green curve in Fig. 7e. This feature gives an indirect indication in the clustering differences between theory and simulations. We believe that the HNC closure tends to exaggerate near neighbour correlations, as can be noted for the hard sphere fluid [START_REF] Hansen | Theory of Simple Liquids[END_REF], but also for orientational ordering [40]. Following this feature of the theory, and the slanting of the no-solution curve in the low density region of Fig. 2, it seems reasonable to suppose that the clustering feature found in the theory happens at a lower density for the same temperature. This finding points towards again towards clustering dominating the lower part of the phase diagram. For this low density ρ = 0.01 in the left panel of Fig. 7e, we has shown the correlation functions hij (k) instead of the S ij (k), since the low density ρ = 0.01 damps all the features because of the definition Eq.( 8). All these figures show minor disagreements here and there, but these can easily be accounted for the missing bridge term in Eq.( 7)

The figures also show that the very good agreement between the 2 approaches, holds both in real and reciprocal space. Such agreement was equally noticed in 3D in Ref. [7]. It was attributed to low fluctuations of the short range order and the subsequent homogeneity enforced by the strong charge ordering. Indeed, approximate IET tend to be less accurate when fluctuations are present, either as reflecting the great possibilities of positional order, such as in simple liquids, or long range correlations such as in the vicinity of phase transitions. In contrast, the strong local order imposed by the charge ordering reduces fluctuations and enforce the type of agreement we observe for all the region where solutions could be found. The fact that the worse agreement is precisely found for the correlations involving the neutral site X, which have more disorder in their positioning, further enforces the homogeneity argument presented here.

The fact that the agreement holds all the way until the no-solution line is hit, is a very strong indication about the nature of the state below the no-solution line. It suggests that this state is not due to some mechanical instability of the upper homogeneous phase. If it was, then we would see a progressive loss of agreement as we near the no-solution line from above. On the contrary, this second phase is due to clustering and not fluctuations, as illustrated by the snapshots in Fig. 3. The passage from homogeneous phase to cluster phase is not made through any thermodynamic signature, such as heat capacity or entropy, nor appearance of critical fluctuations.

The approximate IET cannot account for the clusters which appears below the no-solution line, because they miss high order correlation through the socalled bridge function. However, the rather good agreement with simulations found until the no-solution line, tend to indicate that these bridge diagrams do not play an important role until this line is met from above. It is possible that they become suddenly important below this line, hence explaining why IET cannot get there. This tentative explanation links the cluster "phase" to the raise in importance of high rank correlations, principally through the bridge diagram term. This points requires separate investigations.

Supra-molecular structures and pre-peak

In a previous study of 3D room temperature model ionic liquid by one of us [7], it was found that the presence of neutral sites induced a local segregation of charged and neutral sites, and in agreement with what is observed in realistic such liquids [10]. This segregation reflects itself in the presence of a low-k pre-peak feature, observed both in the cross charge structure factor S +-(k) and in like charge structure factors S ++ (k) and S --(k). By separating out the charge-charge and density-density structure factors through the Bhatia-Thornton transformation [START_REF] Bhatia | Fig.3d -Snapshots for very low density ρ = 0.01 at 3 different temperatures T = 3500K, T = 2500K and T = 500K . . .Fig.4 -Energy (top) and heat capacity (bottom) as a function of temperature (divided by 1000), as obtained from the simulations[END_REF], it was found that only the density-density structure factor retained this pre-peak [7], indicating that it is indeed related to hetereogeneity in the spatial density distribution.

Fig. 8 shows the Bhatia-Thornton (BT) structure factors of selected state points lying just above the the no-solution line, and corresponding to some of the structure factors shown in Fig. 7b-e. These BT structure factors are in fact related to a linear transformation from microscopic densities of charged atom to total density and charge density. The resulting structure factors are defined as in our previous work [7]:

S N N (k) = S ++ (k) + S --(k) + 2S +-(k) (9)
S ZZ (k) = 1 4 z 2 + S ++ (k) + z 2 -S --(k) + 2z + z -S +-(k)
S N N (k) represent the structure factor related to total density fluctuations, while charge fluctuations are represented by S ZZ (k). We observe again that the agreement between the simulation and IET data is very good. Similarly to Fig. 7e, and again because of the low density ρ = 0.01, the lower right panel shows the h(k) corresponding to the BT structure factor, with hcc (k) shifted by 1 in order to enforce the resemblance with the other plots. These plots indicate that only the case ρ = 0.4 shows marked pre-peak feature in S cc (k), whereas the case for ρ = 0.7 shows only a shoulder, and the low density cases show mostly k=0 density fluctuations. These finding are consistent with the snapshots shown in Fig. 3a-d. At high density, we observe a segregation of charged and neutral groups, but the dimensionality does not allow for a marked segregation as in the 3D case. The marked pre-peak for the medium density case is possibly due to the clear clustering (Fig. 3b), which enforces the local heterogeneity. In Fig. 3b, for T=1000K we observe clear chain-like clusters with evident +-chain formation. In other words, we see that domain segregation is strongly affected by dimensionality, and that it is stronger in 3D than in 2D. This finding could be of relevance for the observation of charge segregation in adsorbed realistic RTILs.

Discussion and Conclusion

Although we use a theoretical approach, with simulations and integral equations, this paper is similar to an experimental paper, in the sense that we present only outcome of calculations, which in turn suggest some properties of the system. In particular, we point to the existence of a cluster phase, which is not accompanied by any of the usual signatures for thermodynamical phase transitions. Since we compare approximate theory with simulations, we cannot attest that the nosolution line represents the actual boundary between the homogeneous and the cluster phase. If it was possible to detect this line from simulations alone, the evidence presented here suggests that it could possibly lie below the approximate no-solution line predicted by IET, but close to it.

Although our results concern only screened version of the 3D Coulomb interaction, we do not think that incorporating true 2D Coulomb interaction would modify the conclusions reached here. This is because it is the strong short range order is similar in both type of interactions, and we conjecture that it rules the structural properties of this type of systems. In particular, smectic layer-like clusters would most probably form even in the case of the log-Coulomb interaction.

Clustering plays an important role, even in simple liquids [41]. It is usually related to fluctuations, which concerns principally the k = 0 part of the structure factor, as far as the stability of the system is considered [START_REF] Hansen | Theory of Simple Liquids[END_REF]. Our experience in studying realistic 3D associating liquids, such as water or alcohols, indicates that fluctuations at k = 0, in addition to being related to clustering [42], play little or no role in the global stability of the system. On the contrary, they enhance stable local heterogeneity [42]. Therefore, since in the present case, simulations indicate that bilayer-like clustering appears in the lower part of the phase diagram, they support the fact that this system is governed by charge ordering induced clustering everywhere in the phase diagram, albeit to various degrees. In view of the remarkable agreement found in the correlation functions obtained from simulation and calculated from the theory, these findings help supporting the hypothesis that the no-solution line found in approximate IET could be a physical line distinguishing between different clustering regimes. The lower part of the phase diagram in Fig. 2 would the be dominated by many body correlations, which cannot be captured by two-body level description. This line of argument would also explain why such IET are unable to provide solutions for well mixed but micro-heterogeneous aqueous mixtures, which could equally require explicit many body correlation description. Subsequent investigation along these lines are in progress.

Finally, the polar-nonpolar domain segregation is found to be diminished by dimensional reduction. This is very apparent for dense surface coverage, but the segregation seems to be restored when particle confinement conditions are decreased by lowering the surface coverage density. This finding could have some relevance to 2D adsorption of realistic 3D ionic liquids.

Figure Captions

Fig. 1 Ionic liquid model, with 1 site anion (red) and dimer cation with 1 charged (blue) and 1 neutral (magenta). All sites have same diameter and same 1/r^12 dispersive repulsion (see text). .

Fig. 2 (

 2 Fig.2(Density, temperature) no-solution "phase" diagram from the IET. Yellow dots correspond to the lowest temperature for which IET could be solved, and the blue line is connecting them. The dotted line represent the no-solution line for Model 2 from Ref.[3]. The inset shows the same diagram but on log scale for the density.

Fig. 3a

 3a Fig.3aSnapshots for high density ρ = 0.7 at 3 different temperatures T = 2000K, T = 1000K and T = 500K. The anion is red, cation blue and attached neutral site is magenta.

Fig. 3b Fig. 4

 3b4 Fig.3b Snapshots for medium density ρ = 0.4 at 3 different temperatures T = 1000K, T = 800K and T = 500K Fig.3c Snapshots for low density ρ = 0.1 at 3 different temperatures T = 2500K, T = 1500K and T = 500KFig.3d Snapshots for very low density ρ = 0.01 at 3 different temperatures T = 3500K, T = 2500K and T = 500K Fig.4 Excess energy (top) and heat capacity (bottom) as a function of temperature (divided by 1000), as obtained from the simulations. Each line correspond to a density in the range ρ = 0.76, 0.70, 0.60, 0.50, 0.40, 0.30, 0.20, 0.10, 0.05, 0.02, 0.01, 0.005 and 0.002. Curves for very low densities below 0.1 are in dotted lines, as well as those for very high densities above 0.65 (cyan). The line in red is for ρ = 0.4, which corresponds to the minimum of the no-solution line of the IET in Fig2. Points corresponding to this no-solution line are indicated in orange dots. The green line connecting the orange dots is indicative (see text).

Fig. 5

 5 Fig.5Comparison of the excess energies between simulations (blue curves) and integral equations (magenta curves), as a function of temperatures, and for 3 densities, ρ * = 0.7,0.4 and 0.1. For clarity, the curves for ρ * = 0.7 and ρ * = 0.1 have been shifted by +5 and -5, respectively (the numbers shown in orange next to the labels).

Fig. 6

 6 Fig.6Cluster size distributions for 4 densities (as indicated above each panel) and temperatures (from top to bottom) T = 5000K (red curve), 3500K, 3200K, 3000K, 2800K, 2500K, 2000K, 1500K, 1000K, 900K, 800K, 700K, 600K and 500K (green curve). For each density, the cyan line corresponds to the temperature the closest to the nosolution line of Fig.2. All size distributions below this temperature

Fig. 7a

 7a Fig.7aCorrelation functions (left) and corresponding structure factors (right) for high density ρ = 0.7 and high temperature T = 3000K. X designates the neutral site of the cation in Fig.1.

Fig. 7b

 7b Fig.7bSame as Fig.7a, but for ρ = 0.7 and temperature T = 1500K closer to the no-solution line in Fig.2.

Fig. 7c

 7c Fig.7cSame as Fig.7a, but for medium density ρ = 0.4 and temperature T = 1001K Fig.7d Same as Fig.7a, but for low density ρ = 0.1 and temperature T = 2000K Fig.7e Same as Fig.7a, but for very low density ρ = 0.01 and temperature T = 3500K. The green curve is explained in the text.Note, that it is hij (k) that are plotted in the right panel (see text)

Fig. 8 .Fig. 1 -Fig. 5 -

 815 Fig.8 Bathia-Thornton structure factors S cc (k) and S zz (k) for the state points corresponding to Figs.7b-e. S cc (k) is shown in blue for IET and dotted green for simulations. S zz (k) is shown in red for IET and dotted gray for simulations. The lower right panel shows hcc and hzz (see text)..

Acknowledgments

The authors thank the partenariat Hubert Curien (PHC) from Campus France for financial support under the bilateral PROTEUS PHC project 35120VG.