Bond-counting potentials -- A classical many-body model of covalent bonding with exact solutions in one dimension

TL;DR

This paper introduces bond-counting potentials as a simple yet effective model for covalent bonding, providing exact solutions in one dimension to study phase transitions like LLPT in elemental melts.

Contribution

It presents the first analytic solution of a one-dimensional bond-counting potential system, capturing phase transitions and structural properties relevant to covalent materials.

Findings

System with valency 1 shows molecular dominance at low T.

Valency 2 system exhibits first-order LLPT from molecular to polymeric phase.

Finite systems show phase separation and non-equilibrium transformations.

Abstract

We introduce "bond-counting" potentials, which provide an elementary description of covalent bonding. These simplistic potentials are intended for studies of the mechanisms behind a variety of phase transitions in elemental melts, including the liquid-liquid phase transitions (LLPT) in phosphorus and bismuth. As a first study employing such potentials, an analytic solution of a one-dimensional model system is presented, including its thermodynamic properties and its structure factor. In the simplest case, the chemical valency of each atom is $1$, and either single atoms or diatomic molecules are present. At low temperatures and moderate pressures, the system consists almost exclusively of molecules, and single atoms act as topological defects. A slightly more complicated case involves a valency of $2$, with either single or double bonding. This system exhibits a first-order LLPT from a molecular to a polymeric phase, as in phosphorus. In this case, the one-dimensional model system exhibits phase separation for finite-sized systems at low temperatures. A variant of this system also exhibits a non-equilibrium phase transformation upon heating the molecular condensed phase, qualitatively similar to boiling in white phosphorus.

Equations

$$q_i=\sum _{k\ne i}\theta (a-r_{ik}),$$

$$\mathcal{V}(r_{ij},q_i,q_j)=\{\begin{array}{ll}v(r_{ij})& q_i,q_j\le z\\ \infty & q_i>z\mathrm{or}{q}_j>z\end{array},$$

$$v(r)=\{\begin{array}{ll}\infty & r<d\\ v_{\mathrm{bond}}(r)& d<r<a\\ v_{\mathrm{non}-\mathrm{bond}}(r)& r>a\end{array}.$$

$$q_i=\sum _{k\ne i}\left[\theta (a-r_{ik})+\theta (\tilde{a}-r_{ik})\right],$$

$$v(r)=\{\begin{array}{ll}\infty & r<d\\ v_{\mathrm{double}-\mathrm{bond}}(r)& d<r<\tilde{a}\\ v_{\mathrm{single}-\mathrm{bond}}(r)& \tilde{a}<r<a\\ v_{\mathrm{non}-\mathrm{bond}}(r)& r>a\end{array}.$$

$${\mathcal{V}}_{\mathrm{EAM}}({\mathbf{r}}_1,$$

$$\sum _{i=1}^{N}F\left(\sum _{j\ne i}\rho (r_{ij})\right)+\sum _{i=2}^N\sum _{j=1}^{i-1}U(r_{ij}),$$

$$F(x)=\{\begin{array}{ll}0& x\le z\\ \infty & x>z\end{array},$$

$$U(r)=v(r).$$

$$\mathcal{H}(\mathbf{x},\mathbf{p})=\frac{1}{2m}\sum _{i=1}^N{p}_i^2+\sum _{⟨i,j⟩}u(|x_i-x_j|).$$

$$Z_N(L,T)$$

$$=\frac{1}{N!{\lambda }_T^N}\int d^{N}x{e}^{-\beta {\sum }_{⟨i,j⟩}u(|x_i-x_j|)},$$

$$Z_N(L,T)=\frac{1}{{\lambda }_T^N}{\int }_{\mathcal{D}}{d}^{N}x{e}^{-\beta {\sum }_{i=1}^{N-1}u(x_{i+1}-x_i)}.$$

$$Z_N(P,T)=\frac{1}{l_0}{\int }_0^{\infty }dL{Z}_N(L,T)e^{-\beta PL},$$

$$Z_N(P,T)=\frac{1}{l_0{\lambda }_T^N}{\int }_0^{\infty }d{y}_{1}e$$

$${\left({\int }_0^{\infty }dy{e}^{-\beta \left[u(y)+Py\right]}\right)}^{N-1}$$

$$g(P,T)$$

$$=-\frac{1}{\beta }\ln \left(\frac{1}{{\lambda }_T}{\int }_0^{\infty }dy{e}^{-\beta \left[u(y)+Py\right]}\right).$$

$$u_{\mathrm{HS}}(x)=\{\begin{array}{ll}\infty & x<a\\ 0& x>a\end{array},$$

$$Z_N^{(\mathrm{HS})}(P,T)=\frac{1}{\beta P{\lambda }_T}{\left(\frac{e^{-\beta Pa}}{\beta P{\lambda }_T}\right)}^{N-1}.$$

$$g_{\mathrm{HS}}(P,T)=Pa+\frac{1}{\beta }\ln \left(\beta P{\lambda }_T\right).$$

$$A={\int }_d^{a}dy{e}^{-\beta \left[v_{\mathrm{bond}}(y)+Py\right]},$$

$$B={\int }_a^{\infty }dy{e}^{-\beta \left[v_{\mathrm{non}-\mathrm{bond}}(y)+Py\right]}.$$

$$Z_N(P,T)=\frac{1}{{\lambda }_T^N}\sum _{\{{\sigma }_i\}}\prod _{i=1}^N{T}_{{\sigma }_i{\sigma }_{i+1}}=\frac{1}{{\lambda }_T^N}\mathrm{Tr}\left[T^N\right],$$

$$T=\left(\begin{array}{cc}B& B\\ A& 0\end{array}\right).$$

$$g(P,T)$$

$$=g_{\mathrm{Tak}}(P,T)-\frac{1}{\beta }\ln \left[\frac{1}{2}\left(1+\sqrt{1+\frac{4A}{B}}\right)\right],$$

$$u(y)=\{\begin{array}{ll}\infty & y<a\\ v_{\mathrm{non}-\mathrm{bond}}(y)& y>a\end{array}.$$

$$\ell ={\left(\frac{\partial g}{\partial P}\right)}_T={\ell }_{\mathrm{Tak}}-\frac{\frac{2}{\beta }\frac{\partial }{\partial P}{\left(\frac{A}{B}\right)}_T}{\left(1+\sqrt{1+\frac{4A}{B}}\right)\sqrt{1+\frac{4A}{B}}},$$

$$s=-{\left(\frac{\partial g}{\partial T}\right)}_P=s_{\mathrm{Tak}}+k_{\mathrm{B}}\ln \left[\frac{1}{2}\left(1+\sqrt{1+\frac{4A}{B}}\right)\right]+$$

$$\frac{\frac{2}{\beta }\frac{\partial }{\partial T}{\left(\frac{A}{B}\right)}_P}{\left(1+\sqrt{1+\frac{4A}{B}}\right)\sqrt{1+\frac{4A}{B}}},$$

$${\ell }_{\mathrm{Tak}}={\ell }_{\mathrm{HS}}=a+\frac{1}{\beta P},$$

$$s_{\mathrm{Tak}}=s_{\mathrm{HS}}=k_{\mathrm{B}}\left[\frac{3}{2}-\ln \left(\beta P{\lambda }_T\right)\right].$$

$$\frac{A}{B}=e^{\beta ϵ}\left[e^{\beta P\left(a-d\right)}-1\right].$$

$$\ell =\frac{a+d}{2}+\frac{1}{\beta P}+(a-d)\mathcal{O}(e^{-\beta P\left(a-d\right)}),$$

Full text

I Introduction

Many liquids and solids are well-described within the van der Waals picture: their structure and dynamics are determined primarily by strong repulsive forces that give rise to severe constraints on the particles’ motion, whereas the overall density is largely determined by weaker long-range attractive forces. Hansen&McDonald ; CD83 The packing constraints dominate the short-range order characteristic of the liquid state, and lead to a quasi-universal structure of simple liquid metals and noble gases, at least in the vicinity of the melting curve, with a structure close to that of a system of hard spheres, and coordination numbers of $z=9-11$. Hansen&McDonald ; ANW66 ; YJL73 ; DJC16

On the other hand, covalent-bonded elements with open crystalline structures — mainly group IV, V, and VI elements — exhibit liquid structures characterized by low coordination numbers. The fact that such a liquid structure may be highly non-trivial is reflected by intricate phase diagrams, which include thermodynamic anomalies and/or liquid-liquid phase transitions (LLPTs). BVV97 ; BVV02 ; BVV03 ; MPF07 Indeed, a growing number of LLPTs have been observed in elemental liquids during the last three decades, e.g., in phosphorus, KY00 ; KY02 ; MG03 selenium, BVV89 ; BVV91 sulfur, BVV91 ; LL14 ; ZG14 ; HL18 , bismuth UAG92 ; GY09 ; SY17 ; EM18 and tellurium. TY85 ; BVV92

In some cases, there exists a general understanding of the mechanism behind an LLPT. The most well-known example is phosphorus, which at temperatures of about $1300\,$K and moderate pressures exhibits a fluid state consisting of pyramid-shaped $P_{4}$ molecules. As the pressure is increased to about $1\,$GPa, the volume occupied by the pyramidal molecules becomes prohibitive, and the liquid collapses to a much denser polymeric structure. KY00 ; KY02 This LLPT is first-order and reversible, with a coexistence line ranging from $(P,T)\approx(0.9\,\mathrm{GPa},1300\,\mathrm{K})$ to $(P,T)\approx(0.3\,\mathrm{GPa},2450\,\mathrm{K})$. MG03

A further consequence of strong covalent bonds is the observation of non-equilibrium phase transitions between metastable states. BVV06 A prominent example of an element with such states is, again, phosphorus. Phase Diagrams The $P_{4}$ molecules can be condensed from the vapor, and the resulting material is known as white phosphorus. White phosphorus exists in several different molecular phases as a function of temperature, and undergoes solid-solid, melting and boiling transitions between these metastable phases, despite the fact that the thermodynamically stable phase under these conditions is the non-molecular black phosphorus state. In the vicinity of room temperature, the lifetimes of the metastable phases are exceedingly long, but near boiling, at $\approx 550\,\mathrm{K}$, polymerization of phosphorus begins to occur, DTW46 and the original $P_{4}$ molecules do not reform upon subsequent cooling to room temperature conditions. The present description of the phase behavior of phosphorus does not do justice to the complexity of the topic (e.g., we have not even mentioned red phosphorus). Instead, it focuses exclusively on those aspects — the LLPT and the existence of a non-equilibrium boiling transformation — which, as we will see below, can be qualitatively reproduced within a simplistic model for covalent bonding in one dimension (1D).

In other cases, the mechanisms behind observed LLPTs have not yet been identified. For example, different types of evidence for several LLPTs in bismuth have been reported, both at high pressures UAG92 ; EM18 and at ambient pressure. LXF07 ; GY09 A recent set of experiments presents further challenges for an understanding of the behavior of this elemental liquid.SY17 ; SY16a ; SY16b In these studies, bismuth is compressed to $\sim 2\,$GPa and heated to $\sim 2000\,$K, and then cooled and decompressed back to ambient conditions. The resulting solid contains structural defects and displays a variety of anomalies, including the following peculiarity at melting: upon heating at ambient pressure, this material undergoes a transformation similar to standard melting of Bi I (as evidenced by an endothermic peak which has the appropriate magnitude and temperature, and is largely reversible), but the material does not flow and instead retains its shape like a solid [see the insets in Fig. 1 of Ref. SY16b, ]. Identifying the mechanism behind such a phenomenon is a challenge which requires significant experimental and theoretical work. We will present some speculative ideas on this in Sec. IV, and explain in what sense the present work may serve as a preliminary step in a 3D-simulation approach to this problem.

Several theoretical approaches have been used to investigate LLPTs. On the phenomenological level, the two-state model SS65 ; RE67 was used to give a simple explanation of LLPTs, for both liquids and amorphous solids. BVV97 ; TH00 ; PEG03 According to this model, an LLPT is a result of competition between two kinds of clusters which differ in their short-range order. One liquid phase consists predominantly of entropically favorable clusters of high density, and the second phase consists of a large concentration of energetically favorable clusters of low density (i.e., with a more open structure). twostatevariants The model assumes an energy cost $J$ for mixing the two types of clusters, and yields a first-order phase transition line, terminating at a critical point at a temperature of $J/(2k_{\mathrm{B}})$. At lower temperatures, the two types of clusters exhibit phase separation despite being made of the same substance (i.e., exhibit a miscibility gap). This type of modeling can describe the thermodynamics of the system without requiring a microscopic understanding of the mechanisms of formation and interaction of the two distinct clusters involved.

More elaborate theoretical analyses employ molecular dynamics simulations, either combined with ab-initio density functional theory calculations MT01 ; GM05 or based on empirical potentials. PHP03 ; SS03 Most of these simulations strive to achieve a quantitative understanding of the phenomena in specific materials, and are often successful, but are limited by computational power for the first type of simulation, and by the accuracy of the potentials for the latter type. For example, the predictions of LLPTs by such simulations sometimes contradict experimental observations [see, e.g., the case of nitrogen in Ref. WG17, ], indicating possible gaps in our understanding. A complementary approach in seeking a better understanding of LLPTs, is based on simplistic potentials with very few parameters, where the goal is to reproduce observed phenomena qualitatively, rather than quantitatively.

The canonical example of such a simplistic approach to the liquid state of matter is the hard sphere potential, with only one parameter (the spheres’ diameter), which provides an elementary prototype of the melting transition. ABJ57 ; WWW57 ; HWG68 More recently, a variety of simple soft-core isotropic pair potentials have been investigated extensively, and have been shown to exhibit a wide range of non-trivial phenomena, including LLPTs, polyamorphism in glasses, anomalous melting and water-like liquid anomalies [see Ref. BSV09, and references therein]. In some cases, such as the density maximum in water, the phenomenon of interest could be qualitatively reproduced even within a simplified 1D model, by analytic solution of the thermodynamics of the system. CCH96 ; SLMR99 ; BNA08

In the present work we introduce “bond-counting” many-body potentials, in order to provide a simple description of covalent bonding. In the simple cases studied in detail, each atom is characterized by two size parameters — a “core diameter” which is impenetrable, and a larger “bonding zone diameter” which may be penetrated by up to $z$ neighbors, each forming a bond with a binding energy of one unit. The case $z=1$ corresponds to atoms which form diatomic molecules, such as hydrogen or fluorine; $z=2$ describes atoms which can polymerize, such as sulfur, and $z=3$ may describe bismuth and phosphorous (in 3D). The present contribution is limited to 1D, where analytical solutions are possible. The thermodynamics and structure of the simplest model with $z=1$ are analyzed in detail. Subsequently, a qualitative description of the above-mentioned phenomena observed in phosphorus is obtained for $z=2$. This is done by allowing for double bonds between identical atoms, which opens the possibility of competing molecular and polymeric configurations even in a very simple 1D system.

The paper is organized as follows. Sec. II provides the general definition of potentials which count the number of bonds. The exact solution of this model in 1D is discussed in Sec. III. In Sec. IV we summarize the results and discuss possible extensions of this work, including higher dimensions.

II Bond-counting model potentials for covalent bonding

Before providing the definition of a new family of potentials, it is appropriate to mention a few more general facts regarding existing empirical potentials. Hansen&McDonald First, it is noteworthy that both qualitative and remarkable quantitative accuracy in the description of simple liquids, consisting of atoms with van der Waals, metallic or ionic interactions can be achieved using simple pair potentials, such as the Lennard-Jones or (possibly screened) Coulomb potentials. For molecular liquids, models which combine pair potentials between atoms in different molecules with rigid connections within each molecule are also successful, as long as breaking and reformation of the covalent bonds is not relevant. Of the very many multi-parameter many-body empirical potentials developed to improve the quantitative accuracy of simulations, we will mention only two: the embedded-atom method DMS93 (EAM), to which we will return below, and reactive force fields, STP16 which were developed expressly in order to allow simulations of covalent bond breaking and formation.

Liquids of the latter type, for which the bonds within each molecule are evolving (as in polymerization), are referred to as “associating liquids”. A simple model for these describes each atom or monomer as a hard sphere, augmented by several much smaller off-center spheres with attractive square-well potentials. WMS84 ; JG88 If the attractive potential is sufficiently short-ranged, the formation of more than one bond at each off-center site is excluded. Approximations appropriate for such models have been developed, WMS84 and have seen considerable success in the framework of the statistical associating fluid theory (SAFT). CWG The bond-counting potentials introduced here express the idea of chemical valency directly, rather than relying on geometric constraints and pair potentials. In particular, the extension to double or triple bonds is straightforward within the bond-counting approach, see below.

In order to define bond-counting potentials, consider a system of $N$ classical particles (”atoms”) of mass $m$ and positions $\mathbf{r}_{i}$. Distinguishing distances shorter than the bond-range $a$ from longer distances, each atom is associated not only with its position, but also with the number of its bonds

[TABLE]

where $r_{ik}=|\mathbf{r}_{i}-\mathbf{r}_{k}|$ is the distance between atoms $i$ and $k$ and $\theta(x)$ is the step function. The valency $z$ of each atom cannot be exceeded, namely $q_{i}\leq z$, making this a many-body potential. In addition, consider a potential $\mathcal{V}_{ij}$ between each pair of atoms $i$ and $j$ which may, in general, depend on the number of bonds of both, namely $\mathcal{V}_{ij}=\mathcal{V}(r_{ij},q_{i},q_{j})$. In the simple cases to be studied here, the potential is taken to be independent of $q_{i}$ and $q_{j}$, provided that both are not larger than $z$:

[TABLE]

where

[TABLE]

Here, $d$ is an impenetrable core diameter (which is optional in the sense that $d=0$ is allowed), and $v_{\mathrm{bond}}(r)$, $v_{\mathrm{non-bond}}(r)$ are the interaction potentials inside and outside the bonding zone, respectively.

Double bonds can be taken into account in a straightforward manner within this approach, by introducing an additional diameter $\tilde{a}<a$, such that the region $d<r<\tilde{a}$ corresponds to a double bond. The number of bonds of atom $i$ is then

[TABLE]

and the potential $v(r)$ in Eq. (2) is replaced by

[TABLE]

Similarly, triple bonds can be described by introducing yet another diameter.

We note that just as hard spheres can be viewed as a limiting form of pair potentials, the bond-counting potentials can be represented as limiting forms of more elaborate many-body potentials. For example, using the EAM with the notation

[TABLE]

one obtains the above potentials by using the choice $\rho(r_{ij})=\theta(a-r_{ij})$ for a model with single bonds ($\rho(r_{ij})=\theta(a-r_{ij})+\theta(\tilde{a}-r_{ij})$ for a model with double bonds), with the embedding energy $F$ and pair potential $U$ given by

[TABLE]

Thus, simulation software with an EAM option can be used to obtain results for the present family of potentials, provided that numerical issues associated with the use of step functions are resolved. Note that this limiting case is unconventional in the sense that the embedding energy is purely repulsive, with all the attraction in the pair potential. This is the opposite of the customary use of the EAM, which was originally developed for metals.

It should be emphasized that the family of potentials defined above is isotropic — there is no explicit dependence on bond angles. As appropriate for liquids, it is suited for identifying mechanisms which rely exclusively on the notion of chemical valency, and are not sensitive to the details of the geometry. Adding angular dependencies, as in the development of the modified EAM, Baskes92 is a natural possibility for future studies.

III Exact solution in 1D

In this section we briefly review the Takahashi solution for a classical system with nearest-neighbor interactions in 1D (Sec. III.1), and then use the transfer matrix method to generalize and apply the solution to two models of bond-counting potentials: (1) a model with $z=1$, for which we study in detail the equilibrium statistical mechanics and the resulting structure (Sec. III.2); (2) a model with $z=2$ which allows for double bonds, exhibiting a polymerization transition (Sec. III.3).

III.1 The Takahashi solution for 1D systems with nearest-neighbor interactions

We summarize the Takahashi solution TH42 ; Lieb&Mattis for a 1D classical system of $N$ particles of mass $m$, with Hamiltonian of the form

[TABLE]

Here, $x_{i}$ and $p_{i}$ are the position and momentum of the $i$th particle, and $u(x)$ is the nearest-neighbor interaction potential (the sum in the last term is over nearest neighbors $i$, $j$). Comment1

The canonical partition function for a system of length $L$ at temperature $T$ reads

[TABLE]

where $\beta=1/k_{\mathrm{B}}T$ and $\lambda_{T}=h/\sqrt{2\pi mk_{\mathrm{B}}T}$, with $k_{\mathrm{B}}$ and $h$ the Boltzmann and Planck constants. The crucial step in deriving the exact solution is that in 1D the last integral in (III.1) is $N!$ times the integral over the domain $\mathcal{D}$, defined by $0\leq x_{1}\leq x_{2}\leq\ldots\leq x_{N}\leq L$. The canonical partition function then takes the form

[TABLE]

The evaluation of the partition function is simple in the isobaric-isothermal ($PTN$) ensemble, for which

[TABLE]

where $l_{0}$ is a basic unit of length required to render the partition function $Z_{N}(P,T)$ dimensionless. In the thermodynamic limit the choice of $l_{0}$ is not important; below $l_{0}=\left(\beta P\right)^{-1}$ is used, so that the condition $Z_{N=0}(L,T)=1$ implies $Z_{N=0}(P,T)=1$. Substituting Eq. (10) into (11) and changing variables to $y_{1}=x_{1}$, $y_{i}=x_{i}-x_{i-1}$ for $i=2,\ldots,N$ and $y_{N+1}=L-x_{N}$ (such that $\sum^{N+1}_{i=1}y_{i}=L$), one obtains

[TABLE]

All thermodynamic properties can then be calculated exactly from the Gibbs free energy per particle (or chemical potential), which in the thermodynamic limit is

[TABLE]

An example, to be used below as a reference, is the hard-sphere model

[TABLE]

for which Eq. (III.1) yields

[TABLE]

The corresponding Gibbs free energy per particle [Eq. (III.1)] is

[TABLE]

III.2 Bond-counting potential with $\mathbf{z=1}$

Using the transfer matrix method, the Takahashi solution may be generalized to solve for a system with a $z=1$ bond-counting potential in 1D. The relevant potentials for $z=1$ are given by Eqs. (1)-(3). Each atom can have at most a single bond with one of its nearest neighbors, and the system consists of a mixture of single atoms and diatomic molecules (i.e., paired atoms separated by a distance $d<y<a$). We note that a similar method was used in Ref. BNA92, to solve a 1D model involving two types of atoms (or rather, two orientations of molecules). In the solution below the index of the transfer matrix specifies the type of bond, rather than the type of atom as in Ref. BNA92, .

III.2.1 Thermodynamics

In the $PTN$ ensemble each molecular bond will contribute to the partition function a factor

[TABLE]

whereas each pair of adjacent atoms which are not bound to each other will contribute a factor

[TABLE]

The partition function for $N$ atoms is then $l_{0}\lambda^{-N}_{T}$ times a sum of all possible products of $N-1$ factors of this type (such as $ABABBA$ for $N=7$), where the only restriction on the products is that two consecutive bond factors $A$ cannot occur. An appropriate recursion relation is available, $Z_{N}=\lambda^{-1}_{T}B\,Z_{N-1}+\lambda^{-2}_{T}AB\,Z_{N-2}$, expressing the fact that the last atom in the chain can be either unbound or bonded to the preceding atom: one can add to shorter chains either a single atom or a molecule. This relation is similar to the recursion relation for the Fibonacci sequence, providing a hint that expressions in terms of powers of a $2\times 2$ matrix may be possible.

The calculation of the partition function can be simplified by closing the 1D chain, such that the $N$th atom becomes a neighbor of the first atom. One may define the bond variables $\{\sigma_{i}\}^{N}_{i=1}$, where $\sigma_{i}=1$ if atoms $i$ and $i+1$ are bonded (i.e., $d<y_{i+1}<a$) and $\sigma_{i}=0$ otherwise. With this definition, the contribution of the $i$th bond to the partition function is $B$ if $\sigma_{i}=0$, and $A$ if $\sigma_{i}=1$. The restriction that two adjacent bonds are not allowed implies that the contribution vanishes if $\sigma_{i}=\sigma_{i+1}=1$ for some value of $i$. The partition function can thus be written in the form

[TABLE]

where $\sigma_{N+1}\equiv\sigma_{1}$ and the transfer matrix $T$ is

[TABLE]

The eigenvalues of this transfer matrix are $\Lambda_{\pm}=\left(B\pm\sqrt{B^{2}+4AB}\right)/2$ and since $\Lambda_{+}>\Lambda_{-}$, in the thermodynamic limit the Gibbs free energy per particle reads

[TABLE]

where $g_{\mathrm{Tak}}(P,T)=-k_{\mathrm{B}}T\ln\left(B/\lambda_{T}\right)$ is the Gibbs free energy of the Takahashi model [Eq. (III.1)] with an interaction potential $u(y)$ between a paired and a single nearest-neighbor atoms,

[TABLE]

The length $\ell(P,T)$ and entropy $s(P,T)$ per particle are therefore

[TABLE]

where $\ell_{\mathrm{Tak}}=\left(\partial g_{\mathrm{Tak}}/\partial P\right)_{T}$ and $s_{\mathrm{Tak}}=-\left(\partial g_{\mathrm{Tak}}/\partial T\right)_{P}$.

As a simple example, let us consider the case $v_{\mathrm{non-bond}}(x)=0$ and $v_{\mathrm{bond}}(x)=-\varepsilon$, where $\varepsilon>0$ is a constant energy. Then $u(y)$ [Eq. (22)] reduces to the hard-sphere potential $u_{\mathrm{HS}}(y)$ [Eq. (14)], and thus $g_{\mathrm{Tak}}=g_{\mathrm{HS}}$. From Eq. (16) we obtain

[TABLE]

Note that Eq. (26) is also the entropy per particle of a 1D ideal gas at fixed pressure $P$ and temperature $T$, and that the hard-sphere diameter in Eq. (25) is the bond range $a$ rather than the core diameter $d$.

The ratio $A/B$ for this system is

[TABLE]

In the high density (i.e., low temperature or high pressure) limit $\beta P\left(a-d\right)\gg 1$ Eqs. (23)-(27) give

[TABLE]

This limit corresponds to a close-packed arrangement of $N/2$ molecules, each of diameter $a+d$. The system then behaves effectively as $N$ hard spheres each of diameter $\left(a+d\right)/2$. In the low density and low temperature limit, $e^{-\beta\varepsilon}\ll\beta P\left(a-d\right)\ll 1$, we obtain

[TABLE]

In this limit the length of each molecule follows the equation of state (25) of hard spheres with diameter $\left(3a+d\right)/2$, corresponding to the average diameter of each molecule. In the low density and high temperature limit, $\beta P(a-d)\ll e^{-\beta\varepsilon}$, we find

[TABLE]

which, if the temperature is also high in the sense that $\beta\varepsilon\ll 1$, corresponds to a system of $N$ hard spheres of diameter $d$.

The length per atom $\ell(P,T)$ is plotted in Fig. 1(a) as a function of pressure at various temperatures, for $d/a=0.5$. The various limits discussed above are more easily observed in a plot of the length difference, $\Delta\ell(P,T)=\ell(P,T)-\ell_{\mathrm{HS}}(P,T)$ as a function of $\beta Pa$. Figure 1(b) shows such a plot for the same data as in Fig. 1(a). For $\beta P\left(a-d\right)\gg 1$ all curves approach the value $\Delta\ell/a=\left(d/a-1\right)/2=-0.25$ independent of the value of $\beta\varepsilon$, in agreement with Eq. (III.2.1). For $e^{-\beta\varepsilon}\ll\beta P\left(a-d\right)\ll 1$ Eq. (29) gives $\Delta\ell/a\approx-\left(2\beta Pa\right)^{-1}$ which appears as a straight line on a log-log scale. This behavior is demonstrated by the black curve in Fig. 1(b). Finally, in the limit $\beta P(a-d)\ll e^{-\beta\varepsilon}$ the correction to Eq. (III.2.1) of order $\beta\varepsilon$ gives $\Delta\ell/a\approx-e^{\beta\varepsilon}\left(1-d/a\right)$. This is the value approached by the curves in Fig. 1(b) for $\beta P(a-d)\ll e^{-\beta\varepsilon}$.

Figure 2 shows the entropy difference per atom, $\Delta s(P,T)=s(P,T)-s_{\mathrm{HS}}(P,T)$, as a function of $\beta Pa$. The entropy difference approaches zero in the limits $\beta P\left(a-d\right)\gg 1$ or $\beta P(a-d)\ll e^{-\beta\varepsilon}$, in accordance with Eqs. (III.2.1) and (III.2.1). It is also interesting to note that for $P=\left(a-d\right)/\varepsilon$ one has $A/B=e^{\beta\varepsilon}\left(e^{\beta\varepsilon}-1\right)$, and $\Delta g=g-g_{\mathrm{HS}}$ [Eq. (III.2.1)] becomes independent of temperature. Thus $\Delta s=0$ for $P=\left(a-d\right)/\varepsilon$, irrespective of the temperature.

Another quantity of interest is the fraction of molecules

[TABLE]

This quantity approaches 1 in the limit $\beta P\left(a-d\right)\gg e^{-\beta\varepsilon}$, and approaches 0 for $\beta P\left(a-d\right)\ll e^{-\beta\varepsilon}$, as shown in Fig. 3. The first case contains the two limits $\beta P\left(a-d\right)\gg 1$ or $e^{-\beta\varepsilon}\ll\beta P\left(a-d\right)\ll 1$ discussed in Eqs. (III.2.1) and (29), in which the system consists of $N/2$ molecules.

III.2.2 The structure factor

The transfer matrix method also enables the calculation of the structure factor of the system,

[TABLE]

The inter-atomic distance can be written as a sum over bond lengths, $x_{i}-x_{j}=\sum^{i}_{k=j+1}y_{k}$, and therefore one can factor the oscillatory term here, and use a variant of the transfer matrix method to calculate this quantity. Adding the appropriate oscillatory factor to the expressions for $A$ and $B$ results in:

[TABLE]

together with an analogue of the transfer matrix:

[TABLE]

In terms of these quantities, one obtains

[TABLE]

The eigenvalues of the matrix $t$ are $\lambda_{\pm}=\left(B_{q}\pm\sqrt{B^{2}_{q}+4A_{q}B_{q}}\right)/2$. Substituting the similarity transformations $T=MWM^{-1}$ and $t=mwm^{-1}$, where

[TABLE]

we obtain from Eq. (III.2.2)

[TABLE]

In the thermodynamic limit $N\rightarrow\infty$ only the mode with eigenvalue $\Lambda_{+}$ survives, and one ends up with

[TABLE]

where

[TABLE]

With the form (III.2.2), the limit $N\rightarrow\infty$ of Eq. (III.2.2) gives

[TABLE]

The corresponding radial distribution function $g(x)$ is given by the Fourier transform

[TABLE]

where $n=\ell^{-1}$ is the particle density.

The structure factor and radial distribution function for the case $d/a=0.5$, $v_{\mathrm{non-bond}}(x)=0$ and $v_{\mathrm{bond}}(x)=-\varepsilon$ discussed above are shown in Fig. 4, for a constant density $na=0.75$ (note that at a constant density the values of $\beta Pa$, which are given in the figure, vary only to a limited extent; a plot of the structure factors at a fixed value $\beta Pa\approx 1$ will be qualitatively similar). At low temperatures the structure factor shown in the upper panel has a pronounced first peak, which decreases in height and develops into a shoulder-like feature as the temperature is increased. The radial distribution function displayed in the lower panel has two main features corresponding to neighbors at distance $x=d$ and $x=a$, with discontinuities reflecting the discontinuities of the potential. At high temperatures these two features are distinct peaks, and the height of the second one is larger, since molecules are dissociated into individual atoms. At lower temperatures, the number of bonded molecules increases and the height of the first peak becomes dominant. Interestingly, at low enough temperatures the discontinuity at $x=a$ changes sign, and there is no longer a peak at this value of $x$. The evolution of this discontinuity with temperature is reflected in the strong oscillations observed in the upper panel for large $q$.

III.2.3 Phases and boundaries

The above, especially Fig. 3, shows how this system transforms continuously between a purely atomic fluid and a molecular fluid. It is intuitively clear that a qualitatively similar crossover will occur also in 2D and 3D. In particular, in this 1D model the $P\to 0$ and the $T\to 0$ limits are different and do not commute, and there is in this sense a discontinuity at $T=P=0$. A similar discontinuity has been identified in previous studies of 1D systems as representing the liquid-gas critical point of 3D fluids (see, e.g., Sec. VII.A of Ref. SLMR99, ). It is worth noting that for a system defined just as above but with a negative value of $\varepsilon$, this critical point occurs at a finite pressure, $P_{\mathrm{C}}=|\varepsilon|/(a-d)$, as follows from Eq. (III.2.1). In this case, at low pressure the atoms repel each other and are separated by distances larger than $a$, but a pressure $P>P_{\mathrm{C}}$ causes pairs of atoms to overcome their mutual repulsion. This transition is similar to that found in a system of particles interacting through soft-core pair potentials, see Refs. BSV09, and SLMR99, , and references therein.

Returning to the present system with a positive binding energy, one can identify additional features related to phase behavior. If the temperature is low enough, the fraction of molecules approaches 100%, and only very few unpaired atoms remain. The system then consists of large domains of molecules, separated by “boundary defects” which are unbound atoms, and which can only be created or annihilated in pairs. In this sense, the single atoms act as topological defects in the system.

The dynamics of these defects may be of interest. If the temperature or pressure is slightly changed, the density of defects will adjust and approach the new equilibrium very slowly, due to the need to break the strong covalent bonds for any change to occur. After the kinetic energy of all the atoms has relaxed to the new temperature, changes in the position of a defect along the chain require $\propto\exp(\beta\varepsilon)$ attempts, because of the large barrier for breaking these bonds. Such bond-breaking events enable just a single step along the chain, but in order to annihilate two defects must approach each other by diffusing significant distances along the chain.

Conversely, it is not enough to break a single molecular bond in order to create a pair of defects, because there is an overwhelming probability for the two single atoms just created to recombine, leaving the system in the original state. Two adjacent molecules in a given domain must simultaneously have broken bonds in order for a single novel molecule to form, nucleating a new domain, and creating two boundary defects. The probability of this happening is $\propto\exp(-2\beta\varepsilon)$ per attempt, and it must further be taken into account that once the two newborn defects begin diffusing, the probability that they will recombine is relatively large, because of their proximity to each other. Only after a long relaxation process will the concentration of defects reach its new equilibrium value, and the creation and annihilation processes reach equal rates (detailed balance).

III.3 Bond-counting potential with $\mathbf{z=2}$ and double bonds

We now extend the treatment of the previous subsection to a model with two bonds per atom, where each atom can have up to two single bonds (one on each side), or one double bond with one of its nearest neighbors. We assume that $d>a/2$, so that next-nearest neighbors cannot penetrate the bonding zone and two single bonds on the same side are forbidden. Comment2 Each single bond has an energy $-\varepsilon$ with bond length $\tilde{a}<x<a$, and a double bond contributes an energy $-\tilde{\varepsilon}$ with bond length $d<x<\tilde{a}$. The appropriate transfer matrix is

[TABLE]

where $A$, $B$, and $C$ are defined in analogy with Eqs. (17) and (18). Powers of this matrix contain terms where double bonds, $A$, are always flanked by non-bonds, $C$, whereas single bonds, $B$, are not limited in this manner (but cannot, of course, be adjacent to a double bond, $A$). The phase in which each pair of atoms has a double bond (i.e., the state whose contribution to the partition function is $\ldots ACACAC\ldots$) corresponds to the molecular phase, whereas the phase in which each atom is bonded to its two neighbors by single bonds (i.e., the state whose contribution to the partition function is $\ldots BBBBBB\ldots$) corresponds to a polymeric liquid.

Consider a possible transition between these two phases at $T=0$. The enthalpies associated with each of these bonds are $H_{A}(P)=-\tilde{\varepsilon}+Pd$, $H_{B}(P)=-\varepsilon+P\tilde{a}$ and $H_{C}(P)=Pa$ (this is the minimum enthalpy obtained for the minimum bond length in each case). The corresponding enthalpies per atom are

[TABLE]

A transition will occur for the case $\tilde{\varepsilon}>2\varepsilon$ and $\left(d+a\right)>2\tilde{a}$, for which the molecular phase is energetically favorable but the polymeric phase becomes favorable as the pressure is increased. The enthalpies of the two phases are equal at the critical pressure

[TABLE]

This model has sufficient detail to have a finite free-energy cost for a domain boundary. In order to calculate this cost, consider replacing a double bond by a single one, namely, replacing one of the $A$’s in the term $\ldots ACACAC\ldots$ by $B$. Introducing a larger domain of $B$’s at $T=0$ and $P=P_{\mathrm{C}}$ will cost the same free energy,

[TABLE]

Thus, the two phases are not mixed — for any length per atom in the range $\tilde{a}<\ell<(d+a)/2$, the system will have a clear tendency to phase-separate at zero temperature.

As per Landau’s argument for the absence of phase transitions in 1D systems, entropic effects will always be dominant in the thermodynamic limit, and the resulting finite concentration of phase-boundary defects will round off the phase transformation. However, as the energy scales for covalent bonds can easily be more than an order of magnitude larger than room temperature, it makes sense to also consider the behavior of finite systems; for a system of $N$ atoms, the whole system will be predominantly in one phase for $k_{\mathrm{B}}T\ll\Delta H/\ln N$ (in principle, at finite temperatures $\Delta H$ should be replaced by $\Delta G$).

This transition between a molecular phase and a polymeric phase, obtained for a 1D model system, is closely analogous to the known LLPT in phosphorus discussed in the introduction. KY00 ; KY02 ; MG03 Due to the phase-separation just discussed, the signature of a first-order phase transition is clearer in this case, when compared to the signatures of phase transitions in Takahashi solution to simpler 1D models with nearest-neighbor pair potentials. SLMR99 The reason for this is that in such solutions, the type of bonding between two neighboring atoms has no effect on their other neighbors, and thus no mechanism for forming a finite domain-boundary energy exists.

A variant of the present model can also reproduce a non-equilibrium phase transition, closely analogous to boiling of white phosphorus. In order to obtain a boiling transition, one must include weak attractive forces between molecules (in $v_{\mathrm{non-bond}}$). The polymerization transition can be dynamically hindered at low temperatures, by preparing the system as a low density molecular gas, and then bringing about condensation, i.e., increasing the density in a non-equilibrium manner. If the binding energy of each doubly-bonded molecule is much larger than the temperature (say, by a factor of 50), one can reach the condensed molecular phase, even if the polymeric phase is the thermodynamically stable one, because nucleation events of the polymeric phase will be exceedingly rare. Next, one may fix pressure and change the temperature. In 3D, one may have different solid, liquid and gaseous molecular phases, but for molecules in 1D, only a fluid phase exists. The transition between a liquid-like and a gas-like system is gradual in 1D, rather than abrupt, but one may still define a characteristic temperature for this transition. If this temperature is no longer very small compared to the molecular binding energy, nucleation of the polymeric phase is expected to occur at a rate which is no longer negligible. This is qualitatively what happens in actual experiments on boiling of white phosphorus. DTW46

IV Summary and discussion

We have suggested a simple type of classical many-body model potential, which imposes constraints on (i.e., counts) the number of bonds allowed for each atom, thus describing the essence of chemical valency. An analytical solution of the equilibrium statistical mechanics of a one-bond-per-atom model ($z=1$) in 1D using the transfer matrix method, including the relative density of diatomic molecules and the structure factor, was presented. Our solution generalizes the Takahashi solution for nearest-neighbor pair potentials. The results for the specific case of a constant attractive potential $-\varepsilon$ ($\varepsilon>0$) inside the bonding zone ($d<x<a$, where $d$ and $a$ are the hard core and bonding zone diameters) and zero outside of it ($x>a$) were displayed in detail.

The case of $z=2$ with double bonds was also discussed. This model exhibits a first-order phase transition (at $T=0$ in the thermodynamic limit) between a molecular phase and a polymeric phase, reminiscent of the LLPT observed in phosphorus. Furthermore, this model can also describe a non-equilibrium phase transformation analogous to the boiling of white phosphorus.

There are many open directions for future work with bond-counting potentials. One may introduce greater variety, by considering different valencies $z$ and also mixtures of several different species of atoms. One may study more realistic separation-dependent potentials rather than constants, or more complicated $q_{i},q_{j}$-dependence of $\mathcal{V}_{ij}$ (see Sec. II). In particular, the potential felt by a single atom approaching another single atom could be attractive and continuous at $y=a$, while at the same time the potential for an atom approaching a molecule could be repulsive already for $y>a$. On the other hand, one may study even simpler systems by considering, e.g., the case $d=\epsilon=0$ (or $\epsilon\to\infty$), resulting in a system which has only the integer $z$ parameter in addition to the length parameter $a$. As for the hard-sphere model, the temperature-dependece in such a system is trivial, and the properties, including phase transitions, depend only on $z$ and the reduced density $na^{3}$.

An extension of the model potentials proposed in this paper to 2D and 3D using molecular dynamics simulations may help to gain further insight into the mechanisms responsible for LLPTs and liquid state anomalies. A particularly interesting question is the mechanism behind the LLPTs in bismuth. In fact, one of the present authors has speculated that the covalently-bonded bilayer sheets of crystalline bismuth become only partly disordered at melting ($545\,$K), and fully disintegrate only at the LLPT at $1010\,$K. AN10 This would mean that the phase obtained by melting of solid bismuth consists of 2D polymers (the bilayers), which are unable to pass through each other, but able to glide in the transverse directions, possibly becoming crumpled, thus adding to the disorder.

This suggested mechanism also provides a reasonable context for explanation of the recently-observed anomalies of Refs. SY17, ; SY16a, ; SY16b, , concerning bismuth samples which underwent a high-pressure and high-temperature treatment. The 2D polymers which could result form the cooling and solidification procedure at $2\,$GPa, a pressure at which there is hardly any difference in density between the liquid and the solid, would perhaps have nontrivial topologies — possibly the sheets are rolled up into tubes which are entangled on larger scales. This assumes that the covalent bonds are strong enough to remain intact during the phase transitions from the melt to the “bismuth II” and then to the “bismuth I” solid phases (the quotes are a reminder that these samples contain significant structural defects; SY16a the character of these defects is consistent with this type of explanation). The stability of the complex topology of these sheets could then also be responsible for the absence of flow in the “melted” samples. SY16b An encouraging hint comes from additional studies of the phase behavior of bismuth at $\sim 2\,$GPa and elevated temperatures, which may be interpreted as indicating that at these pressures and at temperatures near melting the bonds retain a certain degree of stability even while the sample undergoes a phase transition. These studies identified cases in which the type of phase observed in this range of $P$ and $T$ depended on the path previously followed, i.e., on which phase the sample was previously in. LC17 ; PE06

It is natural to study the qualitative aspects of such suggested mechanisms using a $z=3$ bond-counting potential in 3D. The complex topology of crumpled 2D polymers suggests that topological defects will play an important role in molecular dynamics or Monte Carlo simulations of such systems. In this sense, the present 1D work may be viewed as a preliminary step in developing this approach, involving topological defects which are particularly simple.

Acknowledgements.

We acknowledge fruitful discussions with E. Pazy and V. V. Brazhkin.