Entropy of Charge Inversion in DNA including One-Loop Fluctuations

Abstract

The entropy and charge distributions have been calculated for a simple model of polyelectrolytes attached to the surface of DNA using a field-theoretic method that includes fluctuations to the lowest one-loop order beyond mean-field theory. Experiments have revealed correlation-driven behavior of DNA in charged solutions, including charge inversion and condensation. In our model, the condensed polyelectrolytes are taken to be doubly charged dimers of length comparable to the distance between sites along the phosphate chains. Within this lattice gas model, each adsorption site is assumed to have either a vacancy or a positively charged dimer attached with the dimer oriented either parallel or perpendicular to the double-helix DNA chain. We find that the inclusion of the fluctuation terms decreases the entropy by ∼50% in the weak-binding regime. There, the bound dimer concentration is low because the dimers are repelled from the DNA molecule, which competes with the chemical potential driving them from the solution to the DNA surface. Surprisingly, this decrease in entropy due to correlations is so significant that it overcompensates for the entropy increase at the mean-field level, so that the total entropy is even lower than in the absence of interactions between lattice sites. As a bonus, we present a transparent exposition of the methods used that could be useful to students and others wishing to use this formulation to extend this calculation to more realistic models.

1. Introduction

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 the 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 macro-ion/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 macro-ion [8]. This interplay between electrostatics and geometry is particularly important for polyelectrolytes in solution, and they have been 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 one 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 the 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 produced 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 rod-like particles in the presence of polyvalent counterions. Their first model consisted of two rod-like DNA molecules [20], and they found attraction between the rods. When Podgornik and Parsegian suggested that the fluctuation forces between such rod-like polyelectrolytes might not be pairwise additive [24], Ha and Liu extended their calculation to include bundles of rod-like 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 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 negative macro-ion 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 for a valence of the polyvalent cations greater than around three, the spacing of the 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 Equation (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 attached to the DNA either parallel or perpendicular to the DNA surface, and the model confirmed that under the right conditions, the charge on the surface could be inverted. The use of the lattice model allowed them to accommodate 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 polyelectrolyte 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 for 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 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 per site of all species, and the charge per site. 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 Appendices Appendix A and Appendix B.

2. The Noninteracting DNA Lattice Gas

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 is a point charge that 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 spermine 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,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 the 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 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 any of the 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 to 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 “particle” 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 $(\mathsf{\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

$${\mathcal{H}}_{\mathrm{NI}}=\sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}.$$

(1)

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],

$${\Omega }_G(T,V,\left\{{\mu }_{\gamma }\right\})=E-TS-\sum _{\gamma }{\mu }_{\gamma }{N}_{\gamma },$$

(2)

where the number of particles of type $\gamma$ is given by

$$N_{\gamma }=\sum _{\mathsf{\ell },b}{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}.$$

(3)

Then, the entropy can be written in terms of this potential as

$$S=-{\left(\frac{\partial {\Omega }_G}{\partial T}\right)}_{V,\left\{{\mu }_{\gamma }\right\}}=k_B{\beta }^2{\left(\frac{\partial {\Omega }_G}{\partial \beta }\right)}_{V,\left\{{\mu }_{\gamma }\right\}},$$

(4)

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

$$Z_G=\sum _{\mathrm{configurations}}{e}^{-\beta (E-{\sum }_{\gamma }{\mu }_{\gamma }{N}_{\gamma })}=e^{-\beta {\Omega }_G},$$

(5)

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

$${\Omega }_G=-\frac{ln{Z}_G}{\beta },$$

(6)

and this allows us to write an expression for the entropy in terms of $Z_G$ as

$$S=k_{B}ln{Z}_G-\frac{k_B\beta }{Z_G}{\left(\frac{\partial Z_G}{\partial \beta }\right)}_{V,\left\{{\mu }_{\gamma }\right\}}.$$

(7)

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

$$Z_{\mathrm{NI}}=\sum _{\mathrm{configurations}}{e}^{-\beta {\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }\left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)n_{\mathsf{\ell },b}^{\left(\gamma \right)}},$$

(8)

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 },$

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta Z_{\mathrm{NI}}}\frac{\partial Z_{\mathrm{NI}}}{\partial {\mu }_{\gamma }}.$$

(9)

This can be verified by taking the derivative of Equation (8).

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}{Z}_{\mathrm{NI}}}\sum _{\mathrm{configurations}}\sum _{\mathsf{\ell },b}{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}{e}^{-\beta {\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }\left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)n_{\mathsf{\ell },b}^{\left(\gamma \right)}},$$

(10)

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

$$Z_{\mathrm{NI}}=\sum _{\mathrm{configurations}}\prod _{\mathsf{\ell },b}\prod _{\gamma }{\left[e^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}\right]}^{n_{\mathsf{\ell },b}^{\left(\gamma \right)}}.$$

(11)

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

$$a_{\mathrm{NI}}^{\left(\gamma \right)}=e^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}.$$

(12)

Note that this is independent of the lattice site ℓ or chain b. The partition function can then be written in the simple form

$$Z_{\mathrm{NI}}=\sum _{\mathrm{configurations}}\prod _{\gamma }\prod _{\mathsf{\ell },b}{\left[a_{\mathrm{NI}}^{\left(\gamma \right)}\right]}^{n_{\mathsf{\ell },b}^{\left(\gamma \right)}}.$$

(13)

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 performed over each site separately, where there are three possible configurations

$$\{n_{\mathsf{\ell },b}^{(\perp )},n_{\mathsf{\ell },b}^{(\Vert )},n_{\mathsf{\ell },b}^{\left(v\right)}\}=\{1,0,0\},\{0,1,0\},\mathrm{or}\{0,0,1\}.$$

(14)

This is the same as saying that $n_{\mathsf{\ell },b}^{\left(\gamma \right)}=1$ for one and only one of $\gamma =\perp , \Vert$, or v, and is zero otherwise. This means that

$$\sum _{\begin{array}{c}\mathrm{single-site}\\ \mathrm{configurations}\end{array}}\prod _{\gamma =\perp ,\Vert ,v}{\left[a_{\mathrm{NI}}^{\left(\gamma \right)}\right]}^{n_{\mathsf{\ell },b}^{\left(\gamma \right)}}=\sum _{\gamma }{a}_{\mathrm{NI}}^{\left(\gamma \right)},$$

(15)

and so the grand partition function is

$$Z_{\mathrm{NI}}=\prod _{\mathsf{\ell },b}\left(\sum _{\gamma }{a}_{\mathrm{NI}}^{\left(\gamma \right)}\right).$$

(16)

Since every term in the product is the same, the expression in parentheses is simply raised to the power ${\mathcal{N}}_{\mathrm{sites}}$, giving

$$Z_{\mathrm{NI}}={\left(\sum _{\gamma }{a}_{\mathrm{NI}}^{\left(\gamma \right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}.$$

(17)

At this point, we can easily see that we can obtain the average number per site using our derivative form from Equation (9) as

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta Z_{\mathrm{NI}}}\frac{\partial }{\partial {\mu }_{\gamma }}{\left(\sum _{{\gamma }^{\prime }}{a}_{\mathrm{NI}}^{\left({\gamma }^{\prime }\right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}.$$

(18)

Taking the derivative and substituting the expression for $Z_{\mathrm{NI}}$ in the denominator, we have

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}=\frac{{\mathcal{N}}_{\mathrm{sites}}}{{\mathcal{N}}_{\mathrm{sites}}\beta {\left({\sum }_{{\gamma }^{″}}{a}_{\mathrm{NI}}^{\left({\gamma }^{″}\right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}}{\left(\sum _{{\gamma }^{\prime }}{a}_{\mathrm{NI}}^{\left({\gamma }^{\prime }\right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}-1}\frac{\partial }{\partial {\mu }_{\gamma }}{a}_{\mathrm{NI}}^{\left(\gamma \right)},$$

(19)

where the derivative of the activity is

$$\frac{\partial }{\partial {\mu }_{\gamma }}{a}_{\mathrm{NI}}^{\left(\gamma \right)}=\frac{\partial }{\partial {\mu }_{\gamma }}{e}^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}=\beta a_{\mathrm{NI}}^{\left(\gamma \right)}.$$

(20)

Using this in the expression for the mean site occupancy in Equation (19) and canceling ${\mathcal{N}}_{\mathrm{sites}}-1$ factors of the sum over $\gamma$, the mean occupancy for species $\gamma$ is given by

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}=\frac{a_{\mathrm{NI}}^{\left(\gamma \right)}}{{\sum }_{{\gamma }^{\prime }}{a}_{\mathrm{NI}}^{\left({\gamma }^{\prime }\right)}}.$$

(21)

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^{(\perp )}\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_{\mathsf{\ell },b}^{\left(\gamma \right)}$ the occupation of species $\gamma$ on site ℓ of chain b,

$${\rho }_{\mathsf{\ell }b}=\sum _{\gamma }{q}_{\gamma }{n}_{\mathsf{\ell }b}^{\left(\gamma \right)}.$$

(22)

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_{\mathsf{\ell }b}^{\left(\gamma \right)}$ is either 0 or 1. The average charge per site is then the sum of the averages of the individual terms,

$${\rho }_{\mathrm{NI}}\equiv {\langle {\rho }_{\mathsf{\ell }b}\rangle }_{\mathrm{NI}}=\sum _{\gamma }{q}_{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}.$$

(23)

This charge per site 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

$${\sigma }_{\mathrm{NI}}^2={\langle {\rho }^2\rangle }_{\mathrm{NI}}-{\rho }_{\mathrm{NI}}^2,$$

(24)

where

$${\langle {\rho }^2\rangle }_{\mathrm{NI}}={\langle {\rho }_{\mathsf{\ell }b}^2\rangle }_{\mathrm{NI}},$$

(25)

and

$${\rho }_{\mathsf{\ell }b}^2={\left[\sum _{\gamma }{q}_{\gamma }{n}_{\mathsf{\ell }b}^{\left(\gamma \right)}\right]}^2=\sum _{\gamma ,{\gamma }^{\prime }}{q}_{\gamma }{q}_{{\gamma }^{\prime }}{n}_{\mathsf{\ell }b}^{\left(\gamma \right)}{n}_{\mathsf{\ell }b}^{\left({\gamma }^{\prime }\right)}.$$

(26)

The product $n_{\mathsf{\ell }b}^{\left(\gamma \right)}{n}_{\mathsf{\ell }b}^{\left({\gamma }^{\prime }\right)}$ describes a double occupancy of site $(\mathsf{\ell },b)$ by species $\gamma$ and ${\gamma }^{\prime }$. Since we have required single occupancy, then we must have $\gamma ={\gamma }^{\prime }$. Also, since $n_{\mathsf{\ell }b}^{\left(\gamma \right)}$ can only take the values 0 or 1,

$${\left[n_{\mathsf{\ell }b}^{\left(\gamma \right)}\right]}^2=n_{\mathsf{\ell }b}^{\left(\gamma \right)}.$$

(27)

Then, the charge per site reduces to

$${\rho }_{\mathsf{\ell }b,\mathrm{NI}}^2=\sum _{\gamma }{q}_{\gamma }^2{n}_{\mathsf{\ell }b,\mathrm{NI}}^{\left(\gamma \right)}.$$

(28)

Taking the thermal average of both sides, we have

$${\langle {\rho }^2\rangle }_{\mathrm{NI}}\equiv {\langle {\rho }_{\mathsf{\ell }b}^2\rangle }_{\mathrm{NI}}=\sum _{\gamma }{q}_{\gamma }^2{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}.$$

(29)

In Figure 6, we have plotted the standard deviation in the charge per site, ${\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 Equation (7) as

$$S_{\mathrm{NI}}=k_{B}ln\left[{\left(\sum _{\gamma }{a}_{\mathrm{NI}}^{\left(\gamma \right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}\right]-\frac{k_B\beta }{{\left({\sum }_{{\gamma }^{\prime }}{a}_{\mathrm{NI}}^{\left({\gamma }^{\prime }\right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}}\frac{\partial }{\partial \beta }{\left(\sum _{\gamma }{a}_{\mathrm{NI}}^{\left(\gamma \right)}\right)}^{{\mathcal{N}}_{\mathrm{sites}}}.$$

(30)

Because the derivative of the activity with respect to $\beta$ can be written in terms of its logarithm as

$$\beta \frac{\partial }{\partial \beta }{a}_{\mathrm{NI}}^{\left(\gamma \right)}=a_{\mathrm{NI}}^{\left(\gamma \right)}ln{a}_{\mathrm{NI}}^{\left(\gamma \right)},$$

(31)

the entropy can be written in the simple form

$$\frac{S_{\mathrm{NI}}}{k_B{\mathcal{N}}_{\mathrm{sites}}}=-\sum _{\gamma }\frac{a_{\mathrm{NI}}^{\left(\gamma \right)}}{{\sum }_{{\gamma }^{″}}{a}_{\mathrm{NI}}^{\left({\gamma }^{″}\right)}}ln\left(\frac{a_{\mathrm{NI}}^{\left(\gamma \right)}}{{\sum }_{{\gamma }^{\prime }}{a}_{\mathrm{NI}}^{\left({\gamma }^{\prime }\right)}}\right).$$

(32)

Substituting the mean occupancy for the noninteracting model from Equation (21) allows us to write the dimensionless entropy per site for the noninteracting model in a simplified form as

$$\frac{S_{\mathrm{NI}}}{k_B{\mathcal{N}}_{\mathrm{sites}}}=-\sum _{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}ln{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}},$$

(33)

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 Equation (33) to the entropy. The entropy is a maximum when the 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.

3. Geometry of the Charged Double Helix of DNA

While stored in the nucleus of a cell, DNA is wrapped compactly both around histone protein complexes and around itself, but on sufficiently 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. As explained in Bishop and McMullen [28], that structure can be constructed by wrapping a flexible ladder around a cylinder such that the ladder lies flat against the surface of the cylinder. Then, the lattice sites are the intersections between the rungs and the sides of the ladder. The position of each site is the mean position of the pair of protruding oxygens that are attached to each phosphorus, which incorporates the resonant structure represented by one electron being shared by two oxygens. These sites are considered to be the locations of negative point charges ($-e$, where e is the magnitude of the charge of the electron). The locations of the protruding oxygens are those given by the SYBYL molecular modeling software, version 6.9.2 (Tripos, Inc., St. Louis, MO, USA: www.tripos.com accessed on 15 June 2006) [42] for B-form DNA, which uses X-ray data [43,44,45]. This form of DNA was used because it is the most relevant physiologically. In the model shown, we used all adenine–thymine base pairs. The SYBYL program used a period of exactly ten base pairs for one complete turn of the helix, which means that the relative rotation angle between adjacent pairs was exactly $\Delta \varphi ={36}^{\circ }$, or $1/10$ of ${360}^{\circ }$. Each strand of the DNA then 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.

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 for by making the 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 the cylindrical coordinates $(\varphi ,z)$ in matrix form as

$$\left[\begin{array}{c}{R}_{\mathrm{DNA}}\varphi \\ z\end{array}\right]=\left[\begin{array}{cc}cos\psi & cos\alpha \\ sin\psi & sin\alpha \end{array}\right]\left[\begin{array}{c}{x}^{\left(1\right)}\\ x^{\left(2\right)}\end{array}\right].$$

(34)

Inverting this gives definitions of the two coordinates as

$$x^{\left(1\right)}=\left(\frac{sin\alpha }{sin(\alpha -\psi )}\right)R_{\mathrm{DNA}}\varphi -\left(\frac{cos\alpha }{sin(\alpha -\psi )}\right)z$$

(35)

and

$$x^{\left(2\right)}=-\left(\frac{sin\psi }{sin(\alpha -\psi )}\right)R_{\mathrm{DNA}}\varphi +\left(\frac{cos\psi }{sin(\alpha -\psi )}\right)z.$$

(36)

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

$$(x^{\left(1\right)},x^{\left(2\right)})=(\mathsf{\ell }{a}^{\left(1\right)},b{a}^{\left(2\right)}),$$

(37)

where

$$\mathsf{\ell }=0,\pm 1,\pm 2,\dots ;b=\pm \frac{1}{2}.$$

(38)

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 [46]. 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 $\mathsf{\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.

4. DNA Interaction Potential

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. Any thermal fluctuations in the positions of the charged lattice sites are omitted for simplicity. The full Hamiltonian must contain these contributions and can be written as

$$\mathcal{H}=\sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}+U_{\mathrm{int}},$$

(39)

where the interaction energy $U_{\mathrm{int}}$ is the total interaction energy between all the charges on all the sites. If ${\mathcal{V}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}$ is the screened Coulomb energy between a charge $q_{{\gamma }_1}$ at site ${\mathsf{\ell }}_1$ on chain $b_1$ with another charge $q_{{\gamma }_2}$ at site ${\mathsf{\ell }}_2$ on chain $b_2$, then $U_{\mathrm{int}}$ can be written as

$$U_{\mathrm{int}}=\frac{1}{2}\sum _{{\gamma }_1,{\gamma }_2}\sum _{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}{q}_{{\gamma }_1}{n}_{{\mathsf{\ell }}_1,b_1}^{\left({\gamma }_1\right)}{\mathcal{V}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}{n}_{{\mathsf{\ell }}_2,b_2}^{\left({\gamma }_2\right)}{q}_{{\gamma }_2},$$

(40)

where $n_{\mathsf{\ell },b}^{\left(\gamma \right)}$ is the number of particles of species $\gamma$ on the lattice site at $(\mathsf{\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 $({\mathsf{\ell }}_1,b_1)\ne ({\mathsf{\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 $(\mathsf{\ell },b)$ is

$${\rho }_{\mathsf{\ell },b}=\sum _{\gamma }{q}_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}.$$

(41)

The interaction potential can then be written in terms of the total charge on each site as

$$U_{\mathrm{int}}=\frac{1}{2}\sum _{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}{\rho }_{({\mathsf{\ell }}_1,b_1)}{\mathcal{V}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}{\rho }_{({\mathsf{\ell }}_2,b_2)},$$

(42)

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 $({\mathsf{\ell }}_1,b_1)$ and $({\mathsf{\ell }}_2,b_2)$ (in SI units) is

$${\mathcal{V}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}=\frac{e^2}{4\pi ϵ}\frac{e^{-q_s{d}_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2)}}{d_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2)},$$

(43)

where $ϵ$ is the electric permittivity of the medium between the two charges, and $d_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\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}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2)$ depends only on the difference between the lattice site indices ${\mathsf{\ell }}_1-{\mathsf{\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,47]. 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 Equation (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 [48]. 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,49]. Thus, we will treat the interaction energy $U_{\mathrm{int}}$ in Equation (42) using the screened Coulomb potential for the dimer–dimer interactions, with the screening vector $q_s$ given in the Thomas–Fermi model as [48]

$$q_s=\sqrt{2n\frac{\beta e^2}{ϵ}}=\sqrt{8\pi n{\mathsf{\ell }}_B},$$

(44)

where n is the concentration of ions and ${\mathsf{\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. We assume that the DNA is submerged in a 0.1 molar NaCl solution and that the dielectric constant is $\kappa =\frac{ϵ}{{ϵ}_0}=78.5$ and the temperature is T = 37 ${}^{\circ }$C = 310 K. For these numbers, the Bjerrum length is ${\mathsf{\ell }}_B=0.79$ nm and the screening length is $q_s^{-1}=0.98$ nm.

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}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\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 Equation (43) in terms of the straight-line distance function $d_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2)$ between the sites at $({\mathsf{\ell }}_1{a}^{\left(1\right)},b_1{a}^{\left(2\right)})$ and $({\mathsf{\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

$$\mathrm{distance}({\varphi }_1,z_1,{\varphi }_2,z_2)={\left\{2R_{\mathrm{DNA}}^2\left[1-cos\left({\varphi }_1-{\varphi }_2\right)\right]+{\left(z_1-z_2\right)}^2\right\}}^{\frac{1}{2}},$$

(45)

and applying our change in variables (35) and (36) gives the straight-line distance between sites as

$$\begin{array}{ccc}\hfill d_{b_1{b}_2}\left(\mathsf{\ell }\right)& =& R_{\mathrm{DNA}}\left\{2\left[1-cos\left(\frac{a^{\left(1\right)}}{R_{\mathrm{DNA}}}\mathsf{\ell }cos\psi +\frac{a^{\left(2\right)}}{R_{\mathrm{DNA}}}bcos\alpha \right)\right]\right.\hfill \\ & & +{\left.{\left[\frac{a^{\left(1\right)}}{R_{\mathrm{DNA}}}\mathsf{\ell }sin\psi +\frac{a^{\left(2\right)}}{R_{\mathrm{DNA}}}bsin\alpha \right]}^2\right\}}^{\frac{1}{2}},\hfill \\ \hfill \end{array}$$

(46)

where $\mathsf{\ell }={\mathsf{\ell }}_1-{\mathsf{\ell }}_2$ and $b=b_1-b_2$. Because of the translational invariance in the lattice, the distance in Equation (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,

$$\begin{array}{cc}\hfill d_{b_1{b}_2}\left(\mathsf{\ell }\right)& =d_{b_2{b}_1}(-\mathsf{\ell })\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill d_{\uparrow \uparrow }\left(\mathsf{\ell }\right)& =d_{\downarrow \downarrow }(-\mathsf{\ell }).\hfill \end{array}$$

(48)

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 use an up arrow ↑ to indicate $b_{1,2}=\frac{1}{2}$, and a down arrow ↓ to indicate $b_{1,2}=-\frac{1}{2}$. Plots of $d_{\uparrow \uparrow }\left(\mathsf{\ell }\right)/R_{\mathrm{DNA}}$ and $d_{\uparrow \downarrow }\left(\mathsf{\ell }\right)/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 ${\mathsf{\ell }}_1-{\mathsf{\ell }}_2$,

$$V_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2)\equiv {\mathcal{V}}_{({\mathsf{\ell }}_1,b_1),({\mathsf{\ell }}_2,b_2)}.$$

(49)

This potential decays exponentially as a function of the straight-line distance d (inset of Figure 12b). However, as a function of the lattice index ℓ, as in Figure 12b, 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 12b through the symmetry relations

$$\begin{array}{cc}\hfill V_{b_1{b}_2}\left(\mathsf{\ell }\right)& =V_{b_2{b}_1}(-\mathsf{\ell })\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill V_{\uparrow \uparrow }\left(\mathsf{\ell }\right)& =V_{\downarrow \downarrow }(-\mathsf{\ell }).\hfill \end{array}$$

(51)

Fourier transforming 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

$$\tilde{\mathcal{V}}={\mathcal{F}}^{†}\mathcal{V}\mathcal{F},$$

(52)

where the elements of $\mathcal{F}$ are independent of $b_1$ and $b_2$ and are given by

$${\mathcal{F}}_{(\mathsf{\ell },b_1),(k,b_2)}=\frac{1}{\sqrt{2\mathcal{N}+1}}{e}^{ik\mathsf{\ell }},$$

(53)

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

$$k=-\mathcal{N}\Delta k\cdots \mathcal{N}\Delta k,$$

(54)

where

$$\Delta k=\left(\frac{2\pi }{2\mathcal{N}+1}\right).$$

(55)

If $\mathcal{N}$ is large, the range of k is essentially continuous from $-\pi$ to $\pi$.

The matrix elements of this Fourier transform are

$${\tilde{\mathcal{V}}}_{(k_1,b_1),(k_2,b_2)}=\frac{1}{2\mathcal{N}+1}\sum _{{\mathsf{\ell }}_1,{\mathsf{\ell }}_2=-\mathcal{N}}^{\mathcal{N}}{e}^{-i{k}_1{\mathsf{\ell }}_1}{\mathcal{V}}_{\left({\mathsf{\ell }}_1{b}_1\right),({\mathsf{\ell }}_2,b_2)}{e}^{i{k}_2{\mathsf{\ell }}_2}.$$

(56)

Substituting for $\mathcal{V}$ from Equation (49), we have

$${\tilde{\mathcal{V}}}_{(k_1,b_1),(k_2,b_2)}=\frac{1}{2\mathcal{N}+1}\sum _{{\mathsf{\ell }}_1,{\mathsf{\ell }}_2=-\mathcal{N}}^{\mathcal{N}}{e}^{-i(k_1{\mathsf{\ell }}_1-k_2{\mathsf{\ell }}_2)}{V}_{b_1{b}_2}({\mathsf{\ell }}_1-{\mathsf{\ell }}_2).$$

(57)

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

$${\tilde{\mathcal{V}}}_{(k_1,b_1),(k_2,b_2)}=\frac{1}{2\mathcal{N}+1}\sum _{\mathsf{\ell }=-\mathcal{N}}^{\mathcal{N}}{e}^{-i{k}_1\mathsf{\ell }}{V}_{b_1{b}_2}\left(\mathsf{\ell }\right)\sum _{{\mathsf{\ell }}_2=-\mathcal{N}}^{\mathcal{N}}{e}^{-i(k_1-k_2){\mathsf{\ell }}_2},$$

(58)

where the sum over ${\mathsf{\ell }}_2$ becomes

$$\sum _{{\mathsf{\ell }}_2=-\mathcal{N}}^{\mathcal{N}}{e}^{-i(k_1-k_2){\mathsf{\ell }}_2}=(2\mathcal{N}+1){\delta }_{k_1,k_2}.$$

(59)

We can now write the Fourier transform of $V_{b_1{b}_2}\left(\mathsf{\ell }\right)$ as

$${\tilde{V}}_{b_1{b}_2}\left(k\right)=\sum _{\mathsf{\ell }=-\mathcal{N}}^{\mathcal{N}}{e}^{-ik\mathsf{\ell }}{V}_{b_1{b}_2}\left(\mathsf{\ell }\right).$$

(60)

Substituting this back into the matrix, we have

$${\tilde{\mathcal{V}}}_{(k_1,b_1),(k_2,b_2)}={\delta }_{k_1,k_2}{\tilde{V}}_{b_1{b}_2}\left(k_1\right).$$

(61)

Thus, the Fourier transform $\tilde{\mathcal{V}}$ of $\mathcal{V}$ is diagonal in the k indices and contains $2\times 2$ blocks

$$\tilde{V}\left(k\right)=\left[\begin{array}{cc}{\tilde{V}}_{\downarrow \downarrow }\left(k\right)& {\tilde{V}}_{\downarrow \uparrow }\left(k\right)\\ {\tilde{V}}_{\uparrow \downarrow }\left(k\right)& {\tilde{V}}_{\uparrow \uparrow }\left(k\right)\end{array}\right].$$

(62)

The symmetries (50) of the screened Coulomb potential on the DNA lattice are reflected in the Fourier transforms (60) as well:

$$\begin{array}{cc}\hfill {\tilde{V}}_{b_1{b}_2}\left(k\right)& ={\tilde{V}}_{b_2{b}_1}{\left(k\right)}^{*}={\tilde{V}}_{b_2{b}_1}(-k)\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill {\tilde{V}}_{\uparrow \uparrow }\left(k\right)& ={\tilde{V}}_{\downarrow \downarrow }\left(k\right).\hfill \end{array}$$

(64)

Also, because $V_{\uparrow \uparrow }\left(\mathsf{\ell }\right)$ is an even function of ℓ, ${\tilde{V}}_{\uparrow \uparrow }\left(k\right)$ is real. The Fourier transforms (60) 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 as a function of the lattice index ℓ, resulting in the “wiggles” in Figure 12b 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 Equations (63) and (64) obeyed by the Fourier transforms allows us to simplify the $2\times 2$ blocks $\tilde{V}\left(k\right)$ in Equation (62) to

$$\tilde{V}\left(k\right)=\left[\begin{array}{cc}{\tilde{V}}_{\uparrow \uparrow }\left(k\right)& {\tilde{V}}_{\uparrow \downarrow }^{*}\left(k\right)\\ {\tilde{V}}_{\uparrow \downarrow }\left(k\right)& {\tilde{V}}_{\uparrow \uparrow }\left(k\right)\end{array}\right].$$

(65)

This matrix is easily diagonalized, yielding the real eigenvalues

$${\lambda }_{k,\pm }={\tilde{V}}_{\uparrow \uparrow }\left(k\right)\pm \left|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right|,$$

(66)

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

$$\xi \left(k\right)=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ \frac{{\tilde{V}}_{\uparrow \downarrow }\left(k\right)}{|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)|}& -\frac{{\tilde{V}}_{\uparrow \downarrow }\left(k\right)}{|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)|}\end{array}\right].$$

(67)

This diagonalizes the $2\times 2$ matrix $\tilde{V}\left(k\right)$ according to

$$\left[\begin{array}{cc}{\lambda }_{+}\left(k\right)& 0\\ 0& {\lambda }_{-}\left(k\right)\end{array}\right]={\xi }^{†}\left(k\right)\tilde{V}\left(k\right)\xi \left(k\right).$$

(68)

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

$${\mathcal{M}}_{(\mathsf{\ell },b),(k,\sigma )}=\frac{{\xi }_{b,\sigma }\left(k\right)}{\sqrt{2\mathcal{N}+1}}{e}^{ik\mathsf{\ell }},$$

(69)

which completely diagonalizes $\mathcal{V}$ such that

$$\Lambda ={\mathcal{M}}^{†}\mathcal{V}\mathcal{M}.$$

(70)

Here, $\Lambda$ is a diagonal matrix, with each element along the diagonal given by:

$${\Lambda }_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}={\delta }_{k_1,k_2}{\delta }_{{\sigma }_1,{\sigma }_2}{\lambda }_{k_1{\sigma }_1},$$

(71)

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 per site 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

$${\mathcal{W}}_{(\mathsf{\ell },b),(k,\sigma )}=\frac{1}{\sqrt{2}}\left\{{\mathcal{M}}_{(\mathsf{\ell },b),(k,\sigma )}+{\mathcal{M}}_{(\mathsf{\ell },b),(-k,\sigma )}\right\},k>0,$$

(72)

and

$${\mathcal{W}}_{(\mathsf{\ell },b),(k,\sigma )}=i\frac{1}{\sqrt{2}}\left\{{\mathcal{M}}_{(\mathsf{\ell },b),(k,\sigma )}-{\mathcal{M}}_{(\mathsf{\ell },b),(-k,\sigma )}\right\},k<0.$$

(73)

Since ${\xi }_{b\sigma }(-k)={\xi }_{b\sigma }^{*}\left(k\right)$, these can be written in terms of real and imaginary parts as

$${\mathcal{W}}_{(\mathsf{\ell },b),(k,\sigma )}=\frac{\sqrt{2}}{\sqrt{(2\mathcal{N}+1)}}\mathrm{Re}\left\{e^{ik\mathsf{\ell }}{\xi }_{b\sigma }\left(k\right)\right\},k>0.$$

(74)

and

$${\mathcal{W}}_{(\mathsf{\ell },b),(k,\sigma )}=-\frac{\sqrt{2}}{\sqrt{(2\mathcal{N}+1)}}\mathrm{Im}\left\{e^{ik\mathsf{\ell }}{\xi }_{b\sigma }\left(k\right)\right\},k<0.$$

(75)

when $k=0$, the orthogonal transformation is the same as the unitary transformation, that is

$${\mathcal{W}}_{(\mathsf{\ell },b),(0,\sigma )}={\mathcal{M}}_{(\mathsf{\ell },b),(0,\sigma )}=\frac{{\xi }_{b,\sigma }\left(0\right)}{\sqrt{2\mathcal{N}+1}},$$

(76)

where

$$\xi \left(0\right)=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right].$$

(77)

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

$${\mathcal{Z}}_{\left(k\sigma \right),\left(k^{\prime }{\sigma }^{\prime }\right)}={\delta }_{\sigma ,{\sigma }^{\prime }}\left\{\begin{array}{cc}\frac{i}{\sqrt{2}}({\delta }_{k,k^{\prime }}-{\delta }_{k,-k^{\prime }})& k<0\\ {\delta }_{k,k^{\prime }}& k=0\\ \frac{1}{\sqrt{2}}({\delta }_{k,k^{\prime }}+{\delta }_{k,-k^{\prime }})& k>0\end{array}\right..$$

(78)

Then, the orthogonal transformation is

$$\mathcal{W}=\mathcal{M}\mathcal{Z}.$$

(79)

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.

5. Partition Function for Interacting DNA Chains

For the interacting lattice gas, the partition function is similar to the noninteracting one in Equations (5) and (8). However, here the full Hamiltonian $\mathcal{H}$ from Equation (39) is used, and then the partition function is given by

$$Z_G=\sum _{\mathrm{configurations}}{e}^{-\beta \left(\mathcal{H}-{\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }{\mu }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}\right)}.$$

(80)

The full Hamiltonian, including the interaction term from Equation (42) written in matrix form, is

$$\mathcal{H}=\sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}+\frac{1}{2}{\rho }^T\mathcal{V}\rho .$$

(81)

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 $\mathcal{V}$ in the interaction term of the Hamiltonian, which gives

$$\mathcal{H}=\sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}+\frac{1}{2}\left({\rho }^T\mathcal{W}\right)\left({\mathcal{W}}^T\mathcal{V}\mathcal{W}\right)\left({\mathcal{W}}^T\rho \right),$$

(82)

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

$$\begin{array}{ccc}\hfill \mathcal{H}& =& \sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}+\frac{1}{2}{\tilde{\rho }}^2\Lambda \hfill \\ & =& \sum _{\mathsf{\ell },b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\mathsf{\ell },b}^{\left(\gamma \right)}+\frac{1}{2}\sum _{k,\sigma }{\left({\tilde{\rho }}_{k,\sigma }\right)}^2{\lambda }_{k,\sigma },\hfill \end{array}$$

(83)

where $\tilde{\rho }={\mathcal{W}}^T\rho$ is the transformed charge per site. 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

$$Z_G=\sum _{\mathrm{configurations}}{e}^{-\beta {\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{\mathsf{\ell },b}^{\left(\gamma \right)}}{e}^{-\frac{\beta }{2}{\sum }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }^2{\lambda }_{k,\sigma }}.$$

(84)

We can write the exponential of the sum in the interaction term as a product of exponentials as

$$Z_G=\sum _{\mathrm{configurations}}{e}^{-\beta {\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{\mathsf{\ell },b}^{\left(\gamma \right)}}\prod _{k,\sigma }{e}^{-\frac{\beta }{2}{\tilde{\rho }}_{k,\sigma }^2{\lambda }_{k,\sigma }}.$$

(85)

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),

$$\langle n^{\left(\gamma \right)}\rangle =\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta Z_{\mathrm{G}}}\frac{\partial Z_{\mathrm{G}}}{\partial {\mu }_{\gamma }},$$

(86)

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_{\mathsf{\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_{\mathsf{\ell },b}^{\left(\gamma \right)}$ only appeared in its exponent. This is because only linear factors of $n_{\mathsf{\ell },b}^{\left(\gamma \right)}$ were in the exponent in the partition function. While it is still true that $n_{\mathsf{\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_{\mathsf{\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_{\mathsf{\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 achieve 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

$$\int \sqrt{\frac{\beta \lambda }{-2\pi }}d\tilde{\Delta }{e}^{\frac{\beta }{2}\lambda {\tilde{\Delta }}^2-\beta \lambda \tilde{\Delta }\tilde{\rho }}{e}^{\frac{\beta }{2}\lambda {\tilde{\rho }}^2}=1,$$

(87)

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 performed 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 [50]. 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 Equation (87), 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

$$\begin{array}{ccc}\hfill Z_G& =& \sum _{\mathrm{configurations}}\prod _{k,\sigma }\int \sqrt{\frac{\beta {\lambda }_{k,\sigma }}{-2\pi }}d{\tilde{\Delta }}_{k,\sigma }{e}^{\frac{\beta }{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2-\beta {\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }{\tilde{\rho }}_{k,\sigma }}{e}^{\frac{\beta }{2}{\lambda }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }^2}\times \hfill \\ & & \times e^{-\beta {\sum }_{\mathsf{\ell },b}{\sum }_{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{\mathsf{\ell },b}^{\left(\gamma \right)}}{e}^{-\frac{\beta }{2}{\tilde{\rho }}_{k,\sigma }^2{\lambda }_{k,\sigma }}.\hfill \end{array}$$

(88)

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 perform 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

$$\begin{array}{ccc}\hfill Z_G& =& \int \left(\prod _{k,\sigma }\sqrt{\frac{\beta {\lambda }_{k,\sigma }}{-2\pi }}d{\tilde{\Delta }}_{k,\sigma }\right)\sum _{\mathrm{configurations}}exp\left\{-\beta \sum _{\mathsf{\ell },b}\sum _{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{\mathsf{\ell },b}^{\left(\gamma \right)}\right.+\hfill \\ & & \left.+\sum _{k,\sigma }\left[\frac{\beta }{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2-\beta {\lambda }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }\right]\right\}.\hfill \end{array}$$

(89)

It will be useful to define an action $\mathcal{S}$ such that the partition function can be written in the form

$$Z_G=\int \mathcal{D}\left[\tilde{\Delta }\right]\sum _{\mathrm{configurations}}{e}^{-\mathcal{S}\left[\tilde{\Delta }\right]},$$

(90)

where

$$\mathcal{S}\left[\tilde{\Delta }\right]=\beta \left\{\sum _{\mathsf{\ell },b}\sum _{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{\mathsf{\ell },b}^{\left(\gamma \right)}+\sum _{k_{,}\sigma }\left[-\frac{\beta }{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2+\beta {\lambda }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }\right]\right\}$$

(91)

and

$$\mathcal{D}\left[\tilde{\Delta }\right]=\left(\prod _{k,\sigma }\sqrt{\frac{\beta {\lambda }_{k,\sigma }}{-2\pi }}d{\tilde{\Delta }}_{k,\sigma }\right).$$

(92)

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 performed 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 [51].

We see that by diagonalizing first, our auxiliary fields are in the diagonal basis of the potential. If we transform the terms containing ${\epsilon }_{\gamma }$ and ${\mu }_{\gamma }$ to this basis, they will no longer be diagonal. Therefore, we will leave those terms in the position basis. It will also be convenient to have the term containing ${\sum }_{k_{,}\sigma }{\lambda }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }$ in the position basis also. Therefore, the expression for the transformed charge per site is

$$\begin{array}{ccc}\hfill {\tilde{\rho }}_{k,\sigma }& =& \sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\sum _{\gamma }{\rho }_{\mathsf{\ell },b}^{\left(\gamma \right)}\hfill \\ & =& \sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\sum _{\gamma }{q}_{\gamma }{n}_{(\mathsf{\ell },b)}^{\left(\gamma \right)},\hfill \end{array}$$

(93)

in terms of the quantities in real space.

Substituting this into the action, we have

$$\mathcal{S}\left[\tilde{\Delta }\right]=\beta \left\{\sum _{\mathsf{\ell },b}\sum _{\gamma }({\epsilon }_{\gamma }-{\mu }_{\gamma })n_{(\mathsf{\ell },b)}^{\left(\gamma \right)}-\sum _{k_{,}\sigma }\left[\frac{1}{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2+{\tilde{\Delta }}_{k,\sigma }{\lambda }_{k,\sigma }\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\sum _{\gamma }{q}_{\gamma }{n}_{(\mathsf{\ell },b)}^{\left(\gamma \right)}\right]\right\}.$$

(94)

It is convenient to identify the last term in this expression as a self-energy of species $\gamma$,

$${\tilde{\Sigma }}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]=\left(\sum _{k,\sigma }{\tilde{\Delta }}_{k,\sigma }{\lambda }_{k,\sigma }{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\right)q_{\gamma }.$$

(95)

By transforming ${\tilde{\Delta }}_{k,\sigma }$ to the position basis as

$${\tilde{\Delta }}_{k,\sigma }=\sum _{{\mathsf{\ell }}^{\prime },b^{\prime }}{\Delta }_{{\mathsf{\ell }}^{\prime },b^{\prime }}{\mathcal{W}}_{({\mathsf{\ell }}^{\prime },b^{\prime }),(k,\sigma )},$$

(96)

this self-energy can also be written in the position basis as

$$\begin{array}{ccc}\hfill {\Sigma }_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]& =& {\tilde{\Sigma }}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]\hfill \\ & =& q_{\gamma }\sum _{{\mathsf{\ell }}^{\prime },b^{\prime }}{\Delta }_{{\mathsf{\ell }}^{\prime },b^{\prime }}\left(\sum _{k,\sigma }{\mathcal{W}}_{({\mathsf{\ell }}^{\prime },b^{\prime }),(k,\sigma )}{\lambda }_{k,\sigma }{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\right)\hfill \\ & =& q_{\gamma }\sum _{{\mathsf{\ell }}^{\prime },b^{\prime }}{\Delta }_{{\mathsf{\ell }}^{\prime },b^{\prime }}{\mathcal{V}}_{({\mathsf{\ell }}^{\prime },b^{\prime }),(\mathsf{\ell },b)}.\hfill \end{array}$$

(97)

This will be a useful form to use when we discuss the saddle-point, or mean-field, approximation.

Keeping the first term in the position basis and the second term in the diagonal basis, we now write the action in a compact form, combining the self-energy term with the first term in the action to obtain

$$\mathcal{S}\left[\tilde{\Delta }\right]=\beta \left\{\sum _{\mathsf{\ell },b}\sum _{\gamma }({\epsilon }_{\gamma }+{\Sigma }_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]-{\mu }_{\gamma })n_{(\mathsf{\ell },b)}^{\left(\gamma \right)}-\sum _{k_{,}\sigma }\left[\frac{1}{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2\right]\right\}.$$

(98)

Analogous to the noninteracting case, we define the activity as

$${\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]=a_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]=e^{-\beta ({\epsilon }_{\gamma }+{\Sigma }_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]-{\mu }_{\gamma })}.$$

(99)

With this definition, the partition function becomes

$$Z_G=\int \mathcal{D}\left[\tilde{\Delta }\right]\sum _{\mathrm{configurations}}\left(\prod _{\mathsf{\ell },b}\prod _{\gamma }{\left\{a_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]\right\}}^{n_{(\mathsf{\ell },b)}^{\left(\gamma \right)}}\right)e^{\beta {\sum }_{k_{,}\sigma }\left[\frac{1}{2}{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2\right]}.$$

(100)

The trace over configurations is now performed over each site separately, exactly as we did for the noninteracting case, where $n_{(\mathsf{\ell },b)}^{\left(\gamma \right)}=1$ for one and only one of $\gamma =\Vert ,\perp$, or v, and zero otherwise. Thus, performing the trace gives

$$\sum _{\begin{array}{c}\mathrm{single-site}\\ \mathrm{configurations}\end{array}}\prod _{\gamma }{\left\{{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]\right\}}^{n_{(\mathsf{\ell },b)}^{\left(\gamma \right)}}=\sum _{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right],$$

(101)

and the grand partition function becomes

$$Z_G=\int \mathcal{D}\left[\tilde{\Delta }\right]e^{-{\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]},$$

(102)

where the effective action is

$${\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]=-\frac{\beta }{2}\sum _{k_{,}\sigma }{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2-\sum _{\mathsf{\ell },b}ln\left(\sum _{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]\right).$$

(103)

This is the general result, which is in principle exact. To understand what the Hubbard–Stratonovich transformation has accomplished for us, recall that the interaction energy $U_{\mathrm{int}}$ posed two difficulties: the additional configuration dependence and the coupling between the sites. By introducing auxiliary fields through the Hubbard–Stratonovich transformation (87), we managed to separate these two complications so that one factor $exp\left[-\frac{\beta }{2}{\sum }_{k_{,}\sigma }{\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }^2\right]$ contains interactions between field fluctuations but no configuration dependence, and another factor $exp\left[\beta {\sum }_{k,\sigma }{\lambda }_{k,\sigma }{\tilde{\rho }}_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }\right]$ is configuration-dependent (through ${\tilde{\rho }}_{k,\sigma }$) but decoupled. This allowed us to define a modified activity (99) that incorporated all the configuration dependence, making it possible to evaluate the sum over configurations directly, as in the noninteracting case.

What we have achieved in using the Hubbard–Stratonovich transformation is to replace the interaction part of the partition function with its functional integral representation. From the form of Equation (100), we see that the largest contribution to the partition function comes from the values of ${\tilde{\Delta }}_{k,\sigma }^2$ for which the effective action ${\mathcal{S}}_{\mathrm{eff}}$ is stationary, i.e., at a saddle point. Thus, a Taylor series expansion of the effective action (103) about the saddle point will yield an order-by-order approximation to the exact partition function (102).

6. Expansion in Powers of the Auxiliary Field

Since it is not possible to evaluate the integral in the partition function of Equation (102) exactly, we will approximate it by expanding the action in Equation (103) about a mean field ${\Delta }_c$, and include fluctuations to the lowest order about this mean field. This mean field is determined by the saddle point of the integrand in Equation (102), which is the largest contribution to the integral.

We can write the effective action, then, as an expansion about this uniform saddle point. As an aid to keeping track of the order of the expansions, we add a multiplicative factor m in front of the effective action:

$$\begin{array}{ccc}\hfill m{\mathcal{S}}_{\mathrm{eff}}\left(\tilde{\Delta }={\tilde{\Delta }}_c+\delta \tilde{\Delta }/\sqrt{m}\right)& =& m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]+{\left.\sum _{k,\sigma }\frac{\delta {\tilde{\Delta }}_{k,\sigma }}{\sqrt{m}}m\frac{\partial {\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}\right|}_{{\tilde{\Delta }}_c}+\hfill \\ & +& {\left.\frac{1}{2}\sum _{k_1,{\sigma }_1,k_2,{\sigma }_2}\frac{\delta {\tilde{\Delta }}_{k_1,{\sigma }_1}}{\sqrt{m}}m\frac{{\partial }^2{\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k_1,{\sigma }_1}\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}\right|}_{{\tilde{\Delta }}_c}\frac{\delta {\tilde{\Delta }}_{k_2,{\sigma }_2}}{\sqrt{m}}+\dots .\hfill \end{array}$$

(104)

Here, the expansion parameter $\delta \tilde{\Delta }$ is scaled by $m^{-1/2}$ by writing $\delta {\tilde{\Delta }}_{k,\sigma }\equiv \sqrt{m}[{\tilde{\Delta }}_{k,\sigma }-{\tilde{\Delta }}_{k,\sigma ,c}]$. The factors of m cancel in the quadratic term, so that it is of order $m^0$.

In evaluating the integral in the partition function, the largest contribution comes from the region near a saddle point, and the terms beyond the saddle point give corrections to the integral. Physically, the saddle-point solution corresponds to a mean-field approximation, and we begin with that. The saddle point is found by setting the first derivative in the expansion to zero, that is

$${\left.\frac{\partial {\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}\right|}_{{\tilde{\Delta }}_c}=0.$$

(105)

The second derivatives evaluated at the saddle point are proportional to the components of a matrix ${\tilde{D}}^{-1}$, which is the inverse propagator for the fluctuations in the auxiliary Hartree field $\Delta$ that we introduced to decouple the interparticle interactions. Because we want to absorb the normalization factor in the grand differential $\mathcal{D}\left[\tilde{\Delta }\right]$, we have defined each of the elements of ${\tilde{D}}^{-1}$ by the second derivative of ${\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]$ divided by $-\beta {\lambda }_{k,\sigma }$ as

$${\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}=-\left(\frac{1}{\beta {\lambda }_{k,\sigma }}\right){\left.\frac{{\partial }^2{\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k_1,{\sigma }_1}\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}\right|}_{{\tilde{\Delta }}_c}.$$

(106)

We will show in Section 8 that, for a spatially uniform saddle point, this matrix is diagonal, that is, ${\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]$ is diagonal in the $\delta {\tilde{\Delta }}_{k_1,{\sigma }_1}$’s, that is,

$${\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}={\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_1,{\sigma }_1)}{\delta }_{k_1,k_2}{\delta }_{{\sigma }_1,{\sigma }_2}.$$

(107)

Although the Gaussian integral can still be performed if ${\tilde{D}}^{-1}$ is not diagonal, the argument becomes much easier to follow if we assume at this point that ${\tilde{D}}^{-1}$ is diagonal, because then the Gaussian integral can be performed straightforwardly.

With this assumption that ${\tilde{D}}^{-1}$ is diagonal, the expansion of the effective action reduces to

$$m{\mathcal{S}}_{\mathrm{eff}}\left(\tilde{\Delta }={\tilde{\Delta }}_c+\delta \tilde{\Delta }/\sqrt{m}\right)\simeq m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]-\frac{1}{2}\sum _{k,\sigma }\beta {\lambda }_{k,\sigma }{\left({\tilde{D}}^{-1}\right)}_{(k,\sigma ),(k,\sigma )}{\left(\delta {\tilde{\Delta }}_{k,\sigma }\right)}^2.$$

(108)

The grand canonical partition function can then be written as

$$Z_G\simeq e^{-m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}\int \left(\prod _{k,\sigma }\sqrt{\frac{\beta {\lambda }_{k,\sigma }}{-2\pi }}d\left(\delta {\tilde{\Delta }}_{k,\sigma }\right)\right)e^{\frac{1}{2}\beta {\sum }_{k,\sigma }{\lambda }_{k,\sigma }{\left({\tilde{D}}^{-1}\right)}_{(k,\sigma ),(k,\sigma )}{\left(\delta {\tilde{\Delta }}_{k,\sigma }\right)}^2}.$$

(109)

This can be written as a product of Gaussian integrals as

$$Z_G\simeq e^{-m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}\prod _{k,\sigma }\int \left(\sqrt{\frac{\beta {\lambda }_{k,\sigma }}{-2\pi }}d\left(\delta {\tilde{\Delta }}_{k,\sigma }\right)\right)e^{\frac{1}{2}\beta {\lambda }_{k,\sigma }{\left({\tilde{D}}^{-1}\right)}_{(k,\sigma ),(k,\sigma )}{\left(\delta {\tilde{\Delta }}_{k,\sigma }\right)}^2},$$

(110)

and each Gaussian integral can be evaluated using the same change in variable and path of integration as was used in the Hubbard–Stratonovich transformation in Equation (87), that originally was used to define the auxiliary fields, and is discussed in Appendix A. Each Gaussian integral produces a result of the same form, regardless of whether ${\left({\tilde{D}}^{-1}\right)}_{(k,\sigma ),(k,\sigma )}$ is positive or negative. This quantity can never be zero, since we would not have defined an auxiliary field for that case, which is the same as the case in which an eigenvalue of the potential is zero. Therefore, the integral can always be evaluated, and the partition function can be written as

$$\begin{array}{ccc}\hfill Z_G& \simeq & e^{-m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}\prod _{k,\sigma }\frac{1}{\sqrt{{\left({\tilde{D}}^{-1}\right)}_{(k,\sigma ),(k,\sigma )}}}.\hfill \end{array}$$

(111)

Note that all the factors of $\beta {\lambda }_{k,\sigma }$ have canceled out due to our choices of the form of the effective action ${\mathcal{S}}_{\mathrm{eff}}$ and the inverse Hartree-field-fluctuation propagator ${\tilde{D}}^{-1}$. Choosing these quantities carefully allows for the correct normalization of the integral [51]. The product can now be brought inside the square root, and since ${\tilde{D}}^{-1}$ is diagonal, this gives the determinant of the matrix inside the square root, and so the partition function assumes the form

$$\begin{array}{ccc}\hfill Z_G& \simeq & e^{-m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}\frac{1}{\sqrt{det\left({\tilde{D}}^{-1}\right)}},\hfill \end{array}$$

(112)

and the partition function can be written as a single exponential,

$$\begin{array}{ccc}\hfill Z_G& \simeq & e^{-\{m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]+\frac{1}{2}lndet\left({\tilde{D}}^{-1}\right)\}}.\hfill \end{array}$$

(113)

The grand canonical potential is then

$${\Omega }_G=-\frac{1}{\beta }ln{Z}_G\simeq \frac{1}{\beta }\{m{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]+\frac{1}{2}lndet\left({\tilde{D}}^{-1}\right)\}.$$

(114)

If we formally write the expansion in terms of orders in m, we have

$${\Omega }_G\simeq m{\Omega }_{m=1}+{\Omega }_{m=0}+\frac{1}{\sqrt{m}}{\Omega }_{m=-\frac{1}{2}}+\frac{1}{m}{\Omega }_{m=-1}+\cdots .$$

(115)

The leading term in the expansion is, therefore, given by the saddle-point value

$$m{\Omega }_{m=1}=m{\Omega }_c=\frac{m}{\beta }{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right],$$

(116)

and what is known as the “one-loop” correction term is

$${\Omega }_{m=0}=\frac{1}{2}lndet\left({\tilde{D}}^{-1}\right)=\frac{1}{2}\mathrm{Tr}ln{\tilde{D}}^{-1}.$$

(117)

7. The Saddle-Point, or Mean-Field, Approximation

The first task in this section is to find an equation that determines the value of the auxiliary field at the saddle-point, or mean-field, value of the effective action. In order to find this saddle point, we must set the first derivative of the effective action in Equation (103) with respect to the auxiliary field ${\tilde{\Delta }}_{k,\sigma }$ to zero. That derivative is given by

$$\frac{\partial {\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}=-\beta {\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }-\beta \sum _{\mathsf{\ell },b}\frac{1}{\left({\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]\right)}\sum _{\gamma }\frac{\partial {\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}.$$

(118)

The activity $a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]$ for species $\gamma$ is given by Equation (99), and its derivative is

$$\frac{\partial {\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}=-\beta {\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]\frac{\partial {\tilde{\Sigma }}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }},$$

(119)

where from Equation (95), the derivative of the self-energy is

$$\frac{\partial {\tilde{\Sigma }}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}=q_{\gamma }{\lambda }_{k,\sigma }{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T.$$

(120)

Substituting this into the derivative of the action, we have

$$\frac{\partial {\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k,\sigma }}=-\beta {\lambda }_{k,\sigma }{\tilde{\Delta }}_{k,\sigma }+\beta {\lambda }_{k,\sigma }\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\frac{{\sum }_{\gamma }{q}_{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}.$$

(121)

Setting this first derivative to zero, we obtain an equation whose solution gives the auxiliary field in the saddle-point approximation,

$${\tilde{\Delta }}_{(k,\sigma )}=\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T\left(\frac{{\sum }_{\gamma }{q}_{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}\right).$$

(122)

We see that on the right-hand side of the equation is the transformation to the diagonal basis of the potential, as given in Equation (96). We could, therefore, write the auxiliary field in the position basis as

$${\Delta }_{(\mathsf{\ell },b)}=\frac{{\sum }_{\gamma }{q}_{\gamma }{a}_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]}{{\sum }_{{\gamma }^{\prime }}{a}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}[\Delta ]},$$

(123)

where $a_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]$ is the activity with the self-energy in the diagonal basis ${\tilde{\Sigma }}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]$, which by Equation (97) has been replaced by the self-energy in the position basis, ${\Sigma }_{\mathsf{\ell },b}^{\left(\gamma \right)}[\Delta ]$. Either of these two equations, Equation (122) or (123), may be considered the equation giving the saddle-point, or mean-field, value of the auxiliary field, and we will refer to this as the mean-field equation. In this expression, $a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]$ is a function of all the ${\Delta }_{(\mathsf{\ell },b),c}$ elements, and so, for $\mathcal{N}=30$, which we have used in our calculations, there are ${\mathcal{N}}_{\mathrm{sites}}=2(2\mathcal{N}+1)=122$ coupled nonlinear equations. Therefore, for simplicity, we make the assumption that the mean field is spatially uniform, which means that

$${\Delta }_{(\mathsf{\ell },b),c}={\Delta }_c,$$

(124)

where ${\Delta }_c$ is constant. This is still a nonlinear equation in ${\Delta }_c$, and so this must be solved numerically with an iterative approach.

Before solving the mean-field equation, it is useful to have a physical interpretation of the mean auxiliary field ${\Delta }_c$. The first step in this direction is the second task of this section, which is to find the mean site occupancies of species $\gamma$ at the mean-field level. This is achieved as in the noninteracting case, by taking the derivative of the partition function with respect to ${\mu }_{\gamma }$. At this mean-field level, the ${\mu }_{\gamma }$ dependence is contained in ${\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]$ in Equation (113), neglecting the term containing the inverse propagator $D^{-1}$. The expression for the mean occupancy can be obtained using Equation (86) as

$${\langle n^{\left(\gamma \right)}\rangle }_c=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta }\frac{\partial }{\partial {\mu }_{\gamma }}ln\left[e^{-{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}\right].$$

(125)

Taking the derivative of the exponential, we can cancel out the partition function in the denominator and we are left with the derivative of the effective action as

$${\langle n^{\left(\gamma \right)}\rangle }_c=-\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta }\frac{\partial }{\partial {\mu }_{\gamma }}{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right].$$

(126)

Substituting the effective action from Equation (103), we have

$${\langle n^{\left(\gamma \right)}\rangle }_c=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta }\sum _{\mathsf{\ell },b}\frac{1}{\left({\sum }_{{\gamma }^{\prime }}{a}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[{\Delta }_c\right]\right)}\frac{\partial a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]}{\partial {\mu }_{\gamma }}.$$

(127)

The derivative of the activity is

$$\frac{\partial a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]}{\partial {\mu }_{\gamma }}=\beta a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right],$$

(128)

and so the mean occupancy becomes

$${\langle n^{\left(\gamma \right)}\rangle }_c=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}}\sum _{\mathsf{\ell },b}\frac{a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]}{\left({\sum }_{{\gamma }^{\prime }}{a}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[{\Delta }_c\right]\right)}.$$

(129)

The third task of this section is to relate the mean charge per site to the mean occupancy. The charge per site has the same form as in Equation (22) for the noninteracting model, so that at the mean-field level we obtain

$${\rho }_c={\langle \rho \rangle }_c=\sum _{\gamma }{q}_{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_c=\frac{{\sum }_{\gamma }{q}_{\gamma }{a}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]}{{\sum }_{{\gamma }^{\prime }}{a}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[{\Delta }_c\right]},$$

(130)

which, like the mean site occupancy, is independent of the site index because of our choice of spatially uniform mean field. This is equivalent to the equation for the saddle-point value of the auxiliary field because the right-hand side of the equation shows that

$${\rho }_c={\Delta }_c,$$

(131)

on comparison with Equation (123).

In the mean-field, or saddle-point, equation for the spatially uniform saddle point, the auxiliary field is the same for all lattice sites, and the self-energy becomes spatially uniform also and can be written as

$${\Sigma }_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]\equiv {\Sigma }_c^{\left(\gamma \right)}={\rho }_c{q}_{\gamma }{\mathbb{S}}_{\mathrm{lattice}}.$$

(132)

where ${\mathbb{S}}_{\mathrm{lattice}}$ is a lattice sum that is independent of the choice of lattice site $(\mathsf{\ell },b)$ for a long DNA double helix, and is given by

$${\mathbb{S}}_{\mathrm{lattice}}=\sum _{{\mathsf{\ell }}^{\prime },b^{\prime }}{\mathcal{V}}_{({\mathsf{\ell }}^{\prime },b^{\prime }),(\mathsf{\ell },b)}.$$

(133)

For the parameters used in the figures, $\beta {\mathbb{S}}_{\mathrm{lattice}}\approx 2.71$.

Consequently, at the mean-field level, the self-energy simply represents the electrostatic energy of charge $q_{\gamma }$ interacting with the rest of the DNA lattice. This term plays the role of the self-energy $\Sigma$ in many-body quantum mechanics [28,50,51]. Since the self-energy is independent of the lattice site, so is the activity, which can be written as

$$a_c^{\left(\gamma \right)}=a_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[{\Delta }_c\right]=e^{-\beta ({\epsilon }_{\gamma }-{\mu }_{\gamma }+{\Sigma }_c^{\left(\gamma \right)})}.$$

(134)

The simplifications made possible by the assumption of a spatially uniform saddle point allow us to rewrite the mean-field equation given by Equation (123) in an explicit form as

$${\rho }_c=\frac{{\sum }_{\gamma }{q}_{\gamma }{e}^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }+q_{\gamma }{\rho }_c{\mathbb{S}}_{\mathrm{lattice}}\right)}}{{\sum }_{\gamma }{e}^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }+q_{\gamma }{\rho }_c{\mathbb{S}}_{\mathrm{lattice}}\right)}}.$$

(135)

Again, recall that ${\epsilon }_v-{\mu }_v=0$, identically, $q_{\Vert }=0$, and ${\mu }_{\Vert }={\mu }_{\perp }={\mu }_{\mathrm{dimer}}$. Therefore, when ${\rho }_c=0$, we can see from Equation (135) that $q_{\perp }{e}^{-\beta ({\epsilon }_{\perp }-{\mu }_{\mathrm{dimer}})}+q_v=0$, and since $q_{\perp }=+1$, and $q_v=-1$, it follows that ${\epsilon }_{\perp }={\mu }_{\mathrm{dimer}}$. Therefore, this is the point at which the charge per site ${\rho }_c$ goes to zero. When ${\epsilon }_{\Vert }=2{\epsilon }_{\perp }$, the charge per site goes to zero at ${\epsilon }_{\Vert }=2{\mu }_{\mathrm{dimer}}$, which is the same place it went to zero for the noninteracting model, as shown in Figure 6.

Returning to the mean site occupancy in Equation (129), the quantity inside the sum over ℓ and b is constant and the sum yields ${\mathcal{N}}_{\mathrm{sites}}$, and so the average occupancy per species becomes

$${\langle n^{\left(\gamma \right)}\rangle }_c=\frac{a_c^{\left(\gamma \right)}}{\left({\sum }_{{\gamma }^{\prime }}{a}_c^{\left({\gamma }^{\prime }\right)}\right)}.$$

(136)

Using the mean-field value ${\rho }_c={\Delta }_c$ of the charge per site, calculated from Equation (135), we have calculated the average occupation numbers $\langle n^{\left(\gamma \right)}\rangle$ via Equation (136) using Equation (134) for the activities. In Figure 15, we have plotted $\delta {\langle n^{\left(\gamma \right)}\rangle }_c={\langle n^{\left(\gamma \right)}\rangle }_c-{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}$, which are the differences in the mean-field occupancies from the occupancies in the noninteracting case shown in Figure 5. The corresponding change in the charge per site $\delta {\rho }_c={\rho }_c-{\rho }_{\mathrm{NI}}$ is shown as the solid blue curve in Figure 16.

The shifts in occupancy plotted in Figure 15 can be primarily understood as a consequence of the electromagnetic “self-energy” ${\Sigma }_c^{\left(\gamma \right)}={\rho }_c{q}_{\gamma }{\mathbb{S}}_{\mathrm{lattice}}$ of each species $\gamma$ interacting with the total charge ${\rho }_c$ of the lattice. Immediately, this implies that the parallel dimers are hardly modified at all, since they are electrically neutral and not directly affected by the total charge ${\rho }_c$ of the lattice. For the vacancies and perpendicular dimers, which are charged, the mean-field corrections plotted in Figure 15 and Figure 16 act to reduce the overall charge of the lattice. Below $\beta {\epsilon }_{\Vert }\approx 2$, the total charge of the lattice is positive, so the mean-field corrections penalize the positively charged perpendicular dimers and enhance the negatively charged vacancies. Above $\beta {\epsilon }_{\Vert }\approx 2$, the total charge of the lattice is negative, and this effect is reversed. Thus, the mean-field corrections tend to reduce (but not eliminate) the charge inversion in the weak binding region, $\beta {\epsilon }_{\Vert }>2$, seen in Figure 5 for the non-interacting model. Note that the biggest effect of this electrostatic self-energy is to heavily penalize the “naked lattice” of negatively charged vacancies, which was the favored ground state at large positive ${\epsilon }_{\Vert }$ for the noninteracting Hamiltonian.

The fourth task of this section is to determine the variance in the charge per site from the mean field, similar to Equation (24) for the noninteracting model, which gives a measure of the importance of spatial fluctuations in the charge per site, and can be written in the form

$${\sigma }_c^2={\langle {\rho }^2\rangle }_c-{\rho }_c^2,$$

(137)

where the average ${\langle {\rho }^2\rangle }_c$ of the square of the charge per site, similar to Equation (29) for the noninteracting model, is

$${\langle {\rho }^2\rangle }_c=\sum _{\gamma }{q}_{\gamma }^2{\langle n^{\left(\gamma \right)}\rangle }_c=\frac{{\sum }_{\gamma }{q}_{\gamma }^2{a}_c^{\left(\gamma \right)}}{{\sum }_{{\gamma }^{\prime }}{a}_c^{\left({\gamma }^{\prime }\right)}}.$$

(138)

The variance ${\sigma }_c^2$ will be useful for knowing the size of the charge fluctuations about the mean field, and this quantity will appear in the next level of approximation. The difference in the standard deviation $\delta {\sigma }_c={\sigma }_c-{\sigma }_{\mathrm{NI}}$ is plotted as the red dashed curve in Figure 16, where ${\sigma }_{\mathrm{NI}}$ for the noninteracting model is plotted as the red dashed curve in Figure 6.

Using a spatially uniform saddle point, which yields a charge per site ${\rho }_c$ that is equal at the mean-field level to the auxiliary field ${\Delta }_c$ in real space, it is interesting to write down the auxiliary field in the diagonal basis, given by

$${\tilde{\Delta }}_{(k,\sigma ),c}=\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T{\left.\langle {\rho }_{\mathsf{\ell },b}\rangle \right|}_c={\rho }_c\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T.$$

(139)

From Equations (74) and (75), we see that when $k\ne 0$, ${\mathcal{W}}_{(k,\sigma ),(\mathsf{\ell },b)}^T$ has the ℓ dependence of the form $e^{\pm ik\mathsf{\ell }}$, which by Equation (59) sums to zero. The only surviving term is for $k=0$, which when summed over ℓ gives $2\mathcal{N}+1$, and we have from Equations (76) and (139)

$${\tilde{\Delta }}_{(k,\sigma ),c}={\tilde{\rho }}_{(k,\sigma ),c}={\delta }_{k,0}{\rho }_c\sqrt{2\mathcal{N}+1}\sum _b{\xi }_{b,\sigma }\left(0\right).$$

(140)

There are two different results corresponding to the two different eigenvalue labels $\sigma =+$ and $\sigma =-$ in Equation (66). Summing over b, these are

$$\begin{array}{ccc}\hfill {\tilde{\Delta }}_{(0,+),c}& =& {\tilde{\rho }}_{(0,+),c}={\rho }_c\sqrt{{\mathcal{N}}_{\mathrm{sites}}}\hfill \\ \hfill {\tilde{\Delta }}_{(0,-),c}& =& {\tilde{\rho }}_{(0,-),c}=0,\hfill \end{array}$$

(141)

where ${\mathcal{N}}_{\mathrm{sites}}=2(2\mathcal{N}+1)$. The interpretation of Equation (141) can be understood from Equation (77). The sum over b in Equation (140) is a sum over the element rows in Equation (77), which for $\sigma =+$ sums the mean charges per site on the two chains and for $\sigma =-$ subtracts the charges per site on the two chains.

The fifth task of this section is to calculate the entropy for the spatially uniform saddle point, which is our mean-field value of the entropy. This requires the effective action at this saddle point, which becomes

$${\mathcal{S}}_{\mathrm{eff},c}=-{\mathcal{N}}_{\mathrm{sites}}ln\left(\sum _{\gamma }{a}_c^{\left(\gamma \right)}\right)-\frac{{\rho }_c^2}{2}{\mathcal{N}}_{\mathrm{sites}}\beta {\mathbb{S}}_{\mathrm{lattice}}.$$

(142)

The grand canonical potential is then

$${\Omega }_{m=1}={\Omega }_c=\frac{1}{\beta }{\mathcal{S}}_{\mathrm{eff},\mathrm{c}}.$$

(143)

We can now compute the entropy from the derivative of the grand thermodynamic potential with respect to the inverse temperature $\beta$, as in Section 2. This derivative becomes much more complicated with the introduction of $U_{\mathrm{int}}$ because the screened Coulomb potential (43) depends on $\beta$ through the screening vector $q_s$ (44), where

$$\frac{\partial q_s}{\partial \beta }=\frac{q_s}{2\beta }.$$

(144)

This introduces a $\beta$-dependence into both terms of Equation (142), so that the mean-field entropy is given by the more complex expression

$${\overline{S}}_c=\frac{S_c}{k_B{\mathcal{N}}_{\mathrm{sites}}}=\frac{S_c^{\left(A\right)}+S_c^{\left(B\right)}}{k_B{\mathcal{N}}_{\mathrm{sites}}},$$

(145)

where

$${\overline{S}}_c^{\left(A\right)}=\frac{S_c^{\left(A\right)}}{k_B{\mathcal{N}}_{\mathrm{sites}}}=-\sum _{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{c}ln{\langle n^{\left(\gamma \right)}\rangle }_c$$

(146)

is the mixing entropy of the three species using the mean-field occupation numbers and

$${\overline{S}}_c^{\left(B\right)}=\frac{S_c^{\left(B\right)}}{k_B{\mathcal{N}}_{\mathrm{sites}}}=\frac{{\rho }_c^2}{2}{\beta }^2\frac{\partial {\mathbb{S}}_{\mathrm{lattice}}}{\partial \beta }$$

(147)

arises directly from the self-energy $\Sigma \propto {\mathbb{S}}_{\mathrm{lattice}}$. The occupation numbers ${\langle n^{\left(\gamma \right)}\rangle }_c$ used in this expression are those of the mean-field level, given by Equation (136) and shown in Figure 15. The mean-field entropy ${\overline{S}}^{\left(0\right)}$, given in Equation (145), is plotted in Figure 17, together with $-{\langle n^{\left(\gamma \right)}\rangle }_{c}ln{\langle n^{\left(\gamma \right)}\rangle }_c$ for each species. Also shown are the two constituents ${\overline{S}}_c^{\left(A\right)}$ and ${\overline{S}}_c^{\left(B\right)}$ of the total entropy. This plot shows that, as with $S_{NI}$ in Figure 7, the mean-field entropy is greatest for ${\epsilon }_{\Vert }$ near ${\mu }_{\mathrm{dimer}}$ and decreases outside this region. The difference $\delta {\overline{S}}_c=\frac{S_c-S_{\mathrm{NI}}}{k_{B}T}$ between this mean-field entropy for the interacting model and the entropy for the noninteracting model is plotted in Figure 18. Note that this curve has roughly the same shape as that of $\delta {\sigma }_c$ in Figure 16, with a large increase in disorder at large positive ${\epsilon }_{\Vert }$ and small changes elsewhere.

By employing field-theoretic techniques borrowed from quantum mechanics, we have included the interaction energy $U_{\mathrm{int}}$ in the partition function and found a mean-field approximation to the exact partition function. However, these mean-field level calculations assume the same average charge ${\rho }_c$ for all DNA molecules in the same solution. Thus, mean-field calculations would predict a repulsive interaction between two DNA molecules, which would not lead to DNA condensation. Indeed, this repulsive interaction is all that this zeroth-order calculation can predict. It is the fluctuations in the charge from its mean-field value that enable the net attraction observed experimentally; to understand the attractive forces between DNA molecules in solution, it is necessary to extend this treatment by at least one additional order. The next order in the expansion (104) will yield Gaussian-type integrals, which are analytically solvable for the lowest-order fluctuations in the average charge and other parameters. Furthermore, once lowest-order fluctuations have been included, it will be possible to calculate thermodynamic properties of the system from the partition function to meaningful levels, such as the entropy changes due to small fluctuations in the local charge.

8. Inverse Hartree-Field-Fluctuation Propagator

In order to calculate the inverse Hartree-field-fluctuation propagator ${\tilde{D}}^{-1}$, we need to find the second derivatives of the effective action, evaluated at the saddle point. Since ${\tilde{D}}^{-1}$ is diagonal, from Equation (106), we can write

$${\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}=\left(-\frac{1}{\beta {\lambda }_{k_1,{\sigma }_1}}\right){\left.\frac{{\partial }^2{\mathcal{S}}_{\mathrm{eff}}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k_1,{\sigma }_1}\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}\right|}_{{\tilde{\Delta }}_c}.$$

(148)

Since we already calculated the first derivative in Section 7, we need one more derivative, now with respect to ${\tilde{\Delta }}_{k_2,{\sigma }_2}$. Substituting from Equation (121) yields

$$\begin{array}{ccc}\hfill {\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}& =& \left(-\frac{1}{\beta {\lambda }_{k_1,{\sigma }_1}}\right)\frac{\partial }{\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}\left\{-\beta {\lambda }_{k_1,{\sigma }_1}{\tilde{\Delta }}_{k_1,{\sigma }_1}+\right.\hfill \\ & & {\left.\left.+\beta {\lambda }_{k_1,{\sigma }_1}\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k_1,{\sigma }_1),(\mathsf{\ell },b)}^T\frac{{\sum }_{\gamma }{q}_{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}\right\}\right|}_{{\tilde{\Delta }}_c}.\hfill \end{array}$$

(149)

Performing the derivative, we have

$$\begin{array}{ccc}\hfill {\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}& =& \left\{1-\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k_1,{\sigma }_1),(\mathsf{\ell },b)}^T\frac{{\sum }_{\gamma }{q}_{\gamma }\frac{\partial {\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}+\right.\hfill \\ & & {\left.\left.+\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k_1,{\sigma }_1),(\mathsf{\ell },b)}^T\frac{{\sum }_{\gamma }{q}_{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{{\left({\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]\right)}^2}\sum _{{\gamma }^{″}}\frac{\partial {\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{″}\right)}\left[\tilde{\Delta }\right]}{\partial {\tilde{\Delta }}_{k_2,{\sigma }_2}}\right\}\right|}_{{\tilde{\Delta }}_c}.\hfill \end{array}$$

(150)

Substituting the derivative of the activity from Equations (119) and (120), we have

$$\begin{array}{ccc}\hfill {\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}& =& \left\{1+\beta {\lambda }_{k_1,{\sigma }_1}\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k_1,{\sigma }_1),(\mathsf{\ell },b)}^T{\mathcal{W}}_{(k_2,{\sigma }_2),(\mathsf{\ell },b)}^T\times \right.\hfill \\ & & {\left.\left.\times \left[\frac{{\sum }_{\gamma }{q}_{\gamma }^2{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}-{\left(\frac{{\sum }_{\gamma }{q}_{\gamma }{\tilde{a}}_{\mathsf{\ell },b}^{\left(\gamma \right)}\left[\tilde{\Delta }\right]}{{\sum }_{{\gamma }^{\prime }}{\tilde{a}}_{\mathsf{\ell },b}^{\left({\gamma }^{\prime }\right)}\left[\tilde{\Delta }\right]}\right)}^2\right]\right\}\right|}_{{\tilde{\Delta }}_c}.\hfill \end{array}$$

(151)

The quantity in square brackets is just the constant ${\sigma }_c^2$, the variance in the charge per site, given in Equation (137). Since $\mathcal{W}$ from Equations (74)–(76), is an orthogonal matrix,

$$\sum _{\mathsf{\ell },b}{\mathcal{W}}_{(k_1,{\sigma }_1),(\mathsf{\ell },b)}^T{\mathcal{W}}_{(k_2,{\sigma }_2),(\mathsf{\ell },b)}^T={\delta }_{k_1{k}_2}{\delta }_{{\sigma }_1{\sigma }_2},$$

(152)

and this means that the matrix is diagonal. With these substitutions, the inverse Hartree-field-fluctuation propagator becomes

$${\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}={\delta }_{k_1{k}_2}{\delta }_{{\sigma }_1{\sigma }_2}\left(1+{\sigma }_c^2\beta {\lambda }_{k_1,{\sigma }_1}\right).$$

(153)

This is the inverse propagator for the fluctuations of the auxiliary Hartree field that we introduced to decouple the interparticle interactions, expressed in the diagonal basis.

The inverse propagator ${\left({\tilde{D}}^{-1}\right)}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}$ in the diagonal basis is plotted in Figure 19 for the two eigenvalues ${\lambda }_{k\pm }$, shown in Figure 14. Notice that the shapes are those of the eigenvalues, but shifted upward. They both have minima at $k=\pi$, with the one corresponding to ${\lambda }_{k+}$ having its maximum at $k=0$ and the one corresponding to ${\lambda }_{k-}$ having maxima near $k=0$. The corresponding plots of the propagator ${\tilde{D}}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}$ are shown in Figure 20, also in the diagonal basis. The propagator increases with $\beta {\epsilon }_{\Vert }$ near $k=\pm \pi$, and this can be attributed, at least in part, to the increase in ${\sigma }_c$, shown in Figure 16, because of the dependence shown in Equation (151). Although the auxiliary field propagator normally describes the behavior of fluctuations of a physical field like the charge, its definition given by Equation (110) shows that the field here has been rescaled by the square root of the eigenvalue. This was discussed earlier in the context of entropy in the paragraph after Equation (92), and so this form is convenient for the entropy calculation, our focus here, but its interpretation is not quite so direct.

9. Inclusion of Lowest-Order Fluctuations

Going beyond the mean field to one-loop order requires use of the full partition function given by Equations (111) and (113). This has consequences for the site occupancies, the average charge per site, the charge variance, and the entropy. The expression for site occupancy has the same form in terms of the partition function as for the noninteracting and mean-field cases in Equations (9) and (86), and can be written explicitly for the one-loop order of approximation as

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}}=\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta Z_{\mathrm{G}}}\frac{\partial }{\partial {\mu }_{\gamma }}\left[\frac{e^{-{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right]}}{\sqrt{det\left({\tilde{D}}^{-1}\right)}}\right],$$

(154)

which is obtained from Equation (116) with m set equal to one. In that expression, m was the parameter that was introduced in Section 6 as an aid to identifying the various orders in the expansion. This can be written as the sum of two terms,

$${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}}={\langle n^{\left(\gamma \right)}\rangle }_c+\delta \langle n^{\left(\gamma \right)}\rangle ,$$

(155)

where ${\langle n^{\left(\gamma \right)}\rangle }_c$ is the occupancy of species $\gamma$ at the mean-field level, given by Equation (136),

$${\langle n^{\left(\gamma \right)}\rangle }_c=-\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta }\frac{\partial }{\partial {\mu }_{\gamma }}{\mathcal{S}}_{\mathrm{eff}}\left[{\tilde{\Delta }}_c\right],$$

(156)

and $\delta \langle n^{\left(\gamma \right)}\rangle$ is the correction to that mean site occupancy due to fluctuations, which is given by

$$\delta \langle n^{\left(\gamma \right)}\rangle =\frac{1}{{\mathcal{N}}_{\mathrm{sites}}\beta }\frac{\partial }{\partial {\mu }_{\gamma }}ln[det\left({\tilde{D}}^{-1}\right)],$$

(157)

with the inverse propagator ${\tilde{D}}_{(k_1,{\sigma }_1),(k_2,{\sigma }_2)}^{-1}$ from Equation (153) written in matrix form as

$${\tilde{D}}^{-1}=I+{\sigma }_c^2\beta \Lambda ,$$

(158)

where I is the identity matrix and $\Lambda$ is the (diagonal) matrix of eigenvalues of the potential. Since ${\tilde{D}}^{-1}$ is diagonal, its determinant is given by the product of its diagonal elements as

$$det{\tilde{D}}^{-1}=\prod _{k,\sigma }(1+\beta {\lambda }_{k,\sigma }{\sigma }_c^2).$$

(159)

The logarithm of this determinant then gives the sum over logarithms of diagonal elements as

$$ln(det{\tilde{D}}^{-1})=\sum _{k,\sigma }ln(1+\beta {\lambda }_{k,\sigma }{\sigma }_c^2).$$

(160)

The only ${\mu }_{\gamma }$ dependence is contained in ${\sigma }_c^2$, and its derivative is given by

$$\frac{\partial {\sigma }_c^2}{\partial {\mu }_{\gamma }}=\beta {\langle n^{\left(\gamma \right)}\rangle }_c[q_{\gamma }^2-{\langle {\rho }^2\rangle }_c-2{\rho }_c(q_{\gamma }-{\rho }_c)],$$

(161)

so that the derivative of $ln(det{D}^{-1})$ is given by the coefficient

$$\mathbb{C}=\frac{1}{2{\mathcal{N}}_{\mathrm{sites}}}\sum _{k,\sigma }\frac{\beta {\lambda }_{k\sigma }}{[1+\beta {\lambda }_{k,\sigma }{\sigma }_c^2]}.$$

(162)

Then, the mean fluctuation in occupancy $\delta \langle n^{\left(\gamma \right)}\rangle$ becomes

$$\delta \langle n^{\left(\gamma \right)}\rangle =\mathbb{C}{\langle n^{\left(\gamma \right)}\rangle }_c[q_{\gamma }^2-{\langle {\rho }^2\rangle }_c-2{\rho }_c(q_{\gamma }-{\rho }_c)].$$

(163)

These changes in the mean-field occupancies at this one-loop order are plotted in Figure 21. These changes are not large, and the fluctuation corrections amount to no more than about 10% of the total. The most dramatic effects are seen in the region of large positive ${\epsilon }_{\Vert }$, with the corrections to the vacancies and perpendicular dimers going in the opposite direction to the corrections due to the mean field in Figure 15. This can be interpreted as a relaxation of the rigid penalties imposed by the electrostatic self-energies at the mean-field level. An enhancement in parallel dimers is also seen across the whole range of ${\epsilon }_{\Vert }$ that peaks around ${\epsilon }_{\Vert }={\mu }_{\mathrm{dimer}}$.

The mean charge per site at this one-loop-correction level of approximation is given by the same form as Equations (22) and (130) for the noninteracting and mean-field approximations and is given by

$${\rho }_{\mathrm{tot}}={\langle \rho \rangle }_{\mathrm{tot}}=\sum _{\gamma }{q}_{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}},$$

(164)

where now ${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}}$ is taken from Equations (155) and (163). and the charge per site becomes

$${\rho }_{\mathrm{tot}}={\rho }_c+\delta {\rho }_{\mathrm{tot}},$$

(165)

where

$$\delta {\rho }_{\mathrm{tot}}={\rho }_c\mathbb{C}(1-3{\langle {\rho }^2\rangle }_c+2{\rho }_c^2).$$

(166)

The change in the charge per site $\delta {\rho }_{\mathrm{tot}}$ is shown as the solid blue curve in Figure 22. The most notable features are a negative contribution to the total charge above $\beta {\epsilon }_{\Vert }\approx 2$, reflecting the addition of more negatively charged vacancies due to the fluctuations, and a positive contribution to the total charge around ${\epsilon }_{\Vert }\approx {\mu }_{\mathrm{dimer}}$. Interestingly, this positive contribution is not associated with the increase in parallel dimer occupancy, since parallel dimers are electrically neutral. Rather, it is due to the fact that, although both the negatively charged vacancies and positively charged perpendicular dimers are both reduced around ${\epsilon }_{\Vert }\approx {\mu }_{\mathrm{dimer}}$, the vacancies are decreased more.

The variance in the charge is given by

$${\sigma }_{\mathrm{tot}}^2={\langle {\rho }^2\rangle }_{\mathrm{tot}}-{\rho }_{\mathrm{tot}}^2.$$

(167)

The average of the square of the charge per site is given by

$${\langle {\rho }^2\rangle }_{\mathrm{tot}}=\sum _{\gamma }{q}_{\gamma }^2{\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}}$$

(168)

and is written explicitly as

$${\langle {\rho }^2\rangle }_{\mathrm{tot}}={\langle {\rho }^2\rangle }_c+\delta {\langle {\rho }^2\rangle }_{\mathrm{tot}},$$

(169)

where ${\langle {\rho }^2\rangle }_c$ is given by Equation (138), and $\delta {\langle {\rho }^2\rangle }_{\mathrm{tot}}$ works out to be

$$\delta {\langle {\rho }^2\rangle }_{\mathrm{tot}}=\mathbb{C}[{\langle {\rho }^2\rangle }_c(1-{\langle {\rho }^2\rangle }_c-2{\rho }_c^2)-2{\rho }_c^2].$$

(170)

The charge variance can then be written as

$${\sigma }_{\mathrm{tot}}^2={\sigma }_c^2+{\sigma }_f^2,$$

(171)

where the fluctuation correction is

$$\begin{array}{ccc}\hfill {\sigma }_f^2& =& -{\langle {\rho }^2\rangle }_c^2\mathbb{C}(1+9{\rho }_c^2\mathbb{C})+{\langle {\rho }^2\rangle }_c\mathbb{C}[(1+8{\rho }_c^2)+6{\rho }_c^2\mathbb{C}(1+2{\rho }_c^2)]\hfill \\ & -& {\rho }_c^2\mathbb{C}[4(1+{\rho }_c^2)-\mathbb{C}(1+4{\rho }_c^2+4{\rho }_c^4)].\hfill \end{array}$$

(172)

The standard deviation can then be written as

$${\sigma }_{\mathrm{tot}}={\sigma }_c+\delta {\sigma }_{\mathrm{tot}},$$

(173)

The change $\delta {\sigma }_{\mathrm{tot}}$ in the standard deviation of the charges per site due to fluctuations is given in Figure 22, where the fluctuations account for changes typically of the order of 10% of the mean-field value.

Having obtained results for the effects of fluctuations on the mean occupancies and charges per site, we now consider the corresponding corrections to the entropy. For that, we need the grand thermodynamic potential, which was expanded in Section 6 as Equations (115)–(117), with the mean-field and lowest-order fluctuation terms given by

$${\Omega }_G\simeq {\Omega }_c+\frac{1}{2\beta }lndet{D}^{-1}.$$

(174)

The corresponding expression for the interacting entropy is

$$S_{\mathrm{tot}}=S_c+\delta S_{\mathrm{tot}},$$

(175)

where the mean-field entropy $S_c$ was derived and the results presented in Section 7 in Equations (145)–(147). The fluctuation contribution to the entropy $\delta S_{\mathrm{tot}}$ is given by

$$\delta S_{\mathrm{tot}}=k_B{\beta }^2\frac{\partial }{\partial \beta }\left\{\frac{1}{2\beta }ln[det\left(D^{-1}\right)]\right\},$$

(176)

which requires performing another derivative.

To evaluate this contribution to the entropy $\delta S_{\mathrm{tot}}$ in Equation (176), we return to Equation (160) and substitute in the fluctuation contribution to the entropy from Equation (176), which after performing derivatives gives

$$\delta S_{\mathrm{tot}}=S_f^{\left(A\right)}+S_f^{\left(B\right)},$$

(177)

where

$$S_f^{\left(A\right)}=-\frac{k_B}{2}\sum _{k,\sigma }ln(1+\beta {\lambda }_{k,\sigma }{\sigma }_c^2)$$

(178)

and

$$S_f^{\left(B\right)}=\frac{k_B}{2}\sum _{k,\sigma }\frac{1}{(1+\beta {\lambda }_{k,\sigma }{\sigma }_c^2)}(\beta {\lambda }_{k,\sigma }{\sigma }_c^2+{\sigma }_c^2{\beta }^2\frac{\partial {\lambda }_{k,\sigma }}{\partial \beta }+\beta {\lambda }_{k,\sigma }\beta \frac{\partial {\sigma }_c^2}{\partial \beta }).$$

(179)

The derivative of the eigenvalue ${\lambda }_{k,\sigma }$ of the potential is given by

$$\frac{\partial {\lambda }_{k,\sigma }}{\partial \beta }=\frac{\partial }{\partial \beta }{\tilde{V}}_{\uparrow \uparrow }\left(k\right)+\sigma \frac{\partial }{\partial \beta }\left|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right|,$$

(180)

where ${\tilde{V}}_{\uparrow \uparrow }\left(k\right)$ is the Fourier transform of the interaction potential of a chain with itself and ${\tilde{V}}_{\uparrow \downarrow }\left(k\right)$ is the Fourier transform of the interaction potential between the two chains. Both of these potentials depend on the screening vector $q_s$, and the only $\beta$ dependence of the potential is contained in $q_s$, which is proportional to $\sqrt{\beta }$. The derivative we need is then given by Equation (144). As a consequence, the derivative of the Fourier transform of the potential is

$$\frac{\partial {\tilde{V}}_{b_1,b_2}\left(k\right)}{\partial \beta }=\frac{1}{2}\sum _{\mathsf{\ell }=-\mathcal{N}}^{\mathcal{N}}\left[e^{-ik\mathsf{\ell }}\frac{\partial V_{b_1,b_2}\left(\mathsf{\ell }\right)}{\partial \beta }+e^{ik\mathsf{\ell }}\frac{\partial V_{b_2,b_1}\left(\mathsf{\ell }\right)}{\partial \beta }\right],$$

(181)

where the derivative of the potential in the position basis is given by

$$\frac{\partial V_{b_1,b_2}\left(\mathsf{\ell }\right)}{\partial \beta }=V_{b_1,b_2}\left(\mathsf{\ell }\right)\left(-\frac{q_s}{2\beta }{d}_{b_1,b_2}\left(\mathsf{\ell }\right)\right).$$

(182)

Substituting and using the symmetry relations in Equations (47) and (50), the derivative of the Fourier transform of the potential is

$$\frac{\partial {\tilde{V}}_{b_1,b_2}\left(k\right)}{\partial \beta }=-\frac{q_s}{2\beta }\sum _{\mathsf{\ell }=-\mathcal{N}}^{\mathcal{N}}{e}^{-ik\mathsf{\ell }}{V}_{b_1,b_2}\left(\mathsf{\ell }\right)d_{b_1,b_2}\left(\mathsf{\ell }\right).$$

(183)

Using these relations, the derivative of the eigenvalues becomes

$$\begin{array}{ccc}\hfill {\beta }^2\frac{\partial {\lambda }_{k,\sigma }}{\partial \beta }& =& -\frac{q_s}{2}\sum _{\mathsf{\ell }=-\mathcal{N}}^{\mathcal{N}}\left\{\beta V_{\uparrow \uparrow }\left(\mathsf{\ell }\right)d_{\uparrow \uparrow }\left(\mathsf{\ell }\right)cos\left(k\mathsf{\ell }\right)\right.\hfill \\ & & \left.+\sigma V_{\uparrow \downarrow }\left(\mathsf{\ell }\right)d_{\uparrow \downarrow }\left(\mathsf{\ell }\right)\left[\frac{\mathrm{Re}\left[{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right]}{|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)|}cos\left(k\mathsf{\ell }\right)-\frac{\mathrm{Im}\left[{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right]}{|{\tilde{V}}_{\uparrow \downarrow }\left(k\right)|}sin\left(k\mathsf{\ell }\right)\right]\right\},\hfill \end{array}$$

(184)

where $\mathrm{Re}\left[{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right]$ and $\mathrm{Im}\left[{\tilde{V}}_{\uparrow \downarrow }\left(k\right)\right]$ are the real and imaginary parts of ${\tilde{V}}_{\uparrow \downarrow }\left(k\right)$.

Equation (177) also requires the derivative of the variance ${\sigma }_c^2$, where ${\sigma }_c^2$ is given by Equation (137). This contains the average charge per site from Equation (130) and the average of the square of the charge per site from Equation (138). This derivative is

$$\begin{array}{ccc}\hfill \beta \frac{\partial {\sigma }_c^2}{\partial \beta }& =& \frac{{\rho }_c}{1+{\sigma }_c^2\beta {\mathbb{S}}_{\mathrm{lattice}}}\left\{\left({\left[\rho \mathcal{E}\right]}_{\mathrm{av}}-{\rho }_c{\mathcal{E}}_c\right)\left[2+(1-{\langle {\rho }^2\rangle }_c)\beta {\mathbb{S}}_{\mathrm{lattice}}\right]+\right.\hfill \\ & & +\left.{\rho }_c\left(3{\sigma }_c^2+{\rho }_c^2-1\right){\beta }^2\frac{\partial {\mathbb{S}}_{\mathrm{lattice}}}{\partial \beta }\right\}-\left({\left[{\rho }^2\mathcal{E}\right]}_{\mathrm{av}}-{\langle {\rho }^2\rangle }_c{\mathcal{E}}_c\right),\hfill \end{array}$$

(185)

where

$$\mathcal{E}\equiv \beta ({\epsilon }_{\gamma }-{\mu }_{\gamma }+{\rho }_c{q}_{\gamma }{\mathbb{S}}_{\mathrm{lattice}}),$$

(186)

$${\mathcal{E}}_c=\sum _{\gamma }\mathcal{E}{\langle n^{\left(\gamma \right)}\rangle }_c,$$

(187)

$${\left[\rho \mathcal{E}\right]}_{\mathrm{av}}=\sum _{\gamma }{q}_{\gamma }\mathcal{E}{\langle n^{\left(\gamma \right)}\rangle }_c,$$

(188)

and

$${\left[{\rho }^2\mathcal{E}\right]}_{\mathrm{av}}=\sum _{\gamma }{q}_{\gamma }^2\mathcal{E}{\langle n^{\left(\gamma \right)}\rangle }_c.$$

(189)

Here, ${[\cdots ]}_{\mathrm{av}}$ is the weighted average of the specified quantity over the distribution given by the mean-field occupancies ${\langle n^{\left(\gamma \right)}\rangle }_c$ of the three species.

It is interesting to note that the combinations ${\left[{\rho }^2\mathcal{E}\right]}_{\mathrm{av}}-{\langle {\rho }^2\rangle }_c{\mathcal{E}}_c$ and ${\left[\rho \mathcal{E}\right]}_{\mathrm{av}}-{\rho }_c{\mathcal{E}}_c$ can be written in forms similar to that of the zeroth-order entropy. These forms are

$${\left[{\rho }^2\mathcal{E}\right]}_{\mathrm{av}}-{\langle {\rho }^2\rangle }_c{\mathcal{E}}_c=-\sum _{\gamma }\left[{\langle n^{\left(\gamma \right)}\rangle }_{c}ln{\langle n^{\left(\gamma \right)}\rangle }_c\left(1-\sum _{{\gamma }^{\prime }}{\langle n^{\left({\gamma }^{\prime }\right)}\rangle }_{c}ln{\langle n^{\left({\gamma }^{\prime }\right)}\rangle }_c\right)\right]$$

(190)

and

$${\left[\rho \mathcal{E}\right]}_{\mathrm{av}}-{\rho }_c{\mathcal{E}}_c=-\sum _{\gamma }{q}_{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{c}ln{\langle n^{\left(\gamma \right)}\rangle }_c+{\rho }_c\sum _{\gamma }{\langle n^{\left(\gamma \right)}\rangle }_{c}ln{\langle n^{\left(\gamma \right)}\rangle }_c.$$

(191)

The fluctuation entropy contribution $S_f^{\left(A\right)}$ is easy to calculate, given the eigenvalues ${\lambda }_{k,\sigma }$. Equations (180)–(191) simply collect the results that allow us to calculate $S_f^{\left(B\right)}$. The results for the total fluctuation entropy $\delta S_{\mathrm{tot}}$ and the two contributions to it are shown in Figure 23 for ${\epsilon }_{\perp }=\frac{1}{2}{\epsilon }_{\Vert }$. Like the mean-field entropy, the fluctuation entropy becomes larger in magnitude as the binding becomes weak. However, the fluctuation contribution is actually negative in this region, suggesting that the fluctuations are leading to increasing order. These corrections to the entropy are typically an order of magnitude larger than the difference between the mean field and noninteracting values of the entropy.

10. Discussion

Entropy forms an important contribution to the free energies of many biological systems. Here, we summarize the results for the entropies and particle numbers per site for a lattice gas model representing the adsorption of molecules on double-stranded DNA. The entropies and particles per site are presented at three successive levels of approximation, first, for the noninteracting system, then at the mean-field level, and finally with the inclusion of fluctuations to one-loop order. Those fluctuation corrections result from correlations induced by the electrostatic forces among the dimer molecules themselves and with the electrically charged DNA substrate.

Perhaps the most transparent descriptors of the DNA system’s behavior are the average occupation numbers $\langle n^{(\Vert )}\rangle$, $\langle n^{(\perp )}\rangle$, and $\langle n^{\left(v\right)}\rangle$ representing the thermal average of the numbers of parallel-adsorbed dimers, perpendicular-adsorbed dimers, and vacancies, respectively, on each site. These particles per site are proportional to the probability that a given site would be found to be occupied by a particular species. These site-occupation numbers are shown in Figure 24 for the three levels of approximation. The short-dashed curves are for the noninteracting level of approximation ${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{NI}}$, the long-dashed curves for the mean-field level ${\langle n^{\left(\gamma \right)}\rangle }_c$, and the solid curves for the inclusion of fluctuations ${\langle n^{\left(\gamma \right)}\rangle }_{\mathrm{tot}}$. The corresponding charges per site and standard deviations of the charge density are shown in Figure 25.

The site occupancies exhibit similar qualitative behavior in all three levels of approximation. Due to the energetics of binding ${\epsilon }_{\Vert }=2{\epsilon }_{\perp }$, parallel dimers occupy the greatest fraction of sites for all negative binding energies $({\epsilon }_{\Vert },{\epsilon }_{\perp }<0)$, followed by the perpendicular dimers, and then the vacancies, at least for the parameters used in these calculations. At weak binding (${\epsilon }_{\Vert },{\epsilon }_{\perp }\to 0$) in the noninteracting system, the difference between these energies becomes negligible, and both $n^{(\Vert )}$ and $n^{(\perp )}$ approach the same value of $40\%$. When the binding energy ${\epsilon }_{\Vert }$ or ${\epsilon }_{\perp }$ increases and reaches the chemical potential ${\mu }_{\mathrm{dimer}}$ of the dimers in solution, they are repelled from the surface of the DNA and go into solution, leaving a “naked” lattice of negatively charged vacancies, with the parallel dimers being more repelled than perpendicular ones.

In contrast to the case of noninteracting sites, the inclusion of interactions between sites gives a pronounced gap of about $25\%$ by which parallel adsorption exceeds perpendicular adsorption in the weak-binding limit. This enhancement in parallel adsorption in the mean-field approximation results from the inclusion of the charge on the perpendicular dimers in an average or mean-field sense, which penalizes the addition of positively charged perpendicular dimers when the total charge on the DNA lattice is positive. Fluctuations reduce this effect slightly. Recall that the perpendicular dimers are positively charged because one end is not attached to a negative binding site, and this means that it is energetically unfavorable to have a second one nearby. Parallel binding, on the other hand, neutralizes the negative charge on the sites, so that a DNA double strand completely covered in parallel-adsorbed dimers would be electrically neutral. The decrease in perpendicular-dimer occupancy caused by this electrostatic repulsion of like charges leads to an increase in sites occupied by vacancies. Vacancies are negatively charged and are energetically favored when the perpendicular dimers provide a positively charged region. There may, in fact, be correlations resulting from the lower electrostatic energy of a configuration in which sites alternate between perpendicular dimers and vacancies, i.e., between positively and negatively charged sites, in a Wigner-lattice-like state similar to those described in the literature [8,19,26].

The average charge ${\rho }_c$ per site depends only on the perpendicular dimers and vacancies since the parallel dimers neutralize sites, so that the charge per site is written simply as

$$\rho =\sum _{\gamma }{q}_{\gamma }\langle n^{\left(\gamma \right)}\rangle =\langle n^{(\perp )}\rangle -\langle n^{\left(v\right)}\rangle ,$$

(192)

and this connection can be seen in the plots of charge per site $\rho$ in Figure 25, when compared with Figure 24. For all approximations, the average charge is positive at negative and slightly positive binding energies, indicating charge inversion. This is a result of the occupation numbers shown in Figure 24, because parallel binding, which neutralizes the charge, increasingly dominates in the strong-binding limit (${\epsilon }_{\Vert }=2{\epsilon }_{\perp }\ll 0$). Similarly, the decrease in the magnitude of the charge inversion from the noninteracting to the mean-field approximations can be explained by the enhancement in parallel adsorption and reduction of perpendicular adsorption due to like-charge repulsion. The fluctuation corrections to the site occupancies are not large when compared with the mean-field values. They only seem significant at positive binding energies where the site occupancies are shifted in the direction of the noninteracting system. Of course, the site occupancies are single-particle properties and do not measure spatial correlations.

The behavior of the entropy is shown in Figure 26, in which the entropy per site is plotted versus $\beta {\epsilon }_{\Vert }$ for all three levels of approximation. The three curves are practically indistinguishable for strong binding (negative values of $\beta {\epsilon }_{\Vert }$) until the binding energy becomes quite weak (approaching zero), indicating little effect of correlations induced by electrostatic repulsion on the entropy in that region. This is because the dimers are closely packed, leaving little freedom for disorder. The pronounced differences are seen in the region of ${\epsilon }_{\Vert }\approx {\mu }_{\mathrm{dimer}}$ and more positive binding energies (${\epsilon }_{\Vert }=2{\epsilon }_{\perp }\gg 0$), because the system becomes more disordered. There, the dimers tend to leave the surface and enter the solution, and the fluctuation corrections to the entropy are greater in that region and act to reduce the entropy. This shows that the electrostatic correlations increase the order and reduce the entropy. Surprisingly, this decrease in entropy at large positive ${\epsilon }_{\Vert }$ due to fluctuations is so significant that it overcompensates for the entropy generated at the mean-field level, leading to a total entropy $S_{\mathrm{tot}}$ that is even lower than in the noninteracting case $S_{\mathrm{NI}}$.

The form of the noninteracting entropy comes from the lattice gas model, in which parallel adsorption is treated as “local”, only occupying one site. When the geometrical blocking effect is included, in which parallel-adsorbed dimers occupy two sites, we would expect the entropy to differ significantly from the $nlnn$ result even at the noninteracting level. It is unclear how large a contribution to the full entropy the blocking effect provides, but it is certainly lower for dimers than for longer polymers, for which the nonlocality is greater. Because of this, it will be important in future work to develop a way of incorporating the blocking effect into these types of field-theoretic models.

Dimers are certainly the simplest example of polymers, although biological polymers are typically considerably longer. There has been extensive work on the dimer model, beginning with the work of Fisher [52,53]. Some of this work has been motivated by the analogy between dimer models and quantum spin systems [36]. In this context, Fisher has derived expressions for the contributions of hard-core crowding effects to the partition function, including the geometrical blocking effect. Fisher’s results include an expression of the noninteracting partition function for dimers on a 1D linear lattice. In a future publication, we plan to adapt Fisher’s approach to study the statistical consequences of these blocking effects [54].

In addition, it has been suggested that methods can be adapted from analogous problems in particle physics [37] and condensed matter physics [55]. This involves using numerical codes like those of Adams and Chandrasekharan [37] that are used for simulations in lattice quantum chromodynamics.

Modeling the behavior of DNA and other biologically active polymeric molecules under physiological conditions is an area in which quantitative calculation is particularly difficult. This is because it is necessary to accommodate simultaneously the influences of geometry, electrostatic forces, and charge correlations. All of these play roles of varying significance in determining the properties of these molecules. Despite these challenges, considerable progress has been made in extending our physical understanding of these systems, and some surprising new physics, such as charge inversion and condensation, have been discovered in the process.

The work presented here focuses on electrostatic effects and especially on the role they play in the entropy, with the overall goal of determining the entropy as a contribution to the free energy and as a driving force in biochemical reactions. For the DNA-based model that we consider, we find that, particularly in the low-coverage regime, a uniform mean-field theory gives changes ∼20% from the usual noninteracting entropy term. We also find, though, that the fluctuation contributions are both considerably larger and in the opposite direction. This establishes that the fluctuations are significant, and the opposite sign of their contribution suggests that they lead to additional order in the system.

Several aspects of the formalism and calculations presented here have been performed in such a way as to be accessible to students and others wishing to extend these techniques to more realistic models of polyelectrolytes condensing on the surfaces of DNA. We have shown in Section 3 why the DNA helices can be treated as a pair of one-dimensional lattices (Figure 9), albeit with unusual interactions due to the three-dimensional nature of the screened Coulomb potential, as shown in Figure 12. We show that by writing the interaction term in the diagonal basis, as in Equation (71) and Figure 14, the problem can be greatly simplified. In addition, we have chosen a basis in which the eigenvalues and eigenvectors are real, in Equations (74) and (75), which makes the physical interpretation much more transparent. Then, we show in detail in Equation (87) and Appendix A how to perform a Hubbard–Stratonovich transformation of the partition function to decouple the sum over configurations from the sum over lattice sites, which allows one to perform the sum over all configurations of the system exactly. This introduces an integral over an auxiliary field $\tilde{\Delta }$. Often, such an auxiliary field emerges as a complicated complex matrix, but in our diagonal representation it is real and represents the charge per site in the diagonal basis, which was shown explicitly for the mean field by comparing Equations (123) and (130). Also, in this diagonal basis, it becomes clear how to choose a proper integration path to have a converging integral. One must then expand the effective action in a power series in the fluctuations of the auxiliary field $\delta \tilde{\Delta }$, as in Equation (104), which physically is in powers of the fluctuations of the charge per site about the mean field. Setting the first-derivative term in this expansion to zero yields the mean field, which itself includes terms of the interaction to all orders, as can be seen in Equation (122) or (123). Whether the expansion is useful depends on the size of the charge fluctuations. The series is an asymptotic one, which means that the infinite series does not necessarily converge. However, as is true for all asymptotic expansions, it is useful as long as the last term kept in the expansion is smaller than the term before it, which is true in our case. This indicates that the mean field, which is the saddle point, is a reasonable first approximation of the system. The fluctuations about this mean field then give insight into how the actual solution differs from the mean field.

In conclusion, we have shown with this simple model the importance of including fluctuation contributions beyond the mean field. The method we have used has allowed us to include vacancies, and we have found that the inclusion of the fluctuation term decreases the calculated entropy by ∼50% in the weak-binding regime. This is because, due to the presence of vacancies, the bound dimer concentration is low because of the competition between dimers being repelled from the DNA molecule and the chemical potential driving them from the solution onto the DNA surface. It is surprising that this decrease in entropy due to correlations is so significant that it overcompensates for the entropy increase at the mean-field level, so that the total entropy is even lower than in the absence of interactions between lattice sites. This effect of fluctuations could be important for other biological systems as well, and this approach can be useful where correlations are significant.