Ising models of charge storage in multifile metallic nanopores

Ising type models of charging of conductive nanopores with ions have already been proposed and investigated for single file cylindrical or single layer slit nanopores. In such pores, the state of ions, the coulombic interactions of which are exponentially screened by their images in pore walls, was named superionic. In the present work we extend the analysis of the superionic state to nanopores that can accommodate multiple rows of ions. By grouping multiple charges in the same row into ‘supercharges’, we map the arrangement of ions in polarised electrodes on a multi-row Ising model in an external field. We investigate one-, two- and three-row cases, which we solve exactly, using a purpose-built semi-numerical transfer matrix method. For pores of different radii, which can accommodate the corresponding number of ion rows, we calculate the dependence of the electrical capacitance and stored energy density on electrode potential. As in charging the single file pores, we find that in narrower pores higher energy densities can be achieved at low applied potentials, while wider pores perform better as the voltage is increased.


Introduction
Since the dawn of electricity, it has been important to find a way to efficiently store electrical energy. This has led to the development of multiple technologies used for energy storage. One such category of devices is supercapacitors (SCs). SCs are high-capacitance energy storage devices with higher charging/discharging rates than those of batteries. Electrochemical double layer capacitors (EDLCs) are one type of SC being used today. First patented by Becker in 1957 [1], they store potential energy at the electrode/electrolyte interface via charge separation between the electrical double layers (EDLs) in the electrolyte at the cathode and the anode [2]. Due to the nanoscale separation between the charge of each electrode and the centre of mass of the counter charge in the electrolyte, EDLCs are able to achieve much higher capacitances and hence much higher stored energy densities than classical dielectric capacitors [3]. In addition, since the process used to charge EDLCs is non-Faradaic, they have far longer lifetimes than rechargeable batteries, and so can typically be cycled upwards of a million times [4].
The capacitances of the double layers at the cathode and anode, combined in series, are each inversely proportional to the distance from the electrode to the centre-of-mass of the counter charge. Since this distance is of nanoscale, the capacitance per unit surface area reaches tens of microfarads per cm 2 . However, the related increase of stored energy densities comprises only a tiny fraction of the energy density of advanced batteries [5,6]. As the capacitance and energy density are roughly proportional the overall area of the electrode/ electrolyte interface, one way to increase the energy storage in EDLCs is to increase that area. Nanostructured electrodes with volume filling surfaces are therefore widely used to boost the volumetric (and gravimetric) energy density.
In order to achieve the highest surface-to-volume ratio, one uses nanoporous electrodes. As long as the pores in the electrode are still wide enough to endorse a full EDL at the pore walls, the capacitance will indeed scale proportionally to the interface area [7]. Hence such porous electrodes will provide the same normalised capacitance (capacitance per electrode surface area) as a flat electrode, but a higher overall capacitance for a device of the same size. However, the simple proportionality of capacitance to the interface area no longer holds true when the pore sizes become of the order of the size of ions. That interesting nanoscale effect was first observed in [8,9], where it was found that the normalised capacitance of carbide-derived-carbon nanoporous electrodes increases as the average pore width decreases and approaches the size of the ion. This result boosted further research of supercapacitors with ultrananoporous/nanostructured electrodes.
These experimental findings motivated many groups to research and justify the observed 'anomalous' normalised capacitance rise for nanoscale pores. A number of theoretical models were constructed to investigate and predict capacitance changes with respect to factors such as pore width, ion diameter and pore ionophilicity [10][11][12][13][14][15]. One such explanation attributed the effect to the formation of a, so called, 'superionic state', in which interionic interactions between neighbouring ions are screened by electrons in the pore walls, allowing more ions of the same charge to pack into the pore for a given voltage. This results in a higher capacitance and energy density [11,12]. Using this explanation, it was made possible to theoretically replicate a number of experimentally observed trends in equilibria as well as charging dynamics [11,12,16,17].
In the course of these investigations, a simple model of a 'single-file nanopore' was proposed [14]. The model considered a metallic, cylindrical pore of infinite length containing a single row of ions. The approach of [14] may look similar to treating ions as hard spheres with point charges at their centres (as per the restrictive primitive model, shown to be an acceptable approximation when comparing with more realistic ionic structures [18]), but in fact it is more 'phenomenological' and thereby more general.
The model relies on several key assumptions: firstly, that the pore is ionophilic and hence always packed with ions. Secondly, that the electrostatic interactions between ions are screened by image forces, and so decay exponentially in the lateral direction [14] so that only nearest neighbour ions interact with each other. Hence the system can be mapped onto an exactly solvable one row, two-state Ising model, in which the charges of the ions and the applied potential are analogous to the spins and magnetic field of the Ising model, respectively [19]. Thus, the system energy can be modelled via the thermally-scaled Hamiltonian of the 1D Ising model with nearest neighbour interactions in the presence of an external field (applied potential): Here k B T is the thermal energy, the charge of ions i and i + 1 is denoted by s i and s i+1 , respectively, u is the dimensionless electrostatic potential inside the pore, scaled to thermal voltage k B T/e, where e is the elementary charge, and α is the coupling constant between nearest neighbour ions at closest approach, given by [14]: where U(d) is the interaction potential between ions of the same sign and d is the diameter of the ions in the pore, which determines the distance of closest approach between them. For an ideal metal: with a, the pore radius, and L B , the Bjerrum length given by L B = e 2 / k B T (in Gaussian units, used throughout this paper) where is the effective dielectric constant.
For this system, exact solutions of the averaged value of 'spin' (in this case, charge) and its derivative with respect to voltage, the response function χ(u), can be found [20]. This gives the normalised differential capacitance: This model shows that (at least at small deviations from the potential of zero charge of the electrode), for a single file row of ions, maximum normalised capacitance indeed occurs for a pore of the same width as the ion diameter [14]. Furthermore, as the pore widens, the normalised capacitance decreases due to lower screening of the interactions between ions. This result is to be expected; as the interactions between neighbouring ions become stronger, the energetic barrier which needs to be overcome to pack many counter-ions into close proximity increases and hence the capacitance falls. However, while single file pores achieve the highest normalised capacitances at small voltages, as first shown experimentally [8,9], they have shortcomings in other respects that need addressing. One of them relates to their charging dynamics, which directly affect the power density of an EDLC containing such pores. As simple steric arguments would suggest, ultra-narrow pores should lead to slower charging than wider pores. Several studies looking at the charging dynamics of nanopores of varying widths confirmed this [17,21]. This problem is usually mitigated by implementing a hierarchical structure of pores containing both narrow nanopores (for achieving maximum capacitance) and wider pores (for achieving better ion transport and accessibility) [22].
As well as this, one also needs to take into account how the charge storage depends on the applied voltage. A study in [16] found that, while at lower voltages a narrower pore was able to store more energy, as the voltage was increased wider pores performed better in this respect. But that conclusion was proved for single file cylindrical and single ion layer slit pores. Thus, it would be interesting to test this important conclusion for wider pores, containing multiple ion rows.
When considering multiple rows of ions, a number of new factors need to be accounted for. Firstly, extending the Ising model approach adopted for the single file [14] onto the multi-row case, the transfer matrix used in finding the solutions for average charge, response function and normalised capacitance [15,23] becomes much larger. Hence, the previously used method of finding the eigenvalues of this matrix algebraically becomes unfeasible, and a new approach is required. Secondly, whereas in the 1D case only interactions between nearest neighbours along one axis need to be considered, now we need to account for interactions both along and between the rows of ions. Next, whereas for the 1D model it is roughly sufficient to know the interactions between ions along the central axis of the pore, for the multi-row cases we have to account for the electrostatic interactions of the ions positioned eccentrically with respect to the axis, which results in a more complex expression for the coupling constant. Finally, the ions can now be arranged in many distinct ways with respect to each other. For example, for a model with two rows of ions, the ions could be placed in two parallel lines in line with each other, or shifted vertically with respect to each other, or twisted with respect to each other. These countless ways of arranging the ions make the analysis of such models more challenging.
In this investigation, the previous Ising model approach taken in [14] in analysing a single file system is extended to the case of multiple ion rows. This is done by presenting a novel semi-numerical method for dealing with larger transfer matrices and hence evaluating the capacitance and charge storage of wider cylindrical pores. A number of static configurations are also considered to begin comparing stability of different arrangements of ions for a given number of rows. Figure 1 shows a schematic of possible two and three row configurations, as well as a one row model for comparison. As mentioned in the Introduction, unlike in the single file model, nearest neighbour ion interactions now need to be considered both along and between the rows. Therefore, the Hamiltonians for such systems become far more complicated than that of the single file system presented in equation (1). This is made easier to deal with by assigning to each group of inter-row nearest neighbours a single 'supercharge', s i , which contains as many charges as there are rows.

Hamiltonians and transfer matrices for multi-row models
Using supercharge notations, and in units of k B T, the Hamiltonian reads differently for the two and three row cases. Unlike the single file Hamiltonian presented in equation (1), which only considers nearest neighbour interactions along the single row, the multi-row Hamiltonians also consider interactions between ions in neighbouring rows (see figure 1).
For the two row case, s i = p i q i , where p i , q i ∈ {±1} and: where α 1−4 correspond to coupling constants between neighbouring ions; u is the dimensionless electrostatic potential of the pore as in the model of [14]. The set of possible supercharge combinations for this system is: These combinations will be used when building the transfer matrix T for the two row model. The three row model is constructed via the same method, but now each supercharge consists of three charges, such that The set of possible supercharge combinations for this system is: (9) For both models, the transfer matrix is constructed using the supercharge combinations given by equations (7) and (9) for the two and three row models, respectively. Each element of the transfer matrix T is given by: where indices m and n indicate the element in the SC row matrix. They run from 1 to 4 for the two row model and 1 to 8 for the three row model. Thus, the transfer matrix is a 4 × 4 and 8 × 8 matrix for the two and three row cases, respectively. The details of the application of transfer matrix theory are described in section 4.

The quantities to calculate: average charge, capacitance, stored energy
To find the important properties (average charge, response function, normalised capacitance and energy density) of the systems, transfer matrix theory is used. In the limit of an infinite lattice (as for the cases presented here) only the largest eigenvalue of the transfer matrix, λ 1 , and its derivatives with respect to u are needed. See section 4.1 for the explanation as to why. The average charge per supercharge is given by: Therefore, the ionic charge per unit surface area of the pore, σ, is found as in [15]: where p is the pore circumference. Since the response function, χ(u), is the derivative of S i with respect to u, it can be expressed as: By combining equations (4) and (13), an expression for the normalised capacitance in terms of λ 1 and its derivatives can be found: Finally, the energy stored by the pore after charging the electrode from u = 0 to u = u, in units of k B T, is given by [16]: which can be evaluated such that E is expressed in terms of λ 1 and its derivatives: Equation (16) is derived in section 4.3. The novel seminumerical approach for calculating λ 1 and its derivatives for large transfer matrices is explained in section 4.

Coupling constants
As mentioned above, when it comes to calculating the coupling constants between ions in the multi-row models, a different equation from the one used in the single file case, equation (2), is needed. Therefore, the energy of the electrostatic interaction between ions in eccentric position with respect to the axis of the cylindrical pore is given by the following expression. Using cylindrical coordinates, the potential between two point charges (of the same sign) placed at (0,ρ 0 ,0) and (z,ρ,φ) in a metallic pore of radius a is given by [24]: (17) where z is the lateral distance between the charges, ρ and ρ 0 are the radial distances of the charges from the central axis, ϕ is the azimuthal angle between the charges, l m = 1 for m = 0 and l m = 2 for m > 0, J m is the cylindrical Bessel function of the first kind of order m, and k nm is the nth positive root of J m (x). As in equation (3), e is the elementary charge and is the effective permittivity of the pore interior due to polarisability of the ions. Equation (17) simplifies to equation (3) when considering ions on the pore axis and when |z| a knm . The coupling constant α is related to the interionic potential as per equation (2): Equations (17) and (18) allow us to describe interactions of charges of ions in any arrangement with respect to each other and the axis of the pore. For a system with a given number of rows, the same Hamiltonian can be used, altering only the coupling constants α to reflect the specific arrangement of ions.

Approximations
As in the model of [14], the pore is assumed to be ionophilic, so that the ions in the pore are closely packed. The intra-ion charge distributions are modelled with point charges at their centres. The ions are assumed to be of the same size, so that the ion centres in each row occupy the sites of a 1D lattice with a lattice constant equal to their diameter (for a single file model the extension of the case when cations and anions have different sizes was investigated in [23]). The exponential screening of electrostatic interactions of ions in the lateral-direction which facilitates the dense packing of ions of the same charge is still present for the wider row models, as confirmed by the e −|z| term in equation (17). This allows us to account only for the nearest neighbour interactions in our Hamiltonians. But note that equation (17) considers the electrode to be ideally metallic, whereas most modern EDLCs are made up of carbon materials. Therefore it is appropriate to question whether the exponential screening is still applicable in such a case. A study performed in [13], which took into account field penetration into the electrode pore walls, suggested that the screening was not as strong as in ideal metals, but was still exponential. Furthermore, a study in [25] used density functional theory to investigate the nature of the screening for different ions inside gold and carbon nanotubes. Once again, the results showed that, while the screening in gold was stronger, both the gold and carbon nanotubes showed sufficiently strong screening for narrow tubes. The screening was considered to be strong enough to result in the superionic state for pore widths of a similar magnitude to ion diameter. This justifies the use of approximation of exponential screening, on the condition that the pore-walls are not atomically thin.

Results and discussion
For all the plots presented, the ion diameter d was set to 0.7 nm. The temperature T was taken as room temperature, 298.15 K, and the dielectric permittivity as 2 (due to the electronic polarisability of ions for a densely packed pore). For these conditions, the Bjerrum length L B is equal to 28 nm.

Comparison of charge and capacitance for one, two and three row models
From the single row model investigation of [14], it follows that before a second row of ions can fit, maximum capacitance in the low voltage regime is achieved when the pore diameter is the same size as that of the ions and it decreases for wider pores. In the present section the one, two and three row two-state models are investigated in the broad range of electrode polarizations. These cases are specifically analysed for pore diameters of 0.70, 1.40 and 1.51 nm, respectively. The ionic charge per unit surface area and differential capacitance against the voltage drop between the electrode and the bulk of the electrolyte for these cases are shown in figure 2.
As shown by the σ − V plot, the multi-row models require a significantly higher applied voltage to begin replacing coions with counter-ions than the single row model. This is intuitively clear and can be explained by two reasons; firstly, the pore for the one row case is much narrower than for the multi-row cases. This leads to higher screening of interionic interactions. Secondly, in the multi-row models each ion has more ion neighbours to interact with. Combining these two factors, the ions in the multi-row models have stronger and more numerous interactions with other ions than the one row model. This manifests itself in two ways: firstly, a more stable arrangement is adopted by the ions at zero applied voltage. Therefore, the system needs to overcome a higher energy barrier to leave this initial stable state and a higher applied voltage is needed to push it to do so. This explains why the σ value in figure 2(a) remains at zero for larger voltages for the multi-row cases than for the one row model when the voltage is increased. Secondly, even once the arrangement of anions/ cations begins changing in the multi-row models, there will be intermediate stable arrangements before the final state consisting of only counter-ions is reached. This explains why plateaus are observed in the σ − V plots. The width of each plateau indicates how stable that arrangement is; a more stable configuration will require a higher energy step to shift from it and hence the voltage will need to be increased by a larger increment, i.e.: the plateaus will be wider for more stable configurations. From the σ value of a plateau, the stable configurations can also be inferred, as will be shown in section 3.2.
The differential capacitance for these models is shown in figure 2(b). As expected, the number of peaks for each model corresponds to the number of transitions between stable arrangements in the ionic charge per unit surface area case. Hence, the two and three row models have two and three capacitance peaks, respectively. The height of these peaks corresponds to the gradient of σ(V) between the plateaus. Hence, as per the stability explanation offered above, it is unsurprising to see the maximum capacitance decrease as the pore size increases. This is because the higher row systems have more stable states due to the stronger and more numerous interactions, and so change from those states less readily than lower stability states as the voltage is increased.

Comparison of stability for two row configurations
The two row two-state model was tested for three distinct cases, with a different ion configuration in the pore for each case. These cases are sketched in figure 4(a). It is important to note that, although the ions are configured differently in each case, the ion density is the same for all three cases. Each case corresponds to a 1.4 nm wide pore, so exactly two ions can fit side-by-side. The coupling constants α 1−4 for each case are calculated subject to equation (18). Up to two decimal places, they are given in table 1 below. The highest value of α, 9.97, corresponds to the interaction of ions with no lateral separation between them. Clearly, the interionic interaction screening by the pore is weakest for this case.
The three cases were compared by looking at the ionic charge per unit surface area and differential capacitance for each configuration over a potential range. These quantities were calculated using equations (12) and (14), respectively. Corresponding plots can be seen in figure 3.
A lot of useful information can be gained from the plots in figure 3. Looking initially at the σ − V plot, it is clear that each of the three configurations respond to the applied voltage differently. Initially, when the electrode is not polarised, all of the configurations maintain a charge-neutral state. However, when the applied voltage gets high enough, the driving force for counter-ions (for this electrode polarisation: anions) to fill the pore becomes sufficiently strong. Hence coions become replaced by counter-ions, resulting in the average charge inside the pore decreasing. Configuration 1 requires the highest applied voltage for this to happen, suggesting that it has the most stable charge-neutral state of the three.
Surprisingly, however, once the initial state is disrupted, the rate of counter-ions replacing co-ions continues to vary. This manifests itself in intermediate plateaus in the σ-V plots, as seen above. Importantly, the value of ( pd/e)σ at these intermediate arrangements gives an indication as to how these arrangements might look. Configurations 1, 2 and 3 have plateaus occurring at ( pd/e)σ values of −1, − 2 3 and −1, respectively. Hence for configurations 1 and 3, one in every four ions is a cation, whereas for configuration 2, one in every three ions is a cation. Thus, we decipher the arrangements of the ions in these stable regions, as shown in figure 4(b).
The transitions between the regions described above lead to the capacitance peaks shown in figure 3(b). Since differential capacitance, given by equation (14), depends on the derivative of ionic charge density, the height of the capacitance peaks depends on the gradient of σ with respect to voltage. This explains why the capacitance plot has multiple peaks for each configuration, as they correspond to the transitions in the σ − V plot.   The values of the coupling constants play a critical role in the shapes of all three plots; the stronger the interaction between the ions, the higher the applied voltage will be needed to overcome the attraction between cations and anions and the unfavourable ion-ion repulsions when the pore is filled mainly with counter-ions. This is the basis of the superionic state first proposed in [11]. As the pore widens, the screening of electrostatic potential between ions also decreases, resulting in the ions interacting more strongly together and so requiring higher applied voltages to achieve a pore containing only ions of the same charge. The strength of the coupling constants also determines the positions of the plateaus.
The analysis of two row two-state stability is only performed here for three configurations, in which the ions are held in rigid positions with respect to each other. This is a highly idealised model, as one would expect the position of the ions to be flexible, shifting to adopt whichever arrangement had the lowest energy for a given voltage. Adding more configurations and testing their energies, one could reach a plot resembling how a real two row system might act.

Energy density
The energy stored by the pores of different sizes, per their unit surface area, as a function of applied voltage, was also investigated. The configurations for the one, two and three row models were the same as in section 3.1. The results are plotted in figure 5. They show a very similar trend to that found in [16]; the narrower, one row system stores a higher energy at lower voltages while the multi-row models store significantly higher energies when the applied voltage is increased. This is not an unexpected result, and can be justified by the same reasoning as in section 3.1; the multi-row models have more stable initial states, due to the stronger and more numerous interionic interactions present. Therefore a higher voltage needs to be applied to begin filling the pore with predominantly counter-ions. Hence the multi-row pores can only store more energy when a sufficiently large potential is applied (see: the 'pressing-a-spring' concept [26]).
As well as this, the plateaus previously seen in ionic charge per unit surface area are once again observed. This is also to be expected, since the intermediate stable states will be disrupted at higher potentials. In this respect, the correspondence between figures 5 and 2(a) can be clearly seen, with the plateaus occurring in the same voltage ranges.
These results show that, while narrow pores can achieve optimal energy density at low voltages, wider pores perform significantly better when the applied voltage is increased. Clearly, the type of pore used should depend on the specific application. However, it does open the avenue of multi-row pores being used for capacitors with higher applied voltages 2 .

Model limitations
The model presented has a number of limitations. The first of these is that the ions are fixed in rigid positions with respect to each other. In reality, the ions would be able to shift slightly to take on more stable arrangements. This would result in less well defined plateaus in the σ − V plots as the different arrangements would each have different corresponding stable configurations. Secondly, the arrangement of ions is kept the same for each row model until a new row can be added. Depending on the system, the more stable arrangements for wider pores could result in the ions of like-charge trying to reduce repulsions by shifting to the sides of the pore. For example, this would explain a smoother transition between a single row and two row model as the pore is widened.
This model also does not consider a number of other factors which could affect the results; firstly, it does not take into account differences in ion sizes, which would lead to different packing for the cation rich and anion rich cases [23]. It also does not consider van der Waals forces between ions (which would differ for cation/cation, cation/anion and anion/anion neighbours), imperfect screening by electrons in the pore walls [13] or interactions of ions across the pore walls [27], as well as likely other phenomena.
Furthermore, the cases studied thus far here all refer to pores containing a small number of ion rows, with all of the ions being adsorbed to the electrode surface. It would be interesting to increase the width such that non-adsorbed ions could appear, to see whether a bulk phase would form and how this would be reflected in the normalised capacitance. Unfortunately, as the number of rows increases, the size of the transfer matrix for the system grows exponentially, and hence higher computational power would be required. In such cases, direct Monte-Carlo will be a much better way forward.
Finally, experimental testing of the observed trends on well-structured monodispersed samples would be needed to confirm the results found, particularly to approve the existence of the intermediate stable ion arrangements. The geometrical model that we studied here is very idealised, and does not consider many factors, such as pore size dispersion and non-uniform pore width.
However, as a model for densely packed, ionophilic cylindrical structures of straight, narrow, multi-row pores, that could be fabricated in the near future, it presents a good 'starting point'. It be developed over time to go beyond the mentioned constraints.

Transfer matrix theory
The transfer matrix theory is most easily introduced for the simple one row, two-state Ising model. This model is defined by several parameters: The charge at any site i can be ±1 s := (s 1 , s 2 , ..., s N ) This is the configuration of charges in the system Ω = {±1, ±1, ... ±1} N This is the space of charge configurations, of length N The Hamiltonian for this system is given by: where α and u are the coupling constant and dimensionless electrostatic potential, respectively, as given in equation (1). The probability distribution of charge configurations on Ω is given by: where Z N is the partition function and β = 1/k B T . In this analysis, energy is taken in units of k B T, which is formally equivalent to taking β = 1. By definition, the probability of all the individual charge configurations will sum to 1, allowing the partition function to be defined as: For such systems, the average charge per site, S i , is found as: Taking the partition function as a function of the external field, u, the following derivation can be performed: Therefore, the average charge per site is given by: Hence the situation is reduced to the evaluation of the partition function. Since in this model the length of the pore is taken to be much larger than the ion size, analysing the system in the limit N → ∞ gives a good approximation. Equation (19) can be re-written: where h(s k , s k+1 ) = αs k s k+1 + us k .
For this Hamiltonian, the partition function Z N can be evaluated: where T is the transfer matrix defined as: (1,1) .
Computing T N would require significant calculation. However, only its trace needs to be found. Since T is a square matrix, via the eigen decomposition it can be diagonalised: where U is a matrix composed of the eigenvectors of T, U −1 is its inverse and D is a diagonal matrix defined as: where λ 1 and λ 2 are the eigenvalues of T such that |λ 1 | > |λ 2 |.Therefore, the trace of T N can be found by: Hence: Therefore, since |λ 1 | > |λ 2 |, in the limit N → ∞ only the largest eigenvalue of T contributes to the value of the partition function. So, the partition function can be expressed as: By combining equations (23) and (32), an expression for the average charge per site in terms of the largest eigenvalue of the transfer matrix can be obtained: which is the same expression as equation (11) and also allows the derivation of equations (13) and (14) to find the response function and normalised capacitance, respectively. From the analysis performed thus far, it is clear that several key quantities need to be found to calculate the average charge per site, response function and differential capacitance: λ 1 , ∂λ 1 /∂u and ∂ 2 λ 1 /∂u 2 . For the one row model, it is relatively simple to obtain expressions for these terms algebraically. λ 1 is found by evaluating the expression: where T is the transfer matrix defined in equation (27). Solving equation (34) gives a quadratic equation for λ: Solving this expression for the maximum value of λ yields: Hence an expression for the largest eigenvalue of the transfer matrix has been obtained and can be used to find the desired properties of the system. The situation with multi-row systems is, however, more complicated.

Semi-numerical approach for large transfer matrices
The method presented above, in which equation (34) is used to find the largest eigenvalue of the transfer matrix, becomes increasingly difficult to use as the size of the transfer matrix grows. For example, it would be unfeasible to find algebraic expressions for the largest eigenvalues of the transfer matrices presented for the two and three row systems in section 2.1. Should an even higher number of rows be used, the matrices become even larger. Hence a new method is presented, in which the values are instead computed semi-numerically. Via numerical computation, it is trivial to find the largest eigenvalue of a given transfer matrix. However, finding its derivatives is not so simple, and requires the use of Jacobi's formula and its second derivative analogue: where the adjugate adj() of a square matrix is the transpose of its cofactor matrix [28]. For the derivations of these expressions, see the appendix. Now consider a matrix B(t) = A(t) − λ(t) · I, where λ(t) is an eigenvalue of A(t) such that: Therefore: Combining equations (37) and (40) allows the derivation of the expression: This expression allows the calculation of ∂λ 1 /∂u and hence S i . Analogous to this, the solution for the second derivative of λ can be found by combining equations (38) and (41): where λ and A are functions of t. This expression in turn allows the calculation of ∂ 2 λ 1 /∂u 2 and hence χ(u) and normalised differential capacitance, C(u). This semi-numerical method can be used to find the differential capacitance for multi-row Ising models with large transfer matrices. Its strength is that it can be so easily applied to different models, with only the Hamiltonian and corresponding transfer matrix needing to be altered.

Energy calculation
The energy for a system charged to the potential value u = u can also be found in terms of λ 1 and its first derivative with respect to u. This is done by evaluating the following integral [16]: Remembering equation (23), the response function, χ(u), can be written in the form: Finally, the capacitance can now be written as C(u) = −C * χ(u), where C * is a constant defined as C * = LB 2πad . Therefore the integral in equation (44) can now be evaluated via integration by parts: Then, by using equation (32), the final expression for the energy is obtained: Once again the semi-numerical method presented in section 4.2 can be used to find the first derivative of λ 1 at u = u. As noted throughout the paper, Gaussian units are used and energy is given in units of k B T.

Conclusion
A semi-numerical transfer matrix method has been presented to investigate the properties of EDLC nanopores containing multiple rows of ions. This method proposed the idea of 'supercharges' to model multiple charges in a single layer. It can be extended to both multi-row and multi-state models easily.
The behaviour of multi-row models was studied for both two and three row two-state (cations and anions only, no voids or solvent molecules) systems. The results are significantly different from those for single row pores. This can be attributed to the decreased inter-ion interaction screening for wider pores and the presence of more ions interacting with each other. The results reveal a possibility of stable intermediate configurations of anions/cations as the applied voltage is increased. These stable configurations give rise to multiple energy barriers that need to be overcome to fill the pore with counter-ions. The existence of these intermediate states was indicated via plateaus in the charge accumulation, σ, value when increasing applied voltage. The stability and arrangement of the states could be inferred via the width and σ value of the plateaus, respectively.
Once again, single-file pores were found to achieve maximum normalised capacitance at low electrode polarizations. This was justified by the easier packing of like-charge ions into the pore compared to multi-row models. However, the energy density investigation showed that whether a single-or multi-row pore option was preferred depended on the applied voltage, with single-file pores storing more energy at low voltages and multi-row pores becoming preferred as the voltage was increased. This result indicates that, whilst single row pores result in the highest normalised capacitances, the consideration of operating voltage for energy density means that slightly wider pores may be the preferred option for larger voltages 2 .
In summary, a novel approach of analysing multi-row models allowed the investigation of densely packed ionophilic pores. The results obtained indicate that, as the pore width increases, it becomes increasingly difficult to reach the Kondrat and Kornyshev 'superionic' state [11]. The results also open the door to follow-up investigations such as exploring the multi-row models for ionophobic pores. This could be done by adapting the model to consider holes as well as cations/anions, as in [15]. Also, as mentioned in section 3.4, this model holds the ions in fixed positions within the pore. In reality, the ions would likely rearrange to the lowest energy arrangement for a given voltage. Therefore, Monte-Carlo simulations could be used to model systems in which the ions change their packing arrangement in response to voltage changes. Such simulations would also test the qualitative predictions of the present study. Experimental tests would require samples with ideal monodisperse cylindrical pores, but modern nanotechnology could provide such architecture in the near future.
where I is an identity matrix of the same size as B. Now a new term D(t) is defined: det(B(t)). (A.1) Hence, by using the property for square matrices that det(A · B) = det(A) · det(B): A new variable Y is now defined: to (A.4), the following expression is obtained: Hence: Performing the same treatment for the δt 2 terms on either side gives: (A.7) Equations (A.6) and (A.7) are Jacobi's formula and its second derivative analogue respectively. However, these forms of the equations are not suitable due to containing B −1 terms; for the purpose of finding the derivative of λ, B is defined as B = A − λ · I, with det(A − λ · I) = 0, and it is therefore non-invertible. Therefore the expression