Dipolar Dimer Liquid

TL;DR

This paper introduces a model of dipolar dimer liquid on a 2D lattice, revealing a partially ordered glacia phase at high density and low temperature, with insights into entropy and phase transitions relevant to hydrophobic effects.

Contribution

It proposes a novel dipolar dimer liquid model, maps phase transitions to an annealed Ising model, and explores entropy effects related to polymers and biological phenomena.

Findings

Identification of a glacia phase at high density and low temperature.

Critical temperature estimation aligning with Ising model bounds.

Power-law suppression of configurational entropy with polymer length.

Abstract

A model of dipolar dimer liquid (DDL) on a two-dimensional lattice has been proposed. We found that at high density and low temperature, it has a partially ordered phase which we called glacia phase. The glacia phase transition can be understood by mapping the DDL to an annealed Ising model on random graphs. In the high density limit the critical temperature obtained by the Monte Carlo simulation is $k_BT_c^G = (3.5\pm0.1)J$, which agrees with the estimations of the upper and lower bounds of $k_BT_c^G$ with exactly solved Ising models. In the high density and low temperature limit, we further studied configurational entropy of the DDL in the presence of the neutral polymers. The suppression of the configurational entropy scales as a power law of the polymer length $\lambda_p^\alpha$ with $\alpha \geq 1$, which implies that the configurational entropy of water plays essential roles in understanding the hydrophobic effect and the protein folding problem.

Figures

Equations

$$\begin{array}{rl}{H}_Q=& \sum _{e\in E;v_1,v_2\in \partial e}J(e|C)Q(v_1)Q(v_2),\\ J(e|C)=& \{\begin{array}{ll}+\infty & \text{ for e∈C,}\\ J>0& \text{ for e∈C,}\end{array}\end{array}$$

$$\begin{array}{r}{H}_D=\frac{{\mathbf{d}}_1\cdot {\mathbf{d}}_2-3({\mathbf{d}}_1\cdot {\mathbf{e}}_r)({\mathbf{d}}_1\cdot {\mathbf{e}}_r)}{r^3},\end{array}$$

$$\begin{array}{r}\Phi ={\rho }_r-{\rho }_b\end{array}$$

$$\begin{array}{rl}Z=& \sum _C\sum _{\{\sigma \}|C}{\mathrm{e}}^{-\beta H_{\sigma }}=\sum _C{Z}_C,\\ H_{\sigma }=& \sum _{b\in \tilde{E};{\tilde{v}}_1,{\tilde{v}}_2\in \partial b}(-J){\sigma }_{{\tilde{v}}_1}{\sigma }_{{\tilde{v}}_2}.\end{array}$$

$$\begin{array}{r}M=\frac{1}{||L||}\sum _{\tilde{v}\in \tilde{V}}{\sigma }_{\tilde{v}},\end{array}$$

$$\begin{array}{r}{S}_G(L)=k_B\ln W(L),\end{array}$$

$$\begin{array}{rl}& S_{DDL}(L|{\Omega }_p)=S_G(L|{\Omega }_p)+S_{\text{C}}\\ =& k_B\ln W(L|{\Omega }_p)+N_D{k}_B\ln 2,\end{array}$$

$$\begin{array}{rl}{S}_{\text{Config}}({\Omega }_p)=& \underset{L\to \infty }{\lim }(S_G(L)-S_G(L|{\Omega }_p))\\ & +(1-N_D)k_B\ln 2.\end{array}$$

Full text

Introduction.– Dimers on lattices have been constant topics of statistical physics Kasteleyn (1961); Temperley and Fisher (1961); Fisher (1961); Kasteleyn (1963); Lieb (1967); Alet et al. (2005); Dickman (2012); Otsuka (2011) and mathematics Thurston (1990); Elkies et al. (1992). They have been used for understanding various phenomena, e.g., surface adsorptions Retter (1987), liquid-vapor transitions Lebowitz and Gallavotti (1971); Dickman (2012), liquid crystals Helimann and Lieb (1979); Abraham and Heilmann (1980), superconductivity Rokhsar and Kivelson (1988), roughening transitions Blote and Hilborst (1982); Levitov (1990), etc. Nevertheless they are also simple enough so that they are exactly solvable by various mathods Kasteleyn (1961); Temperley and Fisher (1961); Fisher (1961); Kasteleyn (1963); Lieb (1967); Elkies et al. (1992).

It is natural to consider the dimers carrying dipole moments for better describing the molecules of the polar solvent, like water ($H_{2}O$) and hydrofluride ($HF$). In addition to the hard-core repulsion, the leading term of the interactions is the electrostatic Coulomb forces due to the polarized charges. (In chemistry, this effect is often referred to as hydrogen bond.)

In this article, we considered a model of dipolar dimers with electrostatic interactions, which we called dipolar dimer liquid (DDL). We found that it has two liquid phases. In the high density and low temperature phase, the charges of the dimers are ordered while no translational symmetry is broken for the dimer configurations, which we called glacia phase. The phase transition may be understood by mapping the the DDL to an annealed Ising model on random graphs. In the high density and low temperature limit, the system has massively degenerated ground states while the charge degrees of freedom are completely frozen. Therefore it allows to calculate the ground-state entropy exactly Kasteleyn (1961); Lieb (1967); de Tiliere ; Jacobsen (2016).

A particularly interesting property of water is the hydrophobic effect. In the text books, a famous rule of thumb in chemistry, “like dissolves like”, is usually attributed to the disruption of the hydrogen bonds when a nonpolar molecule is present in the water (enthalpy effect). A well known example is that oil does not like to dissolve in water. The hydrophobic effect has also been recognized as a key ingredient for understanding the protein folding problem which remains an important open problem in biochemistry over decades. However, it may not be properly understood without including the entropy Silverstein (1998). The entropy contributions come from both the conformations of the proteins and the configurations of the solvent Kauzmann (1959); Brady and Sharp (1997); Schafer et al. (2001); Miyamoto and Kollman (1993). Much effort has been made on conformation entropy of the proteins Doig and Sternberg (1995); Zhang and Liu (2006); Baxa et al. (2014); Chang et al. (2007) but fewer on the configurations of the solvent Breiten et al. (2013).

We calculated the configurational entropy of the DDL in the high density and low temperature limit in the presence of neutral polymers of different conformations. Surprisingly, we found the suppression of the configurational entropy of the DDL depends strongly on the polymer conformations. The dimers surrounding the polymers are more “ordered” although there are no liquid-solid phase transitions. It provides a quantitative understanding of the picture of forming water “cages” around the polymers in the theory of hydrophobic effect.

The DDL is undoubtedly an over-simplified model for the realistic polar solvent. Nevertheless, in such a simple model, it is much easier to identify and distinguish various effects observed in more complicated systems, e.g., the proteins in the water. As one may map the DDL to the Ising models, and it is exactly solvable in some limits, which provides a bench mark for more realistic models and numerical simulations.

Definition of models.– Let $L=\{V,E\}$ be a two-dimensional lattice, where $V$ and $E$ are vertex set and edge set respectively. Each dimer sits on an edge $e\in E$ of the lattice and it covers two neighboring vertices $v_{1},v_{2}\in V$ linked by the edge $e$. Each dimer is charge neutral in total but its two ends carry an equal amount of charges but of opposite signs, therefore the dimer has a dipole moment $d_{e}=Ql_{e}$, where $l_{e}$ is the length of the dimer. The dimers are not allowed to overlap on top of each other and they may also have electric interactions due to the polarized charges and dipole moments. This serves a very general definition of the DDL. In order to capture some characters of the realistic systems and to be mathematically handleable, we shall specify some features in more details.

Since the microscopic molecules are identical, it is natural to assume that the lattice holding the dimers has its edges of the same length $l_{0}$. We shall also assume that the lattice is bipartite, i.e., there exist two subsets $A,B$ of $V$ such that $A\cup B=V,A\cap B=\emptyset$, and that for each edge $e$, it has one vertex in $A$ and the other in $B$. Particularly, in this article we shall focus on the case of the two dimensional square lattice. Denote $C$ the set of the edges covered by a dimer for a given configuration of dimer-covering. The molecules of the polar solvents have permanent dipole moments. One assumes that each dimer carries a dipole moment of the magnitude $d_{0}$, therefore $Q_{0}=d_{0}/l_{0}$ is the amount of the charges that each vertex of a dimer has.

The leading term of the interactions is the Coulomb interaction. For each $v$ vertex covered by a dimer it has a charge $Q(v)=\pm Q_{0}$; otherwise it has no charge $Q(v)=0$. We neglect the long-range tail of the Coulomb potential and only consider the nearest neighbor interaction. So the Hamiltonian for the Coulomb interaction is

[TABLE]

where $J$ is a positive coupling constant.

This term reflects the partially-covalent nature of the hydrogen bond. The number of the hydrogen bond an electronegative atom can have is determined by its alone pairs, e.g., $O$ in $H_{2}O$ may form two hydrogen bonds and $F$ in $HF$ may form three. For the dimers on the square lattice, an atom on a vertex is covalently bonded to one of its neighbors, and it can form hydrogen bond with three more unbonded neighbors. To take the angled geometry of the water molecules into account, one may improve the model by putting angled trimers on a Lieb lattice  111The oxygen atom sits on a the degree-four vertex and its two orthogonal pods lie on the neighboring edge centers., which has been reported experimentally in Ref. 31

The subleading electrostatic interaction is the dipole-dipole interaction

[TABLE]

where $\mathbf{r}$ is the vector from one dipole to the other and $\mathbf{e}_{r}$ is the unit vector along the direction of $\mathbf{r}$. The dipole interaction tends to align two dimers in some preferred directions relative to each other. Some lattice dimer models with neutral dimers and direction-dependent interactions have been discussed by Heilmann and Lieb Helimann and Lieb (1979) for describing the liquid crystals. In Ref. 13, it has been proved that at high enough fugacity and at low enough temperature there is a nematic phase. In the simple solvents like water or hydrogen fluoride, there is no observed evidence of the existence of such a nematic phase, so we neglect the dipole-dipole interaction.

One may also take the ionization effect into account by introducing charged monomers or trimers ($OH^{-}$, $H_{3}O^{+}$ for water and $HF_{2}^{-}$, $H_{2}F^{+}$ for hydrogen fluoride). As $pK_{w}=14$ 222$K_{w}$ is the dissociation constant of the water self-ionization. around the room temperature, the fraction of the ionized water molecules is negligible. When the salts, like the strong electrolyte $NaCl$, are added, the density of the ions and cations in the solution is no longer negligible. It is then necessary to consider the charged monomer-dimer systems instead of the simple DDL.

Therefore we shall consider a DDL on a square lattice $L$ with the short range Coulomb interaction given by Eq. 1. Although it is such a simplified model compared to the real systems, we believe that it in some sense reasonably approximates the polar solvents. On the other hand, this model is simple enough so that we may solve it exactly in some limits. We will show that the liquid has two phases. In the high density and low temperature limit, it reduces to a dimer-covering problem, which allows us to calculate the configurational entropy even in the presence of a neutral polymer.

Phases of dipolar dimer liquid.– The dimer liquid defined in the previous section has three parameters, the coupling constant $J$, the temperature $T$ and the density of the dimers $\rho$.

When the dimer density $\rho$ is small, the nearest-neighbor Coulomb interaction between the dimers is negligible and the no-overlap constraint can be considered as a short range hard core repulsion. So the system is in the vapor phase as described by the Widom-Rowlinson model Widom and Rowlinson (1970). The existence of a liquid-vapor phase transition has been proved in the perspective of grand canonical ensemble by Lebowitz and Gallavotti Lebowitz and Gallavotti (1971) as well as Ruelle Ruelle (1971) by using a Peierls contour argument. In the high density limit the system is in the liquid phase.

In the high density limit, the existence of a charge order phase transition can be made in analogy to the model of two-component Widom-Rowlinson mixture. 333This picture was pointed to the author by Elliott Lieb. A vertex $v\in V$ of the square lattice $L$ is labelled by its row-column coordinate $(m,n)$, we defined its signature $\sigma(v)=(-1)^{(m+n)}$, which is consistent with $A$-$B$ partition of the lattice. If a dimer with its positive charge on a vertex in set $A$, we call it a dimer of $A$-type and color it in red; if a dimer with its positive charge on a vertex in set $B$, we call it a dimer of $B$-type and color it in blue. As the lattice $L$ is bipartite, any dime must be either $A$-type or $B$-type.

We call two dimers are neighbors, if there exists an edge $e$ whose two ends are one vertex of each dimer respectively. There are four kinds of neighbors: 1) a pair of collinear dimers, 2) a pair of perpendicular dimers, 3) a pair of parallel dimers not on the same square plaquet, 4) a pair of parallel dimers sitting on the opposite edges of one sqaure plaquet. Then the dimers of the same type in a configuration of neighbor kind 1)–3) attract each other with strength $J$, and in a configuration of neighbor kind 4) with $2J$; the dimers of different types in a configuration of neighbor kind 1)–3) or neighbor kind 4) repel each other with $J$ or $2J$ respectively. The no-overlap constraint leads to a hard core repulsion.

The order parameter is defined as the difference of the densities of red and blue dimers

[TABLE]

It follows the Peierls contour argument Lebowitz and Gallavotti (1971); Ruelle (1971) that there is a red-blue symmetry broken phase when the fugacities $z_{A}=z_{B}$ is high enough and the dimensionless coupling $\beta J$ is strong enough. We call this phase glacia phase. 444Glacia means ice in Latin. In the glacia phase, the charge of the dimers are ordered, while no translational symmetry is broken for the dimer configurations. Mapping the DDL to the Ising model shows the system looks similar to a glass model.

An more quantitative way to understand this phase transition is to map the DDL to an annealed Ising model on random graphs. Consider the Lieb lattice $\tilde{L}$ associated to the square lattice $L$, i.e., adding a site $\tilde{v}_{e}$ in the middle of each edge $e$. For a dimer of type $A$ on a edge $e$, we associate a spin pointing upwards to the Lieb site $\tilde{v}_{e}$, and for type $B$, we associate a spin pointing downwards to $\tilde{v}_{e}$. The no-overlap constraint of the dimers is interpreted as no two spins are put on the nearest-neighboring Lieb sites. For two neighboring dimers of kind 1)–3), one links the corresponding Lieb sites with a single bond; for two neighboring dimers of kind 4), one links the corresponding Lieb sites with a double bond. For a given configuration of the dipolar dimers, the Lieb sites assigned with a spin and the bonds linking the Lieb sites form a random graph $G=\{\tilde{V}^{L},\tilde{E}\}$. For each bond, we assign an spin-spin interaction $(-J)\sigma_{\tilde{v}_{1}}\sigma_{\tilde{v}_{2}}$. (For a double bond, one counts the bond multiplicity, i.e., in total $2J$ for the spins connected by the double bond). The partition function is

[TABLE]

In general, one may add a weight factor $W_{C}$ before the partition function $Z_{C}$ of the Ising model on the random graph $G_{C}$ associated to the dimer configuration $C$. The partition function $Z_{C}$ is averaged over configurations of dimers, which we call annealed Ising model on random graphs. 555It is a convention in the spin glass that “annealing” refers to averaging the partition function over the disorders. The order parameter in this picture, i.e., the magnetization for the Ising model

[TABLE]

is the same as the red-blue order parameter $\Phi$ defined for the Widom-Rowlinson mixture. Instead of the fugacity in the picture of Widom-Rowlinson mixture, it is easier to used the parameter of dimer density $\rho$ in the picture of Ising model.

When the density of the dimers is low, the random graphs are sparse, the spins are not coupled, so there is no spontaneous magnetization. When the density of the dimers is higher, the spins on the random graphs are more strongly coupled, although there might not be spontaneous magnetization at a finite temperature. For example, if on average a dimer has less than two neighbors, comparing to the one-dimensional Ising model of which the critical temperature $T_{c}=0$, one would not expect a spontaneous magnetization at a finite temperature. When the density is even higher, there is a finite critical temperature $T_{c}$, below which the DDL is in the glacia phase.

The high density limit means all the vertices of the lattice $L$ are occupied by some dimer, i.e., the dimers completely cover the lattice $L$ without any monomer or vacancy. Then each vertex of a random graph $G_{C}$ for a given dimer configuration has degree $6$ counting bond multiplicity, although the number of its neighbors may various. 666A vertex can have: 1) six neighbors all single-bonded, or 2) five neighbors with one double-bonded, or 3) four neighbors with two double-bonded. It is heuristic to compare this annealed Ising model on random graphs to the two dimensional Ising models on the triangular lattice and square lattice, where each site has six or four neighbors but no double bonds. Each spin on a random graph has at least four neighbors with couplings $J$ or $2J$. So the critical temperature of the glacia transition should be higher than that of the square lattice, and comparable to that of the triangular lattice $T_{c}^{\Box}<T_{c}^{G}\sim T_{c}^{\Delta}$. More precisely, for an excitation of two neighboring spins, on triangular lattice it costs $10J$; on the random graphs, it costs $10J$ or $8J$, which depends on the bond multiplicity. So $T_{c}^{G}\lesssim T_{c}^{\Delta}$. For some graphs with as many double bonds as possible, it should be comparable with anisotropic Ising model on a square lattice with $J_{h}=J$ and $J_{v}=2J$. The critical temperature is given by $\sinh(2\beta_{c}^{\Box hv}J_{v})\sinh(2\beta_{c}^{\Box hv}J_{h})=1$ Kramers and Wannier (1941), which serves as a better estimation of the lower bound of $T_{c}^{G}$, i.e., $T_{c}^{\Box hv}\lesssim T_{c}^{G}$.

The critical temperatures of the square lattice $k_{B}T_{c}^{\Box}=\frac{2J}{\ln(1+\sqrt{2})}\approx 2.2692J$ for isotropic case and $k_{B}T_{c}^{\Box hv}\approx 3.2820J$ for anisotropic case , and that of the triangular lattice $k_{B}T_{c}^{\Delta}=\frac{4J}{\ln 3}\approx 3.6410J$. Fig. 1 shows the Monte Carlo simulation of the glacia transition in the high density limit. The critical temperature $k_{B}T_{c}^{G}\approx(3.5\pm 0.1)J$.

So far we have shown that the DDL has a glacia phase when the dimer density is high and the temperature is low. A particularly interesting limit is the high density and low temperature limit, where the dimers completely cover the lattice $L$ and the charges are ordered. In this limit, as the charge degrees of freedom are frozen, it reduces to the dimer-covering problem which is known to be exactly solvable Lieb (1967); Kasteleyn (1961, 1963); Fisher (1961); Temperley and Fisher (1961). Any dimer-covering configuration admit two ground states, all dimers are of either $A$-type or $B$-type. The ground states are massively degenerate due to the exponentially many possible ways of covering the lattice $L$. In the thermodynamic limit, the number of arrangements per dimer is $\mathrm{e}^{2G/\pi}$, where $G\approx 0.91596$ is the Catalan constant. For a given lattice $L$, we define the ground state configurational entropy (not including charge degrees of freedom)

[TABLE]

where $W$ is the number of the ways that dimers completely cover $L$.

Configurational entropy in the presence of the polymers.– We are also able to calculate the configurational entropy in the presence of the polymers. For polymers with fixed conformations, the configurational entropy is the entropy of the surrounding solvent molecules, which should be distinguished from the conformational entropy of the polymer. In the thermodynamic equilibrium, it is the total entropy of the system that is maximized, which should include both the conformational entropy of the polymers and the configurational entropy of the solvent.

When the density of the polymers is low, we may consider a system consists of only one polymer surrounded by the DDL. By definition, a polymer is a string of dimers $D^{p}=(d^{p}_{1},d_{2}^{p},\dots,d_{N}^{p})$. We call the integer $N$ the length of the polymer denoted as $\lambda_{p}=N$, and we call these dimers constituent dimers. The constituent dimers do not overlap with each other and any two consecutive dimers $(d^{p}_{i},d^{p}_{i+1})$ are neighbors. The constituent dimers may also carry charges and dipole moments, as the hydrophilic residues of the amino acids or nucleobases of the nucleotides. A polymer may form loops, where two dimers that are not consecutive in the series $D^{p}$ sit on neighboring sites and be joint together with a chemical bond, e.g., the disulfide bond in the proteins. For simplicity, we shall constraint ourselves to the neutral polymers.

Suppose we fix a conformation $\Omega_{p}$ of the polymer of length $\lambda_{p}=N$, and put it in the DDL. In the high density and low temperature limit, we may calculate the configurational entropy of the DDL

[TABLE]

where “$|\Omega_{p}$” means in the presence of a polymer in the given conformation $\Omega_{p}$, $W(L|\Omega_{p})$ is the number of the ways that the dimers completely covers $L$ in the presence of the polymer, and $N_{D}$ counts the number of the connected domains, into which the polymer cuts $L$. The loss of the entropy due to the presence of the polymer is measured by

[TABLE]

The limit ${L\rightarrow\infty}$ means the scale of the lattice is much larger than the scale of the polymer, then $S_{\text{Config}}(\Omega_{p})$ should be independent of the lattice but only the conformation of the polymer. By abuse of language, we call $S_{\text{Config}}(\Omega_{p})$ configurational entropy of the polymer for the conformation $\Omega_{p}$. $W(L|\Omega_{p})$ is counted in a modified Kasteleyn’s method Jacobsen (2016); de Tiliere .

Naïvely, one would expect $S_{\text{Config}}(\Omega_{p})$ is proportional to the length of the polymer $\lambda_{p}$. We considered three kinds of conformations: 1) closely packed ($\Omega_{p}=\|$), 2) double-strand ladder ($\Omega_{p}=\#$) and 3) bubble shaped ($\Omega_{p}=O$). It turns out that this is not always the case. Fig. 2 shows the configurational entropies $S_{\text{Config}}(\|)$, $S_{\text{Config}}(\#)-S_{\text{Config}}(\|)$ and $S_{\text{Config}}(O)-S_{\text{Config}}(\|)$ in log-log scale. For closely packed polymer, the configurational entropy scales linearly with respect to the polymer length $\lambda_{p}$, which agrees with the naïve expectation. For double-strand ladder polymer, $S_{\text{Config}}(\#)-S_{\text{Config}}(\|)\propto\lambda_{p}^{2}$ asymptotically. For bubble shaped polymer, $S_{\text{Config}}(O)-S_{\text{Config}}(\|)\propto\lambda_{p}^{\alpha}$ with $\alpha=1.6\pm 0.1$. 777Ref. 10 studied random domino tiling of Aztec diamond, and found a similar effect for the entropy.

If one examines the microscopic configurations of the dimers, one find that the polymer does not only exclude a certain volume of the liquid, but also aligns partially the dimers around it. This seems to be consistent with the empirical picture for the hydrophobic effect in chemistry. The water molecules form the “cages” around the polymers. It is important to note that, since for conformations not closely packed, $S_{\text{Config}}(\Omega_{p})\sim\lambda_{p}^{\alpha}$ with $\alpha>1$, the configurational entropy growths rapidly with the scale of the polymer. In the protein folding problem, it implies that the conformations of the proteins might not be properly understood if the surrounding water molecules were neglected.

Conclusion.–We studied the glacia phase transition of the DDL by mapping it to the annealed Ising model on random graphs. In the high density limit the critical temperature is estimated to be $k_{B}T_{c}^{G}=(3.5\pm 0.1)J$. We also calculated the configurational entropy of the DDL in the presence of the polymers and found a non-trivial scaling of the suppression of the configurational entropy. The DDL may serve as a simple model for understanding the complex phenomena in polar solvents.