Gelation, Clustering and Crowding in the Electrical Double Layer of Ionic Liquids

Understanding the bulk and interfacial properties of super-concentrated electrolytes, such as ionic liquids (ILs), has attracted significant attention lately for their promising applications in supercapacitors and batteries. Recently, McEldrew \textit{et al.} developed a theory for reversible ion associations in bulk ILs, which accounted for the formation of all possible Cayley tree clusters and a percolating ionic network (gel). Here we adopt and develop this approach to understand the associations of ILs in the electrical double layer at electrified interfaces. With increasing charge of the electrode, the theory predicts a transition from a regime dominated by a gelled or clustered state to a regime dominated by free, crowded ions. This transition from gelation to crowding is conceptually similar to the overscreening to crowding transition.

Indeed, owing to the high concentration and lack of high dielectric solvent, the energy of electrostatic interactions between nearest neighbours is much larger than thermal energy [3].
Sometimes, only free ions were explicitly treated, with ion pairs being implicitly treated as 'voids' [59][60][61][62]. This proved to be successful in reproducing the conductivity of ILs based on a Nernst-Einstein relation [62], where 10-20% of ions were free, with the remaining being bound up in immobile aggregates. A dynamic exchange between these two states occurs, with a small activation energy (∼ 1 k B T ) for the excitation from 'bound' state to the 'free' state. Thus the concept of an intrinsic narrow gap 'ionic semiconductor', as introduced in Ref. 3, has received substantiation [62]. That work showed a dynamic picture of the free ions, but the exact nature of the immobile clusters was not discerned. Moreover, free ion approaches were also able to explain the qualitative changes in differential capacitance as a function of temperature [60], where a 'melting' of 'frozen' structures is observed to occur as a function of increasing temperature and voltage. However, the fraction of free ions there was inferred from fitting the differential capacitance, and the nature of the 'frozen' states only discussed on a qualitative level.
Therefore, in such a concentrated system as ILs, the formation of clusters larger than ion pairs should occur [47,48,50,63]. McEldrew et al. [49,[64][65][66] developed a theory, based on the theories of thermoreversible polymers [67][68][69][70][71][72][73][74][75], which was able to systematically describe clusters beyond ion pairs, and even a percolating ionic network of infinite size [76][77][78], referred to as the gel. One of the main parameters of the theory is the number of associations an ion can make, where it is assumed cations can only bond to anions and vice versa, which determines the possible clusters that can be formed. If ions can only form one bond, then ion pairs can only occur. If a cation has four association sites, it can form various clusters; for clusters containing a single cation there can be free cations, neutral ion pairs, and negatively charged triples, quadruplets and quintuplets, as schematically shown in The theory of McEldrew et al. [64] was applied to bulk ILs [49], where a consistent theory of ionic transport was also developed based on vehicular motion of the clusters [49,79], and other concentrated electrolytes [65,66]. While the ionic association theory has given a detailed description of the cluster statistics in the bulk, it has not yet been applied to predict the distributions of clusters near charged interfaces in the electrical double layer (EDL) [80].
Here we develop a theory of aggregation and gelation of ILs in the EDL, and investigate the qualitative behaviour in this system. Before proceeding to the mathematical formalism of this theory, we briefly discuss the principles of the chemical equilibrium which hold in a theory which accounts for thermoreversible associations in the bulk and the EDL of ILs.
This discussion will be followed by a short preview of the results, which shall foreshadow the remaining paper. and m anions associated). This bulk cluster equilibrium, as seen in Fig. 2, was shown to agree with the cluster distribution independently computed from molecular simulations for ILs [49]. To study the EDL, there must be an equilibrium between ions near the interface and the bulk.
If a "standard approach" is taken for the EDL-Bulk equilibrium, i.e. where the ions in the bulk are assumed to accumulate in the EDL based on a Poisson-Boltzmann or Poisson-Fermi distribution from the bulk concentrations [80][81][82][83], the ions are not permitted to reversibly associate in the EDL [56,[59][60][61]. Such approaches bring ions to the EDL based on the charge and volume of the ion, but once the species are in the EDL, no further change of the ionic associations can occur. Moreover, with such an approach, it is not clear how the gel phase, if it forms, should be treated, and assumptions about how it will respond to an electrostatic potential will be required. For example, the gel could be assumed to not respond to the potential/field, and it remains as a constant background (charge), similar to solid electrolytes [84]. The gel is not a rigid solid [64], however, so this approach would have short-comings. Alternatively, an approach similar to how solvent is treated in the dilute electrolyte limit could be taken, i.e. the gel can be taken as some constant background dielectric, which would be trivially pushed out of the EDL based on how ions come to the EDL [80,83].
In the bulk, the gel phase is in equilibrium with free ions, and changes to the physical variables (such as temperature) can cause the balance of the equilibrium to shift [64]. Therefore, to model the cluster equilibrium in the EDL, we enforce that the cluster distribution in the EDL has the same form as in the bulk, but is altered by the differing volume fractions of cations and anions in the EDL. An equilibrium between the bulk and EDL then needs to be established, which keeps both the bulk and EDL cluster equilibria consistent. Note only one of the four EDL-bulk equilibrium (bulk free ions-EDL free ions; bulk free ions-EDL clusters; bulk clusters-EDL free ions; bulk clusters-EDL clusters) shown in Fig. 2  decaying to the right. The same notation for the clusters as Fig. 1 is used, but where each ion can form a maximum of 3 associations. equilibria shown in Fig. 2 being held consistently.

III. PREVIEW
The results of this theory shall be briefly summarised, as there are conceptual advances here which are beneficial to bear in mind before the technical parts are introduced. Firstly, the main advancement this paper makes is the account of the chemical equilibrium between clusters within the EDL. Consideration of this equilibrium permits a natural way to investigate strongly aggregating and gelating electrolytes in the EDL.
This theory predicts the following. For electrolytes which do not contain a gel, the large clusters are broken down / expelled from the EDL in favour of free ions, as schematically shown in Fig. 3. The account of clusters beyond free ions results in a differential capacitance response that is larger than that of just a free ion approach. Overall, this is the expected result for the pre-gel regime, and mainly differs to previous free ion theories ( For even larger potentials, these large aggregates are also broken down, and more and more free ions accumulate in the EDL, until the crowding regime is reached. Here crowding refers to an electrolyte composition where practically only free counter-ions exist, which can be achieved locally in the EDL. Gradually, the excluded volume effect causes the thickness of the EDL to increase. Thus the differential capacitance of the EDL, which is inversely proportional to the EDL thickness, starts decreasing at large potential drops across the EDL.
At small voltages, however, the differential capacitance is increasing [81,82], and therefore, the capacitance curve will be a typical "double hump camel" shape overall [81]. Again we find that accounting for the gel/cluster response increases and smooths out the response in comparison to a free ion theory.
Note that one must not expect the theory to provide the picture with 'molecular resolution', explicitly showing overscreening oscillations of charge density and the spatial structure within clusters. Within the large clusters or the gel, alternation of counterions and coions should give rise to such oscillations. The presented theory is able to provide a 'coarse-grained' structure of charge distribution in the EDL, spread over different clusters, the balance between which shifts with the polarization of the electrode. Needless to say, such theory will be thermodynamic, not covering any dynamic aspects of the EDL charging. Physically, these short-range associations between cations and anions in ILs represent the strong electrostatic correlations (beyond mean-field) between ions of opposite sign. Levy et al. [34] showed that the discrete charges of ions with a given packing (e.g. from molecular dynamics simulations) could be well reproduced by a nearest neighbour spin-glass theory that maximizes their alternating-charge ordering, which implies that our short-ranged model of Coulomb correlations should be accurate for these concentrated electrolytes. Moreover, in Ref. 49 it was found that IL ions have well defined 'hot-spots', where cations (anions) prefer to reside around anions (cations), which also supported these assumptions and permitted the functionality to be determined (from the number of hot-spots).
ILs are the simplest possible super-concentrated electrolyte, and should serve to clearly investigate the EDL of electrolytes which can form associations beyond ion pairs. For further details about the IL limit, see Ref. 49.
The electrolyte is treated on the level of a lattice-gas [49,64]. The number of lattice sites, Ω, is given by where N lm is the number of clusters of rank lm, and N gel +/− is the number of cations/anions in the gel phase. The cations and anions which are not part of the gel, i.e., all of the free cations, free anions and clusters of rank lm, are part of the sol. This is schematically shown in Fig. 4. In the absence of gel, all cations and anions are in the sol. Throughout, free ions and clusters shall not be explicitly referred to as being in the sol, as they are such by definition. Therefore, we do not explicitly include sol superscripts for free ions and clusters for clarity of notation. Equivalently, dividing by the total number of lattice sites yields where c lm = N lm /Ω is the dimensionless concentration (# per lattice site) of rank lm clusters, and the dimensionless concentration of cations/anions in the gel is c gel +/− = N gel +/− /Ω, i.e., the volume fraction c gel +/− = φ gel +/− . The volume fraction of a cluster of rank lm is given by φ lm = (l + m)c lm . The volume fraction of cations and anions in the sol phase is The total volume fraction of cations/anions is φ +/− = φ sol +/− + φ gel +/− . When there is no gel, the superscript sol shall be dropped for clarity of notation. In Fig. 4, the partitioning of the sol and gel is schematically shown.
The free energy of the cluster equilibrium [49,64] is taken to be where β = 1/k B T is inverse thermal energy, ∆ lm is the free energy of formation of a cluster of rank lm from free cations and anions, and ∆ gel +/− is the free energy change of cations /anions associating with the gel.
Establishing chemical equilibrium [49,64,66], as shown in the Supplementary Material (SM), the cluster distribution is expressed as where λ = exp(−β∆f +− ) is the ionic association constant, with ∆f +− = ∆u +− − T ∆s +− denoting the free energy of an association, determined by the binding energy (∆u +− ) and (configurational) entropy of an association (∆s +− ). Here W lm is the number of ways to arrange l cations and m anions in a Cayley tree [49,64], i.e. the combinatorial contribution to the free energy of formation of a cluster, which is given by The cluster distribution, c lm , is expressed in terms of φ 10 and φ 01 , but, in principle, these are unknown quantities. In the bulk, the volume fraction of cations and anions is, however, known. Introducing association probabilities, p ij , the probability that an association site of species i is bound to species j, permits the volume fraction of free cations to be written as The association probabilities can be determined through conservation of associations and mass action law to give the probability of cations binding to anions and anions binding to cations Note that the volume fractions have been explicitly retained in these equations, as these probabilities can be investigated as a function of the volume fraction of cations and anions.
When the probability of a cation (anion) having another cation (anion) connected through an anion (cation), a percolating ionic network can form. This condition is given by 1 = where the stars are used to denote the critical probabilities for formation of the percolating ionic network. When the probabilities reach this condition, the volume fractions of cations and anions in the gel and sol must be determined. Flory's treatment of the post-gel regime is employed, in which the volume fraction of free ions can be written equivalently in terms of overall association probabilities, p ij , and association probabilities taking into account only the species residing in the sol, p sol Using Eqs.

Bulk Composition
In the bulk, the critical value for the association constant for the formation of the gel is Refs. 49, 64 for more details) which drive the system towards a mixture rather than a pure state (of free ions). This has also been found by Lee et al. [85] when considering ion pair formation in ILs, where 2/3 of ions were found to be free. This was because the Debye screening length was much shorter than the radius of an ion, which meant the entropy of mixing decided the fraction of free ions, i.e. the IL was effectively treated as an ideal solution. Alternatively, this can be understood from the mass action law, Eq. (9), which shows that the association probabilities only vanish when λ = 0. The association constant is, effectively, the equilibrium constant of the formation of an association, which means that when λ = 1 it is equally favourable to form an association as not to, and therefore, entropy maximisation dictates that a mixture shall form. As the entropy of an association was found to be positive in Ref. 49, the association constant tends to 0 at large temperatures, which is the expected limit of no associations. Note that associations between cations-anions still occur, despite ∆f +− > 0 due to the negative binding energy.
Having described the extent of associations for different λ, we turn to understanding the cluster distribution in more detail. In Fig. 5 the volume fraction of each cluster of rank lm is shown as a function of the size of each cluster (l + m) for the bulk φ + = φ − = 1/2 with f = 3. This is plotted for various overall charges of clusters, with charges up to ±3 being considered. As it is the symmetric case, the volume fraction of a +q charged cluster is equal to the −q charged clusters. Therefore, these are summed together, and each different symbol represents a |l − m| = q, for q = 0, 1, 2, 3.
probability of anions binding to cations (p −+ ) increases. The opposite is true for a decreasing volume fraction of cations. This is a consequence of the conservation of associations, which For λ = 10, the association probabilities are sufficiently large to create a percolating ionic network in the bulk, as seen by p +− p +− being larger than the dashed horizontal line in Fig. 6(right). Again, p +− decreases with increasing φ + , and p −+ increases with increasing φ + . There comes a point, for both small and large φ + , that the volume fraction of cations (or anions) cannot sustain a percolating ionic network, and p +− p +− drops underneath the threshold.
Next, how these association probabilities influence the cluster distribution of the electrolyte is shown in Fig. 7 as a function of the volume fraction of cations. For λ = 1, the volume fraction of free cations (φ 10 ) increases with the volume fraction of cations (φ + ), and the volume fraction of free anions (φ 01 ) decreases with increasing φ + . At φ + = 1, φ 10 = 1 and φ 01 = 0, and vice versa for φ + = 0. This is a reflection of the changing association probabilities and volume fractions, as just discussed. The volume fraction of ion pairs (φ 11 ) and aggregates beyond ion pairs (φ lm>11 ) decreases for increasing and decreasing volume fractions of cations (relative to its bulk value). This is because of the average association For λ = 10, the volume fraction of free cations and free anions behave in an analogous manner to λ = 1, and the volume fraction of ion pairs is always extremely small. The aggregates beyond ion pairs behave qualitatively differently, however, with a highly nonmonotonic dependence on φ + . To understand this, the behaviour of the gel must first be described. In the bulk, the system is significantly gelled. As the volume fraction of cations increases (φ + ), the volume fraction of cations in the gel (φ gel + ) also increases and the volume fraction of anions in the gel (φ gel − ) decreases. This is because the prevalence of cations shifts the free ion-gel equilibrium to accommodate more cations. Therefore, the gel becomes charged and decreases in total volume fraction, until it reaches the critical volume fraction where it can no longer be sustained. This decreasing volume fraction of gel causes a dramatic increase in the volume fraction of aggregates beyond ion pairs, as closer to the gel point larger and larger clusters dissociate from the gel. When the critical volume fraction to sustain a gel is reached, there is no longer gel to create large clusters, and the decreasing average association probability causes a dramatic decrease in the volume fraction of these larger clusters.
Overall, the phenomenology described in this thought experiment is exactly what is expected for the behaviour of the clusters in the EDL: a transition from gelation/aggregation to crowding. Therefore, what is required is a connection between the volume fractions of cations/anions and the electrostatic potential, which can be linked to the spatial distribution of charges through the Poisson equation.
Inspecting the equations for p ±∓ and φ 10/01 demonstrates that, for a given volume fraction of cations, the volume fractions of free anions and free cations are not independent, but constrained by the associations. Therefore, if EDL-bulk equilibrium was established with "standard approaches", such as Boltzmann or Fermi functions for free cations and free anions, the system of equations (the above equations for the concentrations of each species and the Poisson equation) becomes over-determined and the cluster distribution is no longer consistent. Therefore, an alternative approach to investigate the EDL properties in a consistent way is sought. In the SM, this is shown explicitly from the free energy. Moreover, alternative approaches to make the cluster distribution consistent everywhere is outlined.
It is interesting to note that this observation also applies to approaches, such as Refs. 60, 61, where the associations were only treated through the explicit treatment of the remaining free ions. In the bulk, the number of free ions was determined by the binding free energy of cations and anions. The bulk concentration of free cations and free anions was then utilised in a Poisson-Fermi theory for the EDL [56,58,59,61]. This essentially treated the associations as irreversible in the EDL, but reversible in the bulk, which is also how Ref. 51 investigated strongly associating ILs within a sophisticated classical density functional theory. For thermoreversible associations in the EDL, the number of associations based on the concentration of cations and anions must also be determined in the EDL. This deficiency in such approaches [56,58,59,61] does not, however, qualitatively change the predictions of those theories in terms of the free ions.

C. Boltzmann Closure of Free Ions
To close the system of equations without over-determining, a single relationship is required to connect the electrostatic potential with the volume fractions. To achieve this, the following closure relationship is taken based on the ratio of free cations and free anions in where the bar has been used to indicate quantities in the EDL, and u is the electrostatic potential in units of thermal volts (which is given by 1/βe, with e denoting the elementary charge, i.e. the magnitude of an ion studied here). The factor α is a parameter which describes additional correlations beyond mean-field electrostatics (defined such that 0 ≤ α ≤ 1), as introduced in Ref. 59, and shown to work well in Ref. 86. Since λ accounts for short-ranged attractive interactions between cations-anions, only the short-ranged repulsion between cations-cations/anions-anions are considered in α. This prevents double counting the short-ranged attraction between cations-anions. Moreover, α is not assumed to be based on the free ion fractions, but the total ion fractions, since the gel can also become charged.
This parameter acts to smooth-out the response of ions to the electrostatic potential, as Poisson-Boltzmann (or Fermi) is well known to overestimate the response [80,83]. This is analogous to dressed ion theory, where the charge of an ion is rescaled to smaller values because of correlations between ions [38].
Equation (14) assumes that free ions behave in a Boltzmann way, whilst being consistent with the cluster distribution and incompressibility constraint. This expression can be used to solve for the local volume fraction of cations or anions within the EDL, which then determines the local free cation and free anion volume fractions. These can be used in the to find the concentrations of all clusters. Moreover, when the free cation/anion volume fractions are known, so are the association probabilities. This means that the gel can also be treated in a consistent way.
For λ 1, the association probabilities tend to zero p +− = p −+ ≈ 0, which means Using the incomprehensibility condition, this can be solved to giveφ + − φ − = − tanh αu for the charge density, which is the expected result when there are no associations [81,82].

Linear Response
When there are associations, the approach cannot generally be solved analytically, but at linear response some insight can be gained. Taking a linear expansion of Eq. (14) and introducing a symmetric perturbation of free ions,φ 10 = φ 10 + δφ f andφ 01 = φ 01 − δφ f , yieldsc 10 = c 10 (1 − αu) for the free cation concentration andc 01 = c 01 (1 + αu) for the free anion concentration. Using this in the cluster distribution yields Introducing the EDL volume fraction and probabilities as the bulk perturbed by a small value,φ ± = 1/2 ± δφ andp ±∓ = p ± δp, respectively, the volume fraction of cations changes where p = (1 + f λ − √ 1 + 2f λ)/f λ is the association probability in the bulk.
Therefore, the Poisson equation takes the form where 0 and are the permittivity of free space, and relative permittivity and v is the volume of an ion (i.e. a lattice site). The screening length is based on the ionic strength of the cluster distribution, with κ = v 0 /e 2 β. When there are no associations, i.e. p = 0, the Debye length is recovered. In the opposite limit the association probability tends to 1, and extremely large screening lengths can emerge. This expression for the ionic strength in terms of p was previously derived in Ref. 49 from modifying the expression for the weight average degree of aggregation. Therefore, the Boltzmann closure of free ions appears to be well justified at linear response in the pre-gel regime.
Gebbie et al. [36,37] have suggested that ILs screen the field in a Boltzmann way with few free ions. Thus, this closure relation could also be assumed to hold in the post-gel regime.
which is an analogous expression for the sol quantities. Therefore, the change in the volume fraction of the gel is then given by These results can be used to decompose the total screening length into the contributions from the sol and the gel.
In Fig. 8  only the gel contributes to the screening at linear response. This indicates that approaches based on free ions [60,61], which follow the sol contribution closely [49], underestimate the ability of the electrolyte to screen.

Non-Linear Response
Generally, the differential equation which needs to be solved is given by A constant charge boundary condition at the interface is taken and a zero electrostatic potential in the bulk is used. Numerical results for the solution to this non-linear differential equation shall be shown in the next section. Note a Stern layer is not considered here for clarity of discussing the new results. In the SM, a step-by-step guide of how to implement the equations to obtain a numerical solution is given.

A. EDL Structure
To start, the description of the non-linear solution to Eq. (22) in terms of the electrostatic potential and charge density shall be outlined. The results for these are shown in Fig. 9 as a function of distance from an interface which has a charge of 0.02 Cm −2 for the studied association constants. The electrostatic potential (in units of thermal volts) and charge density (in units of the charge per unit lattice site) are both found to monotonically decay from the interface, as expected from a local density approximation.
For λ = 1, the potential and charge density decay to the bulk values within ∼ 30/κ.
The crowding regime, defined by c 01/10 ≈ 1 and c 10/01 ≈ 0, which has to be taken here because an incompressible IL can only accumulate charge density not number density in the EDL, can clearly be seen, since the dimensionless charge density reaches −1 at the interface. For λ = 10, the electrostatic potential and charge density decay further from the interface, owing to the larger screening lengths from more associations. The crowding

B. Differential Capacitance
The differential capacitance, C, can be calculated numerically within this theory through whereσ is the dimensionless surface charge density, C 0 = 0 κ is the Debye capacitance without associations, and u 0 is the dimensionless potential drop across the entire EDL.
When λ 1, i.e. all free ions, and α = 1 the differential capacitance numerically calculated within the presented theory exactly reproduces that of Refs. 81, 82.
In Fig. 11 the numerical differential capacitance as a function of potential is plotted (solid line) alongside the analytical expression derived in Ref. 59 where γ is the free ion fraction based on a similar free ion fraction (dotted line) to the cluster distribution; where the additional factor of √ αγ comes from the definition of κ taken here.
Note that Ref. 60 is based on just the free ions in an IL, referred to as a free ion theory (FIT), where there are only thermoreversible associations in the bulk.
For λ = 1, the free ion fraction is approximately 0.16, which is used in Eq. (24) of the FIT. Since the free ion fraction is significantly smaller than 1/3, a clear "camel" [81] shaped differential capacitance curve is obtained from Eq. (24), i.e. initially the differential capacitance increases before a maximum is reached, after which the differential capacitance is governed by universal charge conservation laws [81]. In contrast, the numerical solution for the theory presented here only has a slight "camel" shape, with the differential capacitance at zero charge also being larger than the FIT.
For λ = 10, the free ion fraction is found to be approximately 0.01. Again, the FIT has a significant U-shape near the potential of zero charge, which goes through a maximum before the crowding regime is reached and the differential capacitance decreases again (note that the differential capacitance follows C/C 0 ∝ 1/ √ u 0 for u 0 > 40, which is where |ρ| = 1 reaches in Fig. 9, and therefore, the employed definition of the crowding regime is reasonable). The new theory presented here, in contrast, only has a slight "camel"-shape which is stretched out over the potential range, and again has a larger differential capacitance at zero charge.
Overall, there are two features of the theory presented here against that of Ref. 60: (1) - The screening length obtained here is always smaller than Eq. (24) because larger clusters are explicitly accounted for. Therefore, the capacitance at zero charge is always larger. (2) -The break-down of clusters and gel causes the crowding regime to be reached at smaller potentials than that of Eq. (24), which makes the "camel" shape less pronounced.
The two association constants investigated here can be considered to be at two different temperatures for the same IL, since they are related to the free energy of an association through λ = exp{−β∆f +− }. Therefore, larger temperatures correspond to smaller λ. As found in Ref. 60, we also observe a transition from a camel to bell shape with increasing temperature, owing to the dissociation of ions. Moreover, the capacitance at zero charge increases with decreasing associations because there are more charge carriers to screen the potential.

Comparison to Experiments
In Ref. 49 In Fig. 12 the differential capacitance curve from the new theory and the free ion theory (FIT), i.e., Eq. (24), is shown. The new theory can reproduce the differential capacitance curve reasonably well with λ = 0.6, as seen in Fig. 12. For λ = 0.6, the free ion fraction is 0.12, which is close to that of Ref. 88. For the FIT with a free ion fraction of γ = 0.12, the differential capacitance curve cannot reproduce experimental data well, being too "camel" shaped. The FIT with the fitted value of γ = 0.27 from Jitvisate and Seddon [87] is also shown for comparison. Therefore, the new theory appears to be able to rationalise the results which were not in quantitative agreement with the theory before.
There are still a couple of free parameters, other than the free ion fraction which was independently taken from Ref. 88, evaluated to obtain the differential capacitance curve in Fig. 12. One of these parameters is α = 0.07, which was the value fitted by Jitvisate and Seddon [87]. Using this parameter, and setting the dielectric constant to the bulk value (for [Emim][TFSI] this is = 12 [87]) obtains a capacitance at zero charge which is too large. Therefore, to obtain a similar capacitance at zero charge, the dielectric constant is set to = 1. As such value is unphysical (at least the polarisability of electronic degrees of freedom ions would contribute the value of 2), such fitting of demonstrates a short-coming of the simple, local theory presented here, and motivates further improvement. Note, to permit a transparent comparison with Jitvisate and Seddon [87], a Stern layer has not been included, which would have reduced the capacitance at zero charge.

A. Underscreening
The presented theory has implications for the interpretation of the experiments in Refs. [35][36][37][38][39][40][41][42]. Inverting the expression for the screening length in terms of the association probability, we obtain In experiments [35], the screening length multiplied by the inverse Debye length is κ ≈ 100 for ILs. Using this value yields an association probability of 0.9997 for a functionality of 3.
This association probability produces ∼ 3 × 10 −11 free ions, with the rest being gel, which corresponds to a λ in excess of 1 × 10 6 and ∆f +− ≈ −14/β. At large distances from the charged interface, we have found that it is actually the gel which screens the long-tail in electrode charge, not the free ions. This could suggest that ILs are not dilute electrolytes, but a gel which can screen electrode charge. In fact, Jurado et al. [41,42] found 1-Hexyl-3-methyl-imidazolium ethylsulfate forms lamellar-like structures on charged mica surfaces, which are conceptually similar to the gel, albeit with more order.
Despite this, the capacitance at zero charge is (in this theory) still fixed by the screening length , which means the underscreening paradox remains [45]. Further development of the presented theory, such as through taking into account the spatial structure of clusters/gel or accounting for electrostatics beyond that investigated here [30], could gain further insight into these puzzling experiments.

B. Overscreening
Ma et al. [51] performed sophisticated classical density functional calculations of a coarsegrained IL with varying degrees of ion pairing, up to that suggested by Gebbie et al. [36]. It was noted that the charged density was somewhat insensitive to the extent of ion pairing, which suggested a link between ion pairing and overscreening. This was further shown by Avni et al. [58], where a link between clusters (ion pairs, triplets and quadruplets) and the linearised overscreening theory of Bazant-Storey-Kornyshev (BSK) [28] occurs in the long-wavelength limit. The BSK theory is known to underestimate overscreening in ILs, however, and the interpretation of it can be subtle [89][90][91], with some suggesting that short-range cation-anion repulsion occurs. The alternating "co-polymer" structure of the gel/large clusters is conceptually similar to extended overscreening structures (bulk or in the EDL). If one accounts for the spatial structure of larger clusters (than ion pairs, triplets and quadruplets [58]) or the gel, extended overscreening could potentially be obtained.
Here, for large association constants, we found that the EDL structure has two regimes: one where electrode charge is screened by free ions, and another where it is screened by the gel. The free ions accumulate near the interface, with the gel persisting towards the bulk.
Conceptually, this could be considered similar to the overscreening-crowding transition [28].
The reader should also bear in mind that the EDL theory presented here is a simple, local density approximation, which cannot explicitly account for overscreening or the internal charge distribution of clusters. However, while a spatial map between ions in clusters, assumed to be Cayley trees here, and charge density does not yet exist, Cayley trees have well defined shells around an ion in a cluster. One can loosely interpret the number of shells around an ion with the distance from that ion, which can be utilised to further establish the link between associations and overscreening.
Let us consider a cation at a 'central' lattice site and examine the probability of an anion being at any of the lattice sites n shells away from the central cation, similar to the pair correlation function between cations and anions. The anion can be present on that lattice site either by being linked to the central cation as part of the same cluster, p, or by random, p r − . There is also a possibility that a cation could be present on that site by random, p r + . Assuming that unassociated species are uncorrelated, p r − = p r + = p r , and taking into account the incompressibility condition, we obtain p r = (1 − p)/2. Thus, the total probability of observing an anion one shell away is (1 + p)/2. For the next shell, the probability that an anion is there via a link to the central cation is exactly zero, as we require alternating associations. However, there is some probability that an anion (or cation) is there by chance.
In this case, the random probability of an anion is (1 − p 2 )/2, which is in fact the total probability of seeing an anion two shells away. Similarly, the analysis for three shells away yields an anion probability of (1 + p 3 )/2. Thus, we can generalise the probability of seeing an anion n shells away from a central cation to In Fig. 13, P +− is plotted as a function of n, where we see that the correlation probabilities display pronounced overscreening when λ = 10, but only modest overscreening for λ = 1.
Therefore, pronounced overscreening could be more consistent with the gel than clusters.
While the developed theory does not spatially resolve the clusters/gel, further development for the spatial distribution of ions in clusters/gel could shed light on this connection.
In this spirit, one issue of the theory (which presumably needs to be dealt with before this can be achieved) is how can extremely large clusters screen the field if it decays on a significantly shorter length scale than the clusters? This is a well known issue for theories of ILs with simple, local density approximations when no associations are included. Taking into account associations increases the screening length relative to the Debye length, but the clusters can become extremely large in size. Therefore, the issue can remain for finite clusters, and perhaps a finite cut-off needs to be introduced for clusters which can contribute to the screening of electrode charge, as is the case of the conductivity of ions in ILs [49,79].
For the gel, it is conceptually reasonable to think of its composition as locally changing, just as an iceberg melts from its surface. For clusters, their individual structures cannot locally vary as they are discrete. In any case, the theory should work for a large number of free ions and a small number of free ions, when there is significant gel, and should act as an interpolation scheme between these two regimes, with the largest error occurring only right at the gel-point where there are significant numbers of very large clusters.
In addition to the above issue, is the applicability of a local density approximation when there are substantial variations of the charge density, on a smaller length scale than the screening length, i.e. the fact there is internal structure to the clusters. The presented theory cannot capture these short-scale variations, with it only being able to deal with mean-field volume fractions of clusters of different ranks and the gel. By developing theories which more accurately account for spatial correlations from electrostatics and non-electrostatic interactions, and which also account for the internal structures of clusters, the overscreening structures and screening lengths observed in surface-force experiments could, perhaps, be more accurately described.

C. Assumptions of the Theory
The developed theory has some shortcomings, as just described, but how well do we expect it to handle the averaged volume fractions of each cluster? There could potentially be a number of issues with the developed theory at a charged interface which could cause it to break-down. For example: 1. Do Cayley tree clusters exist near an interface?
2. Do specific interactions with the interface dominate the compact layer?
3. Do ion pairs and aggregates rotate in the electrostatic fields?
4. Does the association free energy depend on the electrostatic field?
For each of these questions, we shall state what is expected to happen, and how to confirm/falsify these expectations. Point 1 is the major question to be answered, with 2 and 3 being minor if 1 is found to hold reasonably. In the description of 1, some details of how to test the cluster distribution shall also be given. Overall, it is expected that the cluster distribution in the EDL does not exactly hold, but we believe it is a sufficiently good approximation to understand the qualitative behaviour.

Cayley trees and cluster distribution in the EDL
We have shown that clusters are destroyed in favour of free ions at large electrostatic potentials, and therefore, the assumption of Cayley tree clusters is not important at large fields. What needs to be tested, however, is the low-voltage regime of overscreening. Close to the interface, the presence of the interface could block association sites of the ions and drive the ions to form intra-cluster loops. The assumption of Cayley trees applying in the EDL and the assumed cluster distribution can only be confirmed from molecular simulations.
In bulk ILs, McEldrew et al. [49] found that the spatial distribution function (SDF) of cations (anions) around anions (cations) have well defined "hot spots" (for where the concentration is over twice the bulk value), where the ions of opposite sign prefer to reside. When the cations and anions reside in each others hot spots, an association was defined [49]. This highly directional interaction between cations and anions is a consequence of the complicated shapes of the ions in ILs, and is fundamental to the formation of Cayley tree clusters [77], which ILs have been found to form in the bulk [49]. Spherical cut-offs to define associations or kinetic criteria [62] are not sufficient to test this approximation. Moreover, Lennard-Jones obtain the volume fraction of each cluster, the clusters will have to be assigned to bins. The size of these bins should be larger than the average size of clusters, which might become problematic right near the gel point when the cluster sizes become very large.
The gel will also be evident from using the association criteria of McEldrew et al. [49] as a percolating ionic network throughout the simulation. The theory predicts that the gel continuously exists, but changes composition as a function from the interface. Inspection of the percolating ionic network should be able to reveal if there is an equivalence between overscreening and gel screening. Moreover, testing the connection between the associations and overscreening in the bulk would also be interesting.

Specific interactions with the interface
The ions can interact with the interface 'specifically' through non-electrostatic interactions, which often causes one type of ion to accumulate at the interface without an applied voltage, owing to the specific interaction of one type of ion being more favourable than the other. These specific interactions with the surface are not accounted for in the presented theory, and will further contribute to breaking the assumed cluster distribution. Therefore, the presented theory might only work well for the diffuse part of the double layer, not the compact layer of ions which are in contact with the interface.
In a similar spirit to the presented theory, several theories to describe the voltage dependence of the compact layer with clusters were developed by Damaskin-Frumkin-Parsons [92,93]. In those theories, the equilibrium between the free states and clustered states were selfconsistently determined in an electrostatic potential, which gave rise to a voltage-dependence of the differential capacitance of the compact layer. To account for specific interactions with the interface, one approach could be to have an additional equilibrium between the compact layer, described by the theories of Damaskin-Frumkin-Parsons [92,93], and the diffuse part of the EDL.
Again, molecular dynamics simulations can be utilised to test if specific interactions cause a break-down of the cluster distribution, using the cluster criteria of McEldrew et al. [49].

Orientation of clusters
As overscreening is a spatially ordered structure, it is expected that clusters do orient in an electrostatic field. However, we do not expect them to behave as fluctuating Langevin dipoles. It has been shown in Refs. 49, 53 that the lifetime of associations is of the order of ∼1-10 ps, which is presumably too short for the rotation of an ion pair or cluster. Therefore, the orientation of clusters must be a weaker contribution than fluctuating Langevin dipoles.
In Ref. 94, ion pairs (which were not permitted to dissociate in the EDL) were treated as Langevin dipoles, where an additional bump in concentrations of ion pairs and a higher propensity for a "bell"-shaped differential capacitance curve because of dielectric saturation was found. While ion pairs might not behave as fluctuating Langevin dipoles, the large clusters and gel could have vibrational modes which contribute to the dielectric screening of the electrolyte. This could be introduced in a similar way to Ref. 94, as some effective dielectric constant which is saturated to lower values in the crowding regime. Thus, it is expected that this issue is not substantial for the averaged volume fractions, but could cause a breakdown of the assumed cluster distribution, as outlined in the SM. The issue of orientation of clusters can again be confirmed by molecular simulations, provided the cluster distribution can be calculated as described in the previous point.

Field dependent binding
Related to the previous points, if an ion pair is oriented along the electrostatic field decay, the ion pair will presumably be stretched. This could be captured through an electrostatic field dependent binding energy, and therefore, association constant. We would expect the binding energy to decrease with increasing electrostatic fields. This will cause the associations to break-down more than shown here. The presented theory could be considered to overestimate the extent of clusters in the EDL, then. Again, this issue can be confirmed from computing the cluster distribution in the EDL and comparing it against the presented theory.

VII. CONCLUSION
In summary, we have developed a theory for the EDL in ILs which accounts for thermoreversible associations, based on that of McEldrew et al. in bulk ILs. The developed theory is constructed from assuming the bulk cluster distribution applies in the EDL, but with modulated volume fractions of cations and anions, which is coupled to the Poisson equation through a Boltzmann closure relation of the free ions. This theory was shown to recover the expected linear response behaviour, and the large potential limits. Therefore, it should serve as a qualitative description for the EDL of an IL with thermoreversible associations. For strongly associating ILs, the free ions crowd near the interface, but far from the interface the gel screens the electrode charge. The differential capacitance was found to have a larger capacitance at zero charge than free ion theories, from the fact that the clusters can also screen, and the transition to "camel" shape occurs at smaller free ion fractions than models based on free ions.
ILs provided an extremely useful test ground to develop the theory presented here, since