Entropy is often described as a measure of disorder. Its importance is that the free energy of a system involves a trade-off between the entropy and the energy. The simplest example of this is the Helmholtz free energy $F=E-TS$, which applies to a system with a fixed number of particles and a fixed volume at constant temperature. Because of the complexity of biological systems, there are many instances in which the $TS$ term of the free energy either dominates the energy term, or at least plays an important role in the free energy minimization that determines equilibrium and drives the dynamics of the system.

As a consequence, the determination of the optimal structure of many biologically important systems requires minimization of the free energy, and hence this competition between the internal energy and entropy. Often in biological systems the changes in entropy are greater than the changes in internal energy. For instance, in the hydrophobic effect, an unfolded protein lowers the entropy by ordering the water molecules, and so the protein prefers to be in the folded state[1,2]. In the chloroplast stroma, it has been shown that there is an entropy-driven attraction that determines chloroplast ultrastructure through spontaneous Mg${}^{2+}$-induced stacking of membranes[3].

In sickle hemoglobin, the aggregation of monomers into polymers[4] is also entropy driven, with the internal energy and entropy in a delicate balance. In fact, Cao and Ferrone [5] measured the chemical potential ${\mu }_{PC}$, which is the "glue" that holds the polymer together, and found that the entropic contribution to this chemical potential is around $-11$ kcal/mol, which is more than 100% of the chemical potential itself. Another example is F-actin, whose fibrils stack in cross-linked rafts when positive alkaline earth ions are in the solution[6]. In F-actin counter-ions form one-dimensional charge density waves that have a periodicity equal to twice the actin monomer spacing, coupling to twist distortions of the oppositely-charged actin filaments[7]. These phenomena and others owe their existence to Coulomb interactions between the constituent parts, and the geometry of the underlying structure can also play a vital role.

Thermodynamically, the distribution of charge in the solution is governed by a competition between energy and entropy. The bound state in which the ions are condensed on the surface of the macroion/polyelectrolyte is energetically-favored, and the continuum state, in which the ions are free to drift in any direction, is entropically-favored. The balance between these two factors in minimizing the free energy has been shown to vary significantly based on the geometry of the macroion [8]. This interplay between electrostatics and geometry is particularly important for polyelectrolytes in solution, and they have shown to produce a wide array of intriguing phenomena in many different systems. One example is DNA, which in the presence of a physiological salt solution (a 0.1 molar solution of NaCl) is usually negatively charged, with a double helix DNA molecule strongly repelling another. However, if DNA is in a dilute salt solution in which a positively charged polyelectrolyte, such as spermine or sperimidine, has been added to the solution, it has been shown to roll up into a tightly-packed torus[8,9,10,11]. In fact, DNA is usually in a very compact state in cells and viruses.

Historically, the most common description of charged solutions has been Poisson-Boltzmann theory, but continued research has indicated that charged species in solution can have far more complex and counterintuitive effects than simple charge screening [8]. In 1969 Manning proposed [12,13] that a portion of the ambient ions condense onto (i.e., attach to) the surface of a charged macro-ion, partially neutralizing the bare charge. This occurs up to the point at which the energy required to condense another ion equals the available thermal energy $k_{B}T$ at temperature T, where $k_B$ is the Boltzmann constant. The proposed effect, later termed “Manning condensation," marked a significant departure from linearized Poisson-Boltzmann theory that predicted only exponential screening. In Manning’s treatment of polyelectrolytes, which are long chains of charged subunits, the ion condensation was addressed [8] separately from the ions that remain in solution, which were treated using linearized Poisson-Boltzmann theory [12,13].

Further refinements in the treatment of these ambient-ion solutions have been motivated in part by surprising effects observed in DNA that cannot be accounted for using Poisson-Boltzmann theory. By using multivalent cations, modifying the salts in the solutions, and even using alcohol solvents, it is possible to cause the DNA to undergo a radical structural transition into a variety of new geometries, including rod-like bundles and toroids [14]. The toroidal structures, in particular, have received considerable attention in the biological community because such toroidal packing motifs are employed by spermidine and other molecules to contain their own DNA in small volumes [15]. A variety of techniques, including cryo-transmission electron microscopy [15] and UV spectroscopy [16], have enabled direct observation of the formation of these toroidal condensates, and similar studies on other biological polyelectrolytes like F-actin [6] indicate that condensation is a general phenomenon. If the electrical interactions between the polyelectrolyte, such as DNA, and the ambient-ion solution are purely those predicted by the Poisson-Boltzmann description, then two like-charged polyelectrolytes always repel one another, and such condensates will not form. For these condensates to be stable, the net electrical force between like-charged polyelectrolytes must be $\mathit{attractive}$, which is incompatible with the Poisson-Boltzmann theory of ions in solution.

Many of the approaches to understanding the role of polyelectrolytes in biological systems have adopted models of continuum electrostatics, and some of these have focused on solving the Poisson-Boltzmann equation for a cylindrical model of DNA in the presence of multivalent ions[8,17]. However, it has been shown that including charge correlations is essential to understanding these systems[18,19]. Recognizing the need for including fluctuations beyond the mean field, Ha and Liu used a field theoretic approach that produces a systematic expansion taking these effects into account[20,21,22,23]. They employed a simplified model of a DNA particle as a charged cylindrical rod composed of cylindrical charged segments and considered these rodlike particles in the presence of polyvalent counterions. Their first model consisted of two rodlike DNA molecules[20], and they found attraction between the rods. When Podgornik and Parsegian suggested that the fluctuation forces between such rodlike polylectrolytes might not be pairwise additive[24], Ha and Liu extended their calculation to include bundles of rodlike particles and included a one-dimensional form factor depending on the ionic size in order to incorporate short-ranged correlations along the rod length approximately[21,22]. They found that the breakdown of pairwise additivity can lead to qualitatively new behavior, although they still found that that the rods attracted one another under some conditions. Shklovskii then pointed out that the previous results seemed to be independent of radius, and that, in a highly correlated system of negatively-charged cylindrical rods in the presence of positive counterion charges, those positive charges can form a Wigner crystal, which will be three-dimensional and closely packed for small rod radii and form two-dimensional Wigner crystals on the surfaces of rods for large radii[25]. Ha and Liu saw merit in the arguments for the large radius limit but disagreed with the the two-dimensional crystal lattice limit, which they said was due to Shklovskii incorrectly assuming pairwise additive interactions between various surfaces, which would make that limit invalid[22].

To show how the radius of a DNA chain could affect multivalent-counterion-induced attraction between negatively-charged DNA chains, Ha and Liu used a modified model of DNA[23]. They assumed, as before, that a DNA rod is a stack of short cylinders that have a finite radius and length, but to treat two-dimensional charge fluctuations on the surface of the rod, they divided the counterions into two classes, free and condensed, with the condensed counterions modifying the local charge at the surface of the cylinder, giving rise to charge fluctuations on the two-dimensional surface of each cylinder. With this model, they found that the competition between attractive and repulsive interactions tended to balance one another, resulting in no attraction at all for thick rods. They also found that that for valence of the polyvalent cations greater that around 3, the spacing of chains in a bundle and the size of bundles appears nearly independent of the nature of the bundling agent. This is because the increased valence leads to a stronger screening as well as a stronger attraction.

Recognizing that geometry was playing a crucial role in the attraction, Nguyen and Shklovskii [26,27] proposed a simple theory of charge inversion that considers the structure of the polyelectrolyte together with its electrostatic interaction with the substrate. The idea was that there is fractionalization of charge in which a polyelectrolyte can either neutralize the charge or reverse it depending on how it attaches. In that model, a single double-helix strand of DNA is represented by a rigid cylinder with two one-dimensional lattices of negative charges $-e$ in helices around the surface to represent DNA’s double helix of negatively-charged phosphates. In order to model the polyelectrolyte solution in which the DNA is immersed, they take the positively-charged species to be freely-jointed chains. These adsorbing species have multiple charges that may partially attach to the surface, with excess charge protruding into the solution, as shown in Figure 1. They assumed that all the negative charges on the DNA are neutralized by the polyelectrolytes, so that there will be no vacancies, and because there is excess charge dangling into the solution, the charge on the surface can be not only neutralized but reversed. They also assume that the only role of the background salt solution is to provide screening for the interaction between the polyelectrolyte chains. Since they assumed no vacancies on the one-dimensional chains, they could easily count the possible configurations. In calculating the electrical potential, they took the negative charges of the DNA to be spread uniformly on the surface of the cylindrical DNA surface and used simple electrostatic arguments with screened potentials to calculate the energies involved. For polyelectrolytes in a 0.1M NaCl solution with DNA, Nguyen and Shklovskii[26] give in their Eq. (6) an estimate of the net charge inversion, the charge density on DNA divided by its bare charge density. The largest value of that estimate is about 0.07.

Since the simple calculation of Nguyen and Shklovskii[26] did not include the discrete nature of the phosphate sites in calculating the electrostatics and used a simple version of the entropy, Bishop and McMullen[28] decided to use a more rigorous formalism to calculate the charge reversal in DNA. Starting with the simplifying assumption that the polyelectrolytes could be considered as dimers (two units in length) on DNA, Bishop and McMullen allowed for vacancies and used a model of the discrete location of phosphates on the DNA surface, with interactions between charges located on the various sites. In that model, doubly-charged polyelectrolyte molecules attach to the DNA either parallel or perpendicular to the DNA surface and confirmed that under the right conditions, the charge on the surface could be inverted. The use of the lattice model allowed them to accomdodate vacancies, which is not possible if the charge is treated in a continuum manner. Inspired by the work of Ha and Liu[20,21,22,23], they used a field-theoretic approach that employs a loop expansion, where the mean field contains polyelectolyte ions adsorbed on the surface of the DNA, similar to their approach. While the general formalism for the one-loop expansion was presented in that paper, the numerical calculations stopped at the mean field level, and the current work extends that to include fluctuations at the one-loop method. The emphasis here is on the entropy and the role it plays in the configuration of charged species attached to the lattice. The details of this model will be explained in detail in the next section.

The results of these theories indicate that the correlation effects in the solution are also strongly-dependent on the system geometry. To understand this, consider the electric potential outside spherical, cylindrical, and planar positively-charged surfaces in vacuum. For spheres, the potential decays as the inverse of the distance; for the cylinder, the potential decreases logarithmically with the distance; and for the plane, the potential decreases linearly with the distance. Gelbart et al. [8] assert that, in an ambient-ion solution, the entropy of the point ions varies logarithmically with their concentration, and therefore with the distance from the cylinder. Thus they conclude that, in a crude comparison, for spherical geometry, the $r^{-1}$ potential is dominated by the logarithmic entropy; for planar geometry, the entropy is dominated by the linear potential; and for cylindrical geometry, both the energy and the entropy vary logarithmically.[29] In this paper, we will incorporate the impact of the ions in solution on the free energy through the chemical potential of these ions, as detailed in Appendix B.

The helical geometry of DNA, then, sits precisely balanced on the fulcrum between energy- and entropy-driven processes under physiological conditions. The geometry of the biomolecule plays an even more crucial role when the highly-charged ions in solution are not point charges but have geometries of their own. Such is the case for DNA immersed in a solution of polyelectrolytes. These polyelectrolytes can be proteins or other fragments of biomolecules which are routinely found in the nucleus [30], so that, again, the central biological processes occur in precisely the most difficult regimes to model.

For charge inversion on DNA with polyelectrolytes, “physiological conditions" require incorporating the combined effects of charge correlation, thermodynamic fluctuations, crowding, and geometrical considerations all at once. As we have discussed in this introduction, there has been considerable work in addressing all of these issues. Some approaches treat only the geometry, with no interactions[31]. Others include both geometry and interactions, but use a continuum model for electrostatics that neglects discreteness effects[32]. The lattice gas dimer model is unique in its simultaneous consideration of all these effects, and the thermodynamic and geometrical idealizations it makes can be systematically improved.

While the formalism of Bishop and McMullen[28] included both a mean field theory and corrections due to fluctuations or charge correlations, the computed results were obtained only at the mean field level. In this paper, this simple model is extended to calculate the fluctuation corrections and the entropy of the system. The purpose is to determine the importance of the fluctuation terms for inverting the charge and to see whether these terms have a significant impact on the entropy of the system. A preliminary version of this work was in Sievert’s master’s thesis[33]. In Section 2 of this paper, we outline the model of the charged binding sites on DNA and present the computation of the entropy due only to the hard-core repulsion that prevents multiple binding on the same site. In Section 3 we explain the geometry of the double helix and show how it can be represented as a one-dimensional lattice and in Section 4 derive the form of the potential and determine the orthogonal transformation that diagonalizes it. In Section 5, we use functional integral techniques to derive the partition function, which uses a Gaussian integral identity to perform the sum over configurations exactly, at the expense of an integral over a new auxiliary field. In Sec Section 6, we show how the grand canonical potential can be computed order by order, in which the first two terms are the mean field and the one-loop correction to the mean field. In Section 7, we examine the saddle point equation that defines the mean field, and this is used to calculate the entropy, the number densities of all species, and the charge density. In Section 8, we find the explicit form for the inverse-Hartree-field-fluctuation propagator. In Section 9, we include the one-loop order terms in the calculations to reveal the effects of fluctuations. Finally, in Section 10, we compare the results of the various approximations. Details concerning the Gaussian integral identity and the chemical potential of dimers in solution are given in Appendix A and Appendix B.

Our model system starts with a simplified version of the DNA molecule itself, as shown in Figure 2, with double helical chains of phosphate ions connected by base pairs. The phosphate ions are represented by red balls and blue balls, with the base pairs represented by yellow lines. The two different chains are labeled "up chain" and "down chain", which correspond to the direction of the carbon atoms in the sugar backbone. We will assume that each phosphate site has charge $-e$, where e is magnitude of the charge of an electron, and that there are no charges either between the phosphate sites along the helical chain or in the vicinity of the base pairs. In solution are polyelectrolyte ions of charge $+2e$, which we call dimers. These dimers are assumed to be the length of the spacing between two phosphate ions along a helical chain, which is $a^{\left(1\right)}=0.684$ nm$=6.84$ Å, which is longer than a diatomic molecule. Possible candidates for these species are the diamines, 1,3-diaminopropane (DAP${}^{2+}$), putrescine (Put${}^{2+}$), and cadaverine (Cad${}^{2+}$). An extension of this model could include more highly charged spermidine(Spd${}^{3+}$) or spermine (Spn${}^{4+}$). These polyamines are shown in Figure 3, with the lattice spacing $a^{\left(1\right)}$ shown as the dashed line for comparison. Although we will only consider the doubly charged species in our model, the model and analysis could be extended to these more highly-charged species in the future.

Deng and Bloomfield have shown using Raman difference spectra, for the systems they studied, that in the presence of spermidine or spermidine these polyvalent cations bind electrostatically near the DNA phosphates[34]. van Dam et al. agree with this conclusion and suggest that the A form of DNA is stabilized by polyamines fitting perfectly between the phosphate ions[35]. They conclude that DAP, which has charge $+2e$ probably has the best fit to the phosphate lattice of all the polyamines, and that longer species, like spermine and sperimidine are longer and probably also interact with base pairs. In our model, we are assuming that there is no interaction with the base pairs.

The choice of Bishop and McMullen[28] to model dimers was motivated by the wealth of studies (refs. [36] and [37], among others) in the literature about dimer models in other branches of physics, as well as for geometrical simplicity. We will continue to use this same model. Dimers, the shortest polyelectrolytes, have only two possible orientations when adsorbed on DNA, assuming that the length of the dimer is comparable to the helical spacing between charged sites on DNA, as shown in Figure 4. As we stated earlier, this should not be confused with a diatomic molecule, which would not be long enough to span the space distance between two phosphate ions on one of the helical chains of DNA.

In our model, a dimer lying on the cylindrical surface of the molecule must lie parallel to the helical strands, neutralizing the charge on two adjacent sites. Otherwise, the dimer must adsorb perpendicular to the cylindrical surface, extending one end radially out from the central axis of the cylinder and inverting the charge on a single site. Charge inversion by dimers, then, is quite similar to charge inversion due to two species of point ions: one monovalent species representing parallel-adsorbed dimers and one divalent species representing perpendicular-adsorbed dimers. In the previous work, Bishop and McMullen[28] modeled the adsorption of dimers in a lattice-gas model as a two-component solution of point ions, allowing the possibility of vacancies, and used field-theoretic methods to describe the thermodynamics of the system. They carried their calculations to a mean-field level of approximation of the inverted charge on DNA, but did not calculate the entropy. Their work confirmed the possibility of charge inversion within this model.

Those computations yielded the charge per site on the DNA helix as a function of the chemical potential, or equivalently of the concentration of the polyelectrolyte in solution. While the Nguyen-Shklovskii[26] calculation assumed complete filling of the lattice, the Bishop-McMullen approach allowed for vacancies, represented as negatively-charged sites, in addition to the neutral or positively-charged sites arising from dimer adsorption. However, the lattice-gas model, which assumes all sites are equivalent, does not take into account that the parallel dimer occupies two sites. In this model, the binding energy of the parallel dimer is taken to be ${\epsilon }^{(\Vert )}$ and for the perpendicular dimer ${\epsilon }^{(\perp )}$, and these energies can be adjusted to account for the difference in occupancy. In practice, the occupation of two sites is mimicked by making the binding energy of the parallel dimer twice as large as that of the perpendicular dimer, and we will use the same approach here in the numerical calculations. The assignment of binding energies to these species is a simple way of including the complicated electrostatics that allows the polyelectrolyte ions to bind to the surface of the DNA. This is analogous to the "fermion" model of Ha and Liu[22], which assumes that each site can either be empty or occupied by one counterion. Here, as in the earlier Bishop and McMullen work, we have three species, parallel dimers, perpendicular dimers, and vacancies.

This approach allows us to describe adsorption of either species independently for each site. Any configuration of the DNA-dimer complex in this model can then be described by identifying the type of dimer (parallel, perpendicular, or none) adsorbed on each site on the DNA molecule. Such a model also resembles the lattice gas model of condensed matter physics [38], which treats the ways of distributing particles of different types onto a regular lattice of sites. In this section, we assume that the three different species on the lattice do not interact with one another, and in this case, the structure of the lattice does not matter. We call this the noninteracting model. This does not actually mean that there are no interactions. The interaction of the polyelectrolytes with the DNA can actually be quite strong, depending on the values of ${\epsilon }^{(\Vert )}$ and ${\epsilon }^{(\perp )}$, and their attachment to the lattice is analogous to what Ha and Liu call condensed counterions[22]. Also, we will assume that there is only one species per site, which corresponds a hard-core repulsion between sites.

In developing the formalism, it is convenient to consider a vacancy as a third species $\gamma$ of particle, because then we can impose the constraint that each site is singly occupied, either by a parallel dimer $(\gamma =\Vert )$, a perpendicular dimer ($\gamma =\perp$), or a vacancy ($\gamma =v$). For each our three species of “particles" that can reside on our lattice sites, we define a relative charge $q_{\gamma }$ in units of the magnitude e of the charge of an electron. These relative charges are then $q_v=-1$ for vacancies, $q_{\Vert }=0$ for parallel dimers, and $q_{\perp }=+1$ for perpendicular dimers. We will assume that the binding energy ${\epsilon }_{\gamma }$ of species $\gamma$ to the lattice depends only on the type of species and not the location. We will be specifying each lattice point by its location ℓ along helix b, with the pair $(\ell ,b)$ specifying that lattice position. Then the quantity $n_{\ell ,b}^{\gamma }$ will be the number of particles of species $\gamma$ on the lattice site at ℓ on chain b. For each chain, ℓ extends from $-\mathcal{N}$ to $\mathcal{N}$, such that the total number of sites is ${\mathcal{N}}_{\mathrm{sites}}=2(2\mathcal{N}+1)$, and we employ periodic boundary conditions. There will also be free polyelectrolyte ions in solution, and their influence on the stability of the adsorbed dimers will be through their chemical potential.

We begin by considering the Hamiltonian ${\mathcal{H}}_{\mathrm{NI}}$ for this noninteracting system (that is, with no interactions between different sites aside from the hard-core repulsion that blocks double-occupancy), which is Because we have a system that exchanges particles with its surroundings, specifically the ions in the solution surrounding the DNA, which attach to the surface, we work in the grand canonical ensemble. If the system contains $N_{\gamma }$ particles of type $\gamma$ in equilibrium with its surroundings, with an average internal energy E, the grand canonical potential ${\Omega }_G$, which is a function of the temperature T, volume V, and chemical potentials ${\mu }_{\gamma }$ for each species $\gamma$, is written as[38,39], where the number of particles of type $\gamma$ is given by Then the entropy can be written in terms this potential as where $\beta =\frac{1}{k_{B}T}$ and $k_B$ is Boltzmann’s constant. In this partial derivative, the volume and all the chemical potentials ${\mu }_{\gamma }$ of all species $\gamma$ are held constant. It is convenient to define the grand canonical partition function $Z_G$ in terms of the grand canonical potential ${\Omega }_G$ as where the sum over configurations includes all the accessible states of the system. The grand canonical potential can alternatively be written in terms of the partition function as and this allows us to write an expression for the entropy in terms of $Z_G$ as:

For the noninteracting lattice, the average internal energy is represented here by the Hamiltonian (1), $E={\mathcal{H}}_{NI}$, and the noninteracting grand canonical partition function $Z_{\mathrm{NI}}$ becomes where we have suppressed the G subscript for simplicity.

Since a vacant site does not really correspond to a particle, we recognize that energy and chemical potential of the vacancy must be related such that ${\epsilon }^{\left(v\right)}-{\mu }_v=0$. In addition, because the dimers adsorbed on the surface and those in solution are in equilibrium, ${\mu }_{\Vert }={\mu }_{\perp }={\mu }_{\mathrm{dimer}}$ is the chemical potential of dimers in solution at the appropriate concentration. We thus need an estimate for ${\mu }_{\mathrm{dimer}}$. Approximating the dimer as a uniformly-charged cylinder, this value is shown in Appendix B to be $\beta {\mu }_{\mathrm{dimer}}\simeq 0.79$ at the physiological temperature of $T\simeq 310$ K.

The average occupancy ${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}$ for species $\gamma$ per site is found by taking the derivative of this partition function with respect to ${\mu }_{\gamma },$ This can be verified by taking the derivative of Eq. (8). which is by definition the average number per site of species $\gamma$.

Continuing with the evaluation of the partition function, we note that the exponential of a sum can be written as the product of exponentials, enabling us to rewrite the partition function (8) as In order to simplify and to appreciate the physical meaning of this expression, it is useful to define the relative activity [40] of species $\gamma$ in the noninteracting lattice-gas model, given by Note that this is independent of the lattice site ℓ or chain b. The partition function can then be written in the simple form Because we have assumed that there can be only one species per site, parallel dimer, perpendicular dimer, or vacancy, the sum over configurations can now be done over each site separately, where there are three possible configurations This is the same as saying that $n_{\ell ,b}^{\left(\gamma \right)}=1$ for one and only one of $\gamma =\perp ,\Vert$, or v and is zero otherwise. This means that and so the grand partition function is Since every term in the product is the same, the expression in parentheses is simply raised to the power ${\mathcal{N}}_{\mathrm{sites}}$, giving

At this point, we can easily see that we can obtain the average number per site using our derivative form from Eq. (9) as Taking the derivative and substituting the expression for $Z_{\mathrm{NI}}$ in the denominator, we have where the derivative of the activity is Using this in the expression for the mean site occupancy in Eq. (19) and cancelling ${\mathcal{N}}_{\mathrm{sites}}-1$ factors of the sum over $\gamma$, the mean occupancy for species $\gamma$ is given by In Figure 5, we show the mean occupancies for the three species versus ${\epsilon }_{\Vert }$, with ${\epsilon }_{\Vert }=2{\epsilon }_{\perp }$. Negative ${\epsilon }_{\Vert }$ corresponds to an attraction of dimers to the lattice, and we see that parallel adsorption of dimers dominates because of the stronger binding, followed by perpendicular adsorption, with vacancies becoming nonexistent. Positive ${\epsilon }_{\Vert }$ corresponds to a repulsion of dimers from the lattice, so that at large ${\epsilon }_{\Vert }$, the ordering is reversed, for the same reasons. At small positive ${\epsilon }_{\Vert }$, the dimers are repelled from the lattice, but this competes with the chemical potential, which tries to put dimers back onto the lattice. In this low-coverage regime, perpendicular adsorption dominates, while vacancies dominate for large $\beta {\epsilon }_{\Vert }$ because the dimers would prefer to stay in solution. For ${\epsilon }_{\Vert },{\epsilon }_{\perp }\to 0$, the parallel and perpendicular mean occupancies become the same. Similarly, where ${\epsilon }_{\Vert }-{\mu }_{\mathrm{dimer}}=0$, ${\langle n^{(\Vert )}\rangle }_{\mathrm{NI}}={\langle n^{\left(v\right)}\rangle }_{\mathrm{NI}}$ because, as mentioned earlier, ${\epsilon }_v-{\mu }_{\mathrm{dimer}}$ is always zero. Also, ${\langle n^{(\perp )}\rangle }_{\mathrm{NI}}={\langle n^{\left(v\right)}\rangle }_{\mathrm{NI}}$ where ${\epsilon }_{\Vert }=2{\mu }_{\mathrm{dimer}}$ because that is where ${\epsilon }_{\perp }={\mu }_{\mathrm{dimer}}$. It is at this point where ${\langle n^{(\Vert )}\rangle }_{\mathrm{NI}}$ becomes a maximum.

The average charge per site can now be determined by multiplying $q_{\gamma }$, the charge of species $\gamma$, by $n_{\ell ,b}^{\left(\gamma \right)}$ the occupation of species $\gamma$ on site ℓ of chain b, Since we assume that every site has either a parallel dimer, a perpendicular dimer, or a vacancy, with single occupancy, then for a given site and species, $n_{\ell b}^{\left(\gamma \right)}$ is either 0 or 1. The average charge density is then the sum the averages of the individual terms, This charge density is plotted as the solid blue curve in Figure 6. Because a parallel dimer has no charge, $q_{\Vert }=0$, a perpendicular dimer has a positive unit charge, $q_{\perp }=+1$, and a vacancy has a negative unit charge, $q_v=-1$, the total mean field charge per site ${\rho }_c$ in Figure 6 is the difference between the curves ${\langle n^{(\perp )}\rangle }_{\mathrm{NI}}$ and ${\langle n^v\rangle }_{\mathrm{NI}}$ in Figure 5 and goes to zero where those two curves cross.

A measure of the magnitude of the charge fluctuations is given by the average of the square of the charge on a site minus the square of its average, which is the charge variance ${\sigma }_{NI}^2$, given by where and The product $n_{\ell b}^{\left(\gamma \right)}{n}_{\ell b}^{\left({\gamma }^{\prime }\right)}$ describes a double occupancy of site $(\ell ,b)$ by species $\gamma$ and ${\gamma }^{\prime }$. Since we have required single occupancy, then we must have $\gamma ={\gamma }^{\prime }$. Also, since $n_{\ell b}^{\left(\gamma \right)}$ can only take the values 0 or 1, Then the charge density per site reduces to Taking the thermal average of both sides, we have In Figure 6, we have plotted the standard deviation in the charge density, ${\sigma }_{\mathrm{NI}}$, which is the square root of the charge variance, for the noninteracting model. As we saw in Figure 5, there are three distinct regions, ${\epsilon }_{\Vert }\ll 0$, ${\epsilon }_{\Vert }\gg 0$ and ${\epsilon }_{\Vert }$ near ${\mu }_{\mathrm{dimer}}$. These three regions are also reflected in Figure 6. At large negative ${\epsilon }_{\Vert }$, parallel adsorption dominates, which leads to $\rho \approx 0$. At large positive ${\epsilon }_{\Vert }$, vacancies dominate, leading to $\rho <0$, and for ${\epsilon }_{\Vert }$ near ${\mu }_{\mathrm{dimer}}$, there is a small window where perpendicular adsorption of dimers dominates, leading to positive values of $\rho$. Correspondingly, the fluctuations, represented by the standard deviation ${\sigma }_{NI}$, become small when $|{\epsilon }_{\Vert }|$ becomes large, and the fluctuations are largest at small $|{\epsilon }_{\Vert }|$, when there are comparable numbers of all species. The phenomenon of charge inversion is demonstrated in Figure 6 because the average charge is positive, indicating that sufficiently many dimers adsorb in a perpendicular configuration to invert the charge on the molecule from negative to positive. The magnitude of charge inversion increases in the weak-binding limit ${\epsilon }_{\perp },{\epsilon }_{\Vert }\to 0$. The maximum of this inversion in Figure 6 is about 0.25, which is much larger than the estimate of 0.07 from the work of Nguyen and Shklovskii[26] , as discussed in the Introduction. That work treated the charge on DNA as a continuum, while we have used the discrete lattice of phosphates in a lattice-gas model. We have not included interactions between sites yet, although we have assumed that there can only be one species per site, either a parallel dimer, a perpendicular dimer, or a vacancy. We will see later in this paper that interactions reduce this value somewhat, but it will still be much larger than 0.07.

The noninteracting entropy $S_{\mathrm{NI}}$ can be obtained from the partition function using Eq. (7) as

Because the derivative of the activity with respect to $\beta$ can be written in terms of its logarithm as the entropy can be written in the simple form Substituting the mean occupancy for the noninteracting model from Eq. (21) allows us to write the dimensionless entropy per site for the noninteracting model in a simplified form as which agrees with the standard result for the entropy of mixing of an ideal solution with species $\gamma =(\Vert ,\perp ,v)$ [40].

In Figure 7, we show the entropy $S_{NI}$ of the noninteracting model as a function of the binding energy $\beta {\epsilon }_{\Vert }$, assuming ${\epsilon }_{\Vert }=2{\epsilon }_{\perp }$. We also show the individual contributions of Eq. (33) to the entropy. The entropy is a maximum when disorder is greatest, and this occurs when the numbers of each of the species are as close to equal as possible, which occurs at ${\epsilon }_{\Vert }={\mu }_{\mathrm{dimer}}$, the maximum of $\langle n^{(\perp )}\rangle$ in Figure 5.

While stored in the nucleus of a cell, DNA is wrapped compactly both around histone protein complexes and around itself, but on sufficently small scales (≈ 150 base pairs [41], or 15 turns [30], or 50 nm [8]), DNA’s dominant geometrical structure is the familiar double-helix structure shown in Figure 2. Each strand of the DNA takes the shape of a helix with a characteristic radius $R_{\mathrm{DNA}}=0.946$ nm and pitch angle $\psi =29.6^{\circ }$, as shown in Figure 8. Here the pitch angle $\psi$ denotes the angle with respect to the $xy$-plane that gives the appropriate altitude per unit circumferential winding; in cylindrical coordinates, $tan\psi =\Delta z/R_{\mathrm{DNA}}\Delta \varphi$, where $\Delta z$ is the distance along the z axis. Each strand of DNA has a “direction", identified by a particular carbon on the backbone structure, corresponding to the chirality of the helix. In the DNA double-helix, the two strands are antiparallel and therefore have opposite chiralities. As a consequence, the azimuthal angle between the two helices is always a constant, $\Delta \varphi ={160}^{\circ }$. Because this phase shift is not exactly ${180}^{\circ }$, the chains have unequal separation in the clockwise and counterclockwise azimuthal directions. The larger gap is referred to as the major groove, and the smaller as the minor groove.

When the hydrogen atoms dissociate under physiological conditions, the resulting negative charges ($-e$, where e is the magnitude of the charge of the electron), can be regarded as located at the mean position of the oxygen atoms on the backbone, as shown in Figure 2. This figure was constructed using geometrical data taken from Bishop and McMullen [28], which is based on experimental x-ray diffraction data [42,43] as input to the SYBYL molecular modeling program[44]. The oxygen atoms occur at regular intervals along each strand, separated by a helix segment of arc length $a^{\left(1\right)}=0.684$ nm. It is these sites located at regular intervals along the helix to which dimers will adsorb. These negatively-charged sites do not occur at the same altitudes on both strands, however. Rather, there is a vertical separation $\Delta z=0.023$ nm between corresponding sites on the two strands. With the relative phase of the strands and the vertical separation between sites on those strands taken together, corresponding sites on the two strands may be viewed as connected by a helical segment of arc length $a^{\left(2\right)}=3.34$ nm at a pitch angle of $\alpha =0.{394}^{\circ }$. This geometry is shown in Figure 8.

When positively-charged dimers approach the DNA molecule, they will be attracted to the negative charges at the sites on the double-helix. We consider dimers with a length comparable to the spacing $a^{\left(1\right)}$ between sites on a strand and having positive charges $+e$ at either end. The dimers can then adsorb onto the surface in two possible orientations, either parallel to the strand or perpendicular to the helix axis. (See Figure 4.) If the dimer adsorbs parallel to the strand, the positive charges from the dimer lie directly over the negative charges on the strand, neutralizing the charge on two adjacent sites. If the dimer adsorbs perpendicular to the surface of the bounding cylinder, one end of the dimer sits atop the site, while the other extends radially outward. This perpendicular adsorption effectively inverts the charge on the site from $-e$ to $+e$. In order to use a lattice-gas model, this geometric constraint is loosened by having the parallel dimer block only a single site. This deficiency can be somewhat compensated by making binding energy of the parallel dimer twice as large, ${\epsilon }_{\Vert }\approx 2{\epsilon }_{\perp }$.

Note that because the length of the dimer (equal to the same-chain site spacing $a^{\left(1\right)}=0.684$ nm) is much smaller than both the cross-chain site spacing $a^{\left(2\right)}=3.34$ nm and the vertical separation $\Delta z=3.4$ nm between turns of the helix, other orientations of the dimer are not possible.

The problem thus described is a complex one, but the similarities with the lattice gas models of condensed matter physics provide guidelines for how to proceed. These prescriptions, however, are aimed at the treatment of a periodic crystalline lattice, and, although the DNA sites exhibit helical symmetry, they do not constitute a periodic lattice in the strict sense of the term. However, an appropriate choice of coordinates can take advantage of the helical symmetry, so that, in these new coordinates, the positions of the sites will fall on a regular, one-dimensional lattice.

We will define these coordinates on a cylinder of radius $R_{DNA}$, as shown in Figure 8 and Figure 9. The first coordinate $x^{\left(1\right)}$ traces out a path with pitch angle $\psi$ along a single helical strand, and the other coordinate $x^{\left(2\right)}$ traces out a path with pitch angle $\alpha$ that connects corresponding sites on the two strands. Geometrically, a cylinder can be regarded as flat in the sense that it has no curvature. In Figure 9, we show the way that the cylinder can be cut with scissors and unwrapped so that this lattice can be mapped on a flat surface as shown in Figure 10. If we define the origin of coordinates $x^{\left(1\right)}$ and $x^{\left(2\right)}$ to be halfway along the helical path between the partners on the two chains, as shown in Figure 9 and Figure 10, then the positions of the sites on both strands form a one-dimensional lattice in the coordinates $(x^{\left(1\right)},x^{\left(2\right)})$.

These coordinates can be written simply in terms of cylindrical coordinates $(\varphi ,z)$ in matrix form as Inverting this gives definitions of the two coordinates as and With these definitions, the difference in coordinates between adjacent sites on the same strand is $\Delta x^{\left(1\right)}=a^{\left(1\right)}$, and the difference in coordinates between corresponding sites on the two strands is $\Delta x^{\left(2\right)}=a^{\left(2\right)}$. That is, $a^{\left(1\right)}$ is the distance along the helical path of a single chain from one phosphate ion to the next, and $a^{\left(2\right)}$ is the distance along a helical path from a phosphate ion on one chain to its partner phosphate ion on the other chain.

Next we define a lattice index ℓ, which specifies the cell (altitude on the double-helix) and chain index b, which specifies the basis site, where $b=-\frac{1}{2}$ for the “down"(↓) chain and $b=\frac{1}{2}$ for the “up"(↑) chain, as shown in Figure 11. Using these variables, the coordinates can be written as where Thus, although the sites on the DNA molecule do not constitute a periodic lattice in real space, they do constitute a lattice in an appropriately-defined coordinate space (see Figure 11). As we will see, however, this choice of coordinates will make the form of the interaction potential more complicated as a result.

The use of $(x^{\left(1\right)},x^{\left(2\right)})$ instead of the cylindrical coordinates $(\varphi ,z)$ indicates a more fundamental shift in our description of the DNA double-helix. The two-dimensional surface on which the helices lie is a cylinder of radius $R_{DNA}$, and the helices inherit the cylinder’s geometric properties. The geometry of the cylinder, however, is locally indistinguishable from the geometry of the flat plane. One common consequence of this is that it is possible to smoothly wrap a flat sheet of paper around a cylinder. In contrast, it is not possible to smoothly wrap a sheet of paper around a sphere; this problem is well-known because of the geometrical distortions that occur in flat maps of a spherical Earth. Maps of a cylindrical surface, however, have no such distortions. This geometrical difference is quantified by the Riemannian curvature tensor, which vanishes for both the cylinder and the plane, but not for the sphere[45]. This means that the local geometry of the cylinder behaves in exactly the same way as the local geometry of the plane, so that a helix on a cylinder is geometrically equivalent to a line on a plane.

Our choice of coordinates is simply a map of the cylindrical surface that reduces the double-helix to two parallel lines, as shown in Figure 11. The result of this map is that we have a linear lattice with two phosphate sites per cell, with the cells labeled from $\ell =-\mathcal{N}$ to $\mathcal{N}$. All $2\mathcal{N}+1$ cells are identical, and the lattice can be imagined to satisfy periodic boundary conditions.

The analogy describing DNA as a ladder wrapped around a cylinder of radius $R_{\mathrm{DNA}}$ then has a true mathematical basis, because the local structure of the double-helix is equivalent to the structure of a “ladder" — a one-dimensional chain with a unit cell containing one site from each strand. If interactions are ignored, then the problem is described simply by this one-dimensional lattice.

The geometrical structure of the double helix is important in determining the electrostatic energy $U_{\mathrm{int}}$ of the sites (with and without dimers) interacting with one another. The full Hamiltonian must contain these contributions and can be written as where the interaction energy $U_{\mathrm{int}}$ is the total interaction energy between all the charges on all the sites. If ${\mathcal{V}}_{({\ell }_1,b_1),({\ell }_2,b_2)}$ is the screened Coulomb energy between a charge $q_{{\gamma }_1}$ at site ${\ell }_1$ on chain $b_1$ with another charge $q_{{\gamma }_2}$ at site ${\ell }_2$ on chain $b_2$, then $U_{\mathrm{int}}$ can be written as where $n_{\ell ,b}^{\left(\gamma \right)}$ is the number of particles of species $\gamma$ on the lattice site at $(\ell ,b)$ and $q_{\gamma }$ is the charge of that species, in units of the magnitude e of the charge of an electron. The factor of $1/2$ is to ensure that we are not double counting when we sum over all lattice sites, and we implicitly exclude same-site interactions $({\ell }_1,b_1)\ne ({\ell }_2,b_2)$. The total charge on a given site is given by summing all the charges on that site, so that the total charge on the site $(\ell ,b)$ is The interaction potential can then be written in terms the the total charge on each site as where $\mathcal{V}$ is the electrostatic interaction energy between the sites. This screened electrostatic energy between charges $q_{{\gamma }_1}$ and $q_{{\gamma }_2}$ at positions $({\ell }_1,b_1)$ and $({\ell }_2,b_2)$, (in SI units) is where $ϵ$ is the electric permittivity of the medium between the two charges, and $d_{b_1{b}_2}({\ell }_1-{\ell }_2)$ is the straight-line distance between the two charges. Distances along the chain are invariant under translations by a lattice spacing, and so $d_{b_1{b}_2}({\ell }_1-{\ell }_2)$ depends only on the difference between the lattice site indices ${\ell }_1-{\ell }_2$. As in the noninteracting lattice-gas model, we have three species of “particles" on the lattice sites, vacancies of charge $q_v=-1$, parallel (‖) dimers of charge $q_{\Vert }=0$, and perpendicular (⊥) dimers of charge $q_{\perp }=+1$.

Under physiological conditions, the presence of monovalent salt ions such as Na${}^{+}$ leads to screening of the bare charges. Traditionally, screening is treated by modeling the ions as a continuous density, resulting in the nonlinear Poisson-Boltzmann equation [19,46]. The Poisson-Boltzmann equation is not analytically solvable, so it is often further approximated by linearization. The resulting Thomas-Fermi model is analytically solvable and gives the screened Coulomb (or Yukawa) potential between two charges given in Eq. (43) where $ϵ=78.5{ϵ}_0$ is the permittivity of water (with ${ϵ}_0$ the permittivity of free space) and $q_s$ is the magnitude of the screening wave vector [47]. Various names are ascribed to both the nonlinear and linearized equations, including Debye-Hückel, Thomas-Fermi, and Poisson-Boltzmann, but we will always refer to the Poisson-Boltzmann equation when we mean the nonlinear form and the (linearized) Thomas-Fermi equation when we mean the linearized form.

One must be cautious in using the screened Coulomb potential for screened interactions. The continuous density approximation from which the nonlinear Poisson-Boltzmann equation was derived constitutes a form of mean-field theory[19], which fails to describe any of the effects due to correlations between the ions such as charge inversion and condensation. Thus we cannot model screening of the DNA molecule by the dimers using the Thomas-Fermi model, or even the Poisson-Boltzmann equation. Instead, we must treat the dimers as individual particles and compute their interactions so that we do not ignore correlation effects.

For the screening due to the monovalent salt, however, the valence involved is small enough and the concentration is low enough that a mean-field treatment is generally believed to be a good approximation [19,48]. Thus, we will treat the interaction energy $U_{\mathrm{int}}$ in Eq. (42) using the screened Coulomb potential for the dimer-dimer interactions, with the screening vector $q_s$ given in the Thomas-Fermi model as [47] where $\beta =1/k_{B}T$ is the inverse temperature, n is the concentration of ions, and ${\ell }_B=\beta e^2/\left(4\pi ϵ\right)$ is the Bjerrum length [32]. The Bjerrum length is a characteristic length scale equal to the distance at which two proton charges interact with energy $k_{B}T$. It should be noted that these expressions differ slightly in the literature because of various authors’ choices of units for the electrostatics. Here, we use SI units.

By choosing a coordinate system based on the arc lengths around the surface of the cylinder, we have managed to align the sites into a regular lattice, but the potential ${\mathcal{V}}_{({\ell }_1,b_1),({\ell }_2,b_2)}$ depends on the straight-line distance between the sites rather than the surface arc-length connecting them. We must then express Eq. (43) in terms of the straight-line distance function $d_{b_1{b}_2}({\ell }_1-{\ell }_2)$ between the sites at $({\ell }_1{a}^{\left(1\right)},b_1{a}^{\left(2\right)})$ and $({\ell }_2{a}^{\left(1\right)},b_2{a}^{\left(2\right)})$. The straight-line distance in cylindrical coordinates between two points on the surface of a cylinder of radius $R_{\mathrm{DNA}}$ is and applying our change of variables (35,36) gives the straight-line distance between sites as where $\ell ={\ell }_1-{\ell }_2$ and $b=b_1-b_2$. Because of the translational invariance of the lattice, the distance in Eq. (46) depends only on the differences between the lattice and basis indices, and not their values individually. From this relation, we see the symmetry property of the distance,

Since $b_1$ and $b_2$ are either $+\frac{1}{2}$ or $-\frac{1}{2}$, $b=b_1-b_2$ takes the three values $b=0$ for interactions of one site on a single chain with all other sites on the same chain, $b=1$ for interactions between a site on the “up" chain with with all the sites on the “down" chain, and $b=-1$ for interactions of one site on the “down" chain with all the sites on the “up" chain. We will write an up arrow ↑ will indicate $b_{1,2}=\frac{1}{2}$, and a down arrow ↓ will indicate $b_{1,2}=-\frac{1}{2}$. Plots of $d_{\uparrow \uparrow }(\ell )/R_{\mathrm{DNA}}$ and $d_{\uparrow \downarrow }(\ell )/R_{\mathrm{DNA}}$ are shown in Figure 12a. There are wiggles in the curve because, as the chain wraps around the cylinder (see Figure 4), the distances between lattice sites can become larger or smaller than they would be on a linear chain.

In order to understand the consequences of the screened Coulomb potential (43) depending on the straight-line distances d between sites, we define a new notation for the potential, which is expressed in terms of the difference in lattice site indices ${\ell }_1-{\ell }_2$, This potential decays exponentially as a function of the straight-line distance d (inset of Figure 12(b)). However, as a function of the lattice index ℓ as in Figure 12(b), there are corresponding “wiggles" in the decay. This occurs because the distance between the sites decreases slightly as the helix completes a turn. These wiggles are present for every turn the helix makes, but the effect on the potential becomes negligibly small after about the first turn. Thus, by using the indices ℓ and b to describe the relationships between sites, we have reduced the DNA double-helix structure to a regular one-dimensional lattice with two sites per cell at the cost of an irregular potential due to the geometrical structure of the DNA molecule. The potentials for other combinations of $b_1$ and $b_2$ can be related to those in Figure 12(b) through the symmetry relations

Fourier transforming in the site index ℓ block diagonalizes this matrix into $2\times 2$ blocks, corresponding to the chain indices $b_1$ and $b_2$. Explicitly, the Fourier transform is where the elements of $\mathcal{F}$ are independent of $b_1$ and $b_2$ and are given by and the dagger denotes the adjoint (complex conjugate transpose). We use periodic boundary conditions, so that the wave vector k appearing here, which is dimensionless, takes the $2\mathcal{N}+1$ values where If $\mathcal{N}$ is large, the range of k is essentially continuous from $-\pi$ to $\pi$.

The matrix elements of this Fourier transform are Substituting for $\mathcal{V}$ from Eq. (48), we have

Since we regard the chain as having periodic boundary conditions, changing summation indices to $\ell ={\ell }_1-{\ell }_2$ and ${\ell }_2$ gives where the sum over ${\ell }_2$ becomes

We can now write the Fourier transform of $V_{b_1{b}_2}(\ell )$ as Substituting this back into the matrix, we have Thus, the Fourier transform $\tilde{\mathcal{V}}$ of $\mathcal{V}$ is diagonal in the k indices and contains $2\times 2$ blocks The symmetries (49) of the screened Coulomb potential on the DNA lattice are reflected in the Fourier transforms (58) as well: Also, because $V_{\uparrow \uparrow }(\ell )$ is an even function of ℓ, ${\tilde{V}}_{\uparrow \uparrow }\left(k\right)$ is real. The Fourier transforms (58) are plotted in Figure 13.

As can be seen in Figure 13, the Fourier transforms $\tilde{V}\left(k\right)$ have regions of k that are negative, which is a result of transforming with respect to our helical coordinate system. When the screened Coulomb potential (43) is Fourier transformed with respect to the straight-line distance d, the function is positive definite. We have instead Fourier transformed the potential a function of the lattice index ℓ, resulting in the “wiggles" in Figure 12(b) caused by the turns of the helix. These deviations from the exponential decay of $V\left(d\right)$ result in the deviations here from the positive-definite Fourier transform.

Using the symmetry operations in Eqs. (61a) and (61b) obeyed by the Fourier transforms allows us to simplify the $2\times 2$ blocks $\tilde{V}\left(k\right)$ in Eq. (60) to This matrix is easily diagonalized, yielding the real eigenvalues which are plotted in Figure 14. The transformation that diagonalizes $\tilde{V}\left(k\right)$ is a unitary matrix $\xi \left(k\right)$ whose columns are the eigenvectors of $\tilde{V}\left(k\right)$. Choosing the ${\lambda }_{k+}$ eigenvector to be first, this transformation matrix is given by This diagonalizes the $2\times 2$ matrix $\tilde{V}\left(k\right)$ according to

We can write the combined process of Fourier transforming $\mathcal{V}$ into block-diagonal form and then diagonalizing the $2\times 2$ blocks $\tilde{V}\left(k\right)$ as a single unitary transformation $\mathcal{M}$ given by which completely diagonalizes $\mathcal{V}$ such that Here $\Lambda$ is a diagonal matrix, with each element along the diagonal given by: where $\sigma$ is either + or −.

Although we have diagonalized the potential with a unitary transformation, the fact that the transformation matrix is complex means that the charge density in the diagonal basis could also be complex, which is physically undesirable. However, the potential matrix in the position basis is real and symmetric, which means that it should be diagonalizable by a real orthogonal transformation matrix $\mathcal{W}$. In fact, the transformation matrix is not unique, since there is a degeneracy in the eigenvalues, since ${\lambda }_{k,\sigma }={\lambda }_{-k,\sigma }$. This means that we can choose eigenvectors that are real by taking two orthogonal linear combinations of the eigenvectors of the unitary transformation $\mathcal{M}$ for k and $-k$. Those two new eigenvectors, which represent the columns $(k,\sigma )$ for $k>0$ and $k<0$ of the new transformation matrix $\mathcal{W}$, are given by and Since ${\xi }_{b\sigma }(-k)={\xi }_{b\sigma }^{*}\left(k\right)$, these can be written in terms of real and imaginary parts as and When $k=0$, the orthogonal transformation is the same as the unitary transformation, that is where Here the columns are labeled by the eigenvalue $\sigma =\pm$, and the rows are labeled by the chain $b=-1/2$ for the first row and $b=1/2$ for the second row.

This transformation can be written in a more compact form by introducing the unitary transformation matrix $\mathcal{Z}$, whose elements are defined by Then the orthogonal transformation is

In deriving the thermodynamic quantities for the interacting lattice gas model, it will be convenient to be able to use this orthogonal transformation to diagonalize the potential.

For the interacting lattice gas, the partition function is similar to the noninteracting one in Eqs. (5) and (8). However, here the full Hamiltonian $\mathcal{H}$ from Eq. (39) is used, and then the partition function is given by The full Hamiltonian, including the interaction term from Eq. (42) written in matrix form, is In Section 4, we showed that the orthogonal transformation matrix $\mathcal{W}$ diagonalizes $\mathcal{V}$. In order to streamline the calculation, it will be useful to insert the identity matrix in the form $\mathcal{W}{\mathcal{W}}^T$, to the left and right of of $\mathcal{V}$ in the interaction term of the Hamiltonian, which gives where we have grouped factors to show the transformed quantities.

The Hamiltonian, with the interaction term written in the diagonal basis, can then be written as where $\tilde{\rho }={\mathcal{W}}^T\rho$ is the transformed charge density. We have not transformed the first term in the Hamiltonian to the diagonal basis because it is more convenient later in the calculation to have it in the position basis. Substituting this form into the partition function and writing the resulting expression as two separate exponentials, we have We can write the exponential of the sum in the interaction term as a product of exponentials as From this expression, we can see that the mean site occupancy for the interacting lattice-gas model can be obtained in the same way as in the noninteracting case (9), by differentiating the partition function with respect to the chemical potential.

In the noninteracting case, since all the information about the configuration was contained in $n_{\ell ,b}^{\left(\gamma \right)}$, we were able to decouple the sum over configurations from the sum over lattice sites because the activity was independent of position, and $n_{\ell ,b}^{\left(\gamma \right)}$ only appeared in its exponent. This is because only linear factors of $n_{\ell ,b}^{\left(\gamma \right)}$ were in the exponent in the partition function. While it is still true that $n_{\ell ,b}^{\left(\gamma \right)}$ contains the configuration dependence, it now appears quadratically in the exponent of the partition function, since ${\tilde{\rho }}_{k,\sigma }$ contains linear factors of the $n_{\ell ,b}^{\left(\gamma \right)}$, and the exponent in the interaction term contains ${\tilde{\rho }}_{k,\sigma }^2$. This makes it impossible to decouple the sum over configurations from the sum over lattice sites. However, we can use an integral identity, known as the Hubbard-Stratonovich transformation, to replace the term quadratic in ${\tilde{\rho }}_{k,\sigma }$ in the exponential with a term linear in ${\tilde{\rho }}_{k,\sigma }$, meaning that there will only be linear terms in $n_{\ell ,b}^{\left(\gamma \right)}$. This requires introducing the integral over an auxiliary field ${\tilde{\Delta }}_{k,\sigma }$. Then the sum over configurations is possible to do exactly, although at the cost of having to integrate over the auxiliary fields.

The integral identity, which is a complicated way of writing unity, as is shown in Appendix A, is given by where the integration path is over the real axis from $-\infty$ to ∞ when $\lambda <0$ and over the imaginary axis from $-i\infty$ to $i\infty$ when $\lambda >0$. In the previous work by Bishop and McMullen[28], the problem was done in the position basis, and this integral identity, known as the Hubbard-Stratonovich transformation, was written with the matrix version of the potential in the exponent, as given in Negele and Orland[48]. This form creates a dilemma when there are negative and possibly zero eigenvalues of the potential, and it is difficult to determine the path of integration, since it changes depending on the sign of the eigenvalues. Zero eigenvalues would make the determinant of the matrix in that formula zero, and one then finds that there is a division by zero in the formula. By transforming to the diagonal basis, all these difficulties are avoided, since there is a separate integral for each eigenvalue, and if a particular eigenvalue is zero, the identity is not used at all.

To use the identity in Eq. (84), we make the identifications that $\lambda ={\lambda }_{k,\sigma }$, $\tilde{\rho }={\tilde{\rho }}_{k,\sigma }$, and $\tilde{\Delta }={\tilde{\Delta }}_{k,\sigma }$. Substituting this “1" into the partition function, we have Now we see that the two exponentials containing ${\tilde{\rho }}_{k,\sigma }^2$ cancel, as we planned. Of course, we have gained this convenience by introducing the auxiliary field ${\tilde{\Delta }}_{k,\sigma }$, and we will have to do the integral over this variable at a later stage. However, since ${\tilde{\Delta }}_{k,\sigma }$ does not depend on the configuration of the system, we can bring the sum over configurations and the leading exponential factors in the partition function inside the integral, and the expression for the partition function becomes It will be useful to define an action $\mathcal{S}$ such that the partition function can be written in the form where and Note that the factors of $\beta {\lambda }_{k,\sigma }$ are included here in the grand differential $\mathcal{D}\left[\tilde{\Delta }\right]$, rather than incorporating them in the action $\mathcal{S}\left[\tilde{\Delta }\right]$, as was done in Bishop and McMullen[28] and in Sievert[33]. When these factors appear in the action, they introduce a term of the form $ln\sqrt{det\left(\beta v\right)}$. This adds a large constant value to the mean field results, which is subtracted out when the fluctuation terms are included. It actually has no real physical meaning and is part of the normalization of the integral[49].