text:
# Ab initio quasiharmonic equations of state for dynamically-stabilized soft-mode materials

## I Introduction

The quasiharmonic approximation  provides a means of extracting finite-temperature properties of materials from static calculations. It assumes the vibrational properties can be understood in terms of excitations of non-interacting harmonic normal modes: phonons. Lattice dynamics  can be used to calculate phonon energies by evaluating the eigenvalues of the dynamical matrix, which involves second derivatives of the crystal energy with respect to atomic displacements.

The frequencies of these modes depend on the crystal’s density. Hence they have a temperature-dependence that arises simply because of thermal expansion in the material.

Recent developments in ab initio energy calculations have enabled full phonon dispersion curves to be obtained, leading to a resurgence of interest in the quasiharmonic approach. Difficulties arise, however, when the dynamical matrix has negative eigenvalues, indicating that the crystallographic structure is not a local minimum of energy. Such crystals may be dynamically stabilized: because of their large entropy, they may represent a minimum of free energy at high temperature. Within the conventional assumption of harmonic phonons, the quasiharmonic framework leads to divergent free energies in these cases. Such systems have been treated numerically from first principles using Monte Carlo methods with effective Hamiltonians  or molecular dynamics simulations of a reduced set of modes . Here, we relax the quasiharmonic assumption of harmonic modes while retaining the approximation of non-interacting phonons and show how the intrinsic anharmonicity of such modes can be included in analytic free energy calculations.

The mineral periclase (MgO) is of some geological significance as one of the supposed constituents of the Earth’s lower mantle. It is generally believed that along the geotherm—the conditions of pressure and temperature actually occurring in the mantle—periclase remains in a single phase. Under other conditions, however, previous calculations have suggested that periclase has two phases in its solid state: a sodium chloride-like face-centered cubic phase (B1) and a cesium chloride-like simple cubic phase (B2) that is favored at extremely high pressures .

Imaginary phonon frequencies are found in the B2 phase of periclase and have attracted enormous attention in the other principal constituent of the lower mantle, magnesium silicate perovskite (MgSiO3) .

Equilibrium structures, thermodynamic properties and compositions depend on the free energy. Here we present a calculation of the free energy of periclase as a function of density and temperature. We use the pseudopotential-plane wave approach to evaluate total energies and the method of finite displacements to evaluate pressure-dependent force constants , including effective charges and dielectric constants  for the longitudinal optic modes. Based on calculated phonon frequencies and the quasiharmonic approximation, we present a first-principles calculation of the phase diagram and thermodynamic equation of state of solid periclase: the relationship between pressure, density and temperature.

## II The quasiharmonic method

### A Ab initio calculation of specific Helmholtz free energies

The first stage of the calculation is to obtain the specific (with respect to mass) Helmholtz free energy of each phase as a function of density and temperature. We write the free energy as the sum of the frozen-ion interaction energy and the free energy due to lattice vibrations.

#### 1 The frozen-ion energy

The frozen-ion energy—the interaction energy of the crystal with the ions fixed in their equilibrium positions—is, by definition, temperature-independent. Hence the free energy is simply equal to the internal energy. In order to determine the dependence on density, total energy densiyty functional calculations are carried out for each phase at a range of different lattice parameters.

#### 2 The lattice thermal Helmholtz free energy

We calculate the lattice thermal Helmholtz free energy of each phase as a function of density and temperature within the framework of the harmonic approximation .

For a range of lattice parameters, we evaluate the matrix of force constants for a supercell (several unit cells) of the phase under consideration, subject to periodic boundary conditions. We evaluate the forces on the ions in a crystal when one ion is displaced slightly from its equilibrium position: from such calculations the matrix of force constants may be constructed.

We denote the matrix of force constants by $`\varphi `$, where $`\varphi _{l,n,\alpha ;m,p,\beta }`$ is the component of force in direction $`\alpha `$ on ion $`n`$ in unit cell $`l`$ when ion $`p`$ in unit cell $`m`$ is displaced infinitesimally in direction $`\beta `$, divided by the magnitude of the displacement.

Pairs of density-functional calculations are carried out with an ion displaced from equilibrium along one of the Cartesian axes by a small amount in first a positive and then a negative sense. By averaging the resulting Hellmann-Feynman forces on the ions from the first simulation with the negative of the forces from the second, first-order anharmonic contributions to the force constants are eliminated.

The set of rotations under which the crystal structure is invariant are identified and the rotation matrices, together with the mappings between the ions under the symmetry operations, are evaluated. For a given pair of ions $`(l,n)`$ and $`(m,p)`$, the matrix of force constants $`\varphi _{l,n,\alpha ;m,p,\beta }`$ transforms as a second-rank tensor. However, for a symmetry operation, the transformed matrix must be the same as the (unrotated) matrix of force constants between the pair of ions which are mapped to $`(l,n)`$ and $`(m,p)`$. Hence, new elements of the matrix of force constants can be obtained by application of these point symmetries. Translational symmetries can be identified and exploited in a similar fashion.

The matrix of force constants should be symmetric  and Newton’s third law must be satisfied: if an ion is displaced slightly then the restoring force on that ion must be equal and opposite to the total force on all of the other ions. Hence we must have:

$`\varphi _{l,n,\alpha ;m,p,\beta }`$ $`=`$ $`\varphi _{m,p,\beta ;l,n,\alpha }`$ (1)
$`\varphi _{l,n,\alpha ;l,n,\beta }`$ $`=`$ $`{\displaystyle \underset{(m,p)(l,n)}{}}\varphi _{l,n,\alpha ;m,p,\beta }.`$ (2)
The force constants are obtained from separate numerical calculations; hence small violations of these requirements may occur. These two conditions are therefore alternately imposed on the matrix of force constants until further application leaves the matrix unchanged .

The next stage of the calculation involves the construction of dynamical matrices  for various wavevectors in the first Brillouin zone. Diagonalization of the dynamical matrix for a given wavevector gives the spectrum of corresponding eigenfrequencies. Strictly, these are only exact when the wavelengths are commensurate with the dimensions of the supercell . However, provided that the resulting dispersion curves are smooth, it may be assumed that the interpolation errors are negligible.

Cochran and Cowley  have shown that the elements of the dynamical matrix for an ionic crystal can be written as the sum of a term that behaves analytically as the wavevector tends to zero and a term that is non-analytic at the zone center. The latter term vanishes as the boundary of the Brillouin zone is approached. This term arises because longitudinal optic (LO) phonons cause an electric polarization field to be set up within the crystal as the oppositely-charged ions are displaced in opposite directions. The resulting long-range interactions cannot be calculated within the framework we have described so far because of the limited size of the simulation supercell. At the zone center itself, the LO phonon sets up a uniform electric polarization that is incompatible with the periodic boundary conditions on the supercell.

Cochran and Cowley’s expression for the dynamical matrix is:

$`\stackrel{~}{\varphi }_{n,\alpha ;p,\beta }(𝐤)=`$ $`\stackrel{~}{\varphi }_{n,\alpha ;p,\beta }^N(𝐤)+{\displaystyle \frac{4\pi e^2}{\mathrm{\Omega }|𝐤|^2\sqrt{M_nM_p}}}`$ (3)
$`\times `$ $`\left({\displaystyle \underset{\gamma =1}{\overset{3}{}}}k_\gamma Z_{n,\gamma ,\alpha }(𝐤)\right)^{}ϵ_0^1(𝐤)`$ (4)
$`\times `$ $`\left({\displaystyle \underset{\gamma =1}{\overset{3}{}}}k_\gamma Z_{p,\gamma ,\beta }(𝐤)\right),`$ (5)
where $`𝐤`$ is the wavevector, $`e`$ is the electronic charge, $`\mathrm{\Omega }`$ is the volume of the unit cell, $`M_n`$ is the mass of ion $`n`$, $`Z_{n,\alpha ,\beta }(𝐤)`$ is the Born effective charge tensor for ion $`n`$ and $`ϵ_0(𝐤)`$ is the electronic (frequency dependent) dielectric function. The first term on the right-hand side is the component of the dynamical matrix that is analytic as $`𝐤0`$, while the second term is the non-analytic part due to macroscopic polarization effects. We use our matrix of force constants evaluated using the Hellmann-Feynman theorem in a cubic supercell to evaluate the analytic part as:

$$\stackrel{~}{\varphi }_{n,\alpha ;p,\beta }^N(𝐤)=\frac{1}{\sqrt{M_nM_p}}\underset{m}{}\varphi _{0,n,\alpha ;m,p,\beta }e^{i𝐤(𝐑_0𝐑_m)},$$
(6)
where $`𝐑_m`$ is the position vector of unit cell $`m`$.

Following Parlinski et al we assume that the second term on the right-hand side of Equation 5 falls off from its value at the Brillouin zone center with a Gaussian profile. For wavevectors in the first Brillouin zone this term is:

$`{\displaystyle \frac{4\pi e^2}{\mathrm{\Omega }|𝐤|^2\sqrt{M_nM_p}ϵ_0(\mathrm{𝟎})}}\left({\displaystyle \underset{\gamma =1}{\overset{3}{}}}k_\gamma Z_{n,\gamma ,\alpha }(\mathrm{𝟎})\right)^{}`$ (7)
$`\times `$ $`\left({\displaystyle \underset{\gamma =1}{\overset{3}{}}}k_\gamma Z_{p,\gamma ,\beta }(\mathrm{𝟎})\right)\times e^{\left(\frac{|𝐤|}{\rho _0\kappa ^{1/2}}\right)^2},`$ (8)
where $`\kappa ^{1/2}`$ is the distance from the center to the boundary of the Brillouin zone along the $`k_x`$-, $`k_y`$\- and $`k_z`$-directions and $`\rho _0`$ is a parameter determining the rate at which the term falls off as the edge of the Brillouin zone is approached. Following Parlinski, we set $`\rho _01.2`$.

The frequency density-of-states function is evaluated using the method of Swift  in which the Brillouin zone is sampled using Monte Carlo methods. For a single harmonic mode of frequency $`\omega `$, the Helmholtz free energy is given by:

$$F_1(\omega )=k_BT\mathrm{log}\left(e^{\beta \mathrm{}\omega /2}e^{\beta \mathrm{}\omega /2}\right),$$
(9)
where $`\mathrm{}`$ is the Dirac constant, $`k_B`$ is Boltzmann’s constant, $`T`$ is the temperature and $`\beta =1/k_BT`$. Hence, by numerically integrating the product of the specific density-of-states with the mean free energy of a normal mode, the specific lattice thermal free energy can be calculated for a range of temperatures. The dependence on density is found by interpolating between the results at different lattice parameters.

### B Polymorphism

### C The Gibbs free energy

We have calculated the Helmholtz free energy $`f(v,T)`$ as a function of temperature $`T`$ and specific volume $`v`$ (the reciprocal of the density). However, the appropriate thermodynamic potential for constructing the $`(p,T)`$-phase diagram and evaluating the polymorphic equation of state is the specific Gibbs free energy $`g(p,T)`$, where $`p`$ is the pressure. The Gibbs free energy function for each phase can be evaluated using the Legendre transformation:

$$g(p,T)=f+pv=f\left(\frac{f}{v}\right)_Tv.$$
(10)
### D The phase diagram

Under conditions of fixed pressure and temperature, the system consists entirely of the available phase with the lowest Gibbs free energy. Thus the phase diagram in $`(p,T)`$-space can be evaluated.

### E Combining phases

For each pressure and temperature, we may evaluate the polymorphic Gibbs free energy $`g_{\mathrm{poly}}(p,T)`$ as the lowest of the Gibbs free energies for each phase. Given this, we may carry out a Legendre transformation to the polymorphic Helmholtz free energy:

$$f_{\mathrm{poly}}(v,T)=g_{\mathrm{poly}}pv=g_{\mathrm{poly}}p\left(\frac{g_{\mathrm{poly}}}{p}\right)_T.$$
(11)
Differentiating this, we obtain the pressure as a function of specific volume and temperature: the desired polymorphic equation of state.

## III Extension of the quasiharmonic method to unstable phonons

### A Analytic model of soft-mode phonons

In minerals such as perovskites  it is possible to describe the transition from a high-temperature phase to a low-temperature phase of lesser symmetry as the “freezing in” of a finite amplitude of an unstable phonon of the high-symmetry phase, plus a finite strain on the unit cell. We consider the application of quasiharmonic ideas to these materials.

The simple harmonic model gives a negative energy and divergent free energy arising from the unstable modes. In reality, the soft-mode phonon is best described by a potential double-well with a local maximum at the mean structure, corresponding to the high-symmetry phase.

Let $`x_i`$ be a coordinate describing the structural feature involved in the phase transition at a particular wavevector. The corresponding normal mode can be modeled by considering the dynamics of the set of $`\{x_i\}`$ moving in fixed local potential double-wells . In the harmonic limit, normal modes are uncoupled. However, because we are considering finite displacements there will in general be coupling between our double-well oscillators, this being most pronounced around the phase transition and at high temperatures. Coupling can be approximately treated by renormalization .

Much work has been concentrated on the Landau model in which the double-well is a free energy in the form of a quartic polynomial $`V(x)=Ax^4B(T)x^2`$, where $`B(T)`$ changes sign with temperature through coupling to other modes.

Such a polynomial expansion of the total energy is also possible, perhaps incorporating still higher-order terms . However, analytic terms beyond second order imply phonon coupling. This is inconsistent with the harmonic approximation used to describe non-soft modes: even at high phonon number the normal modes are assumed to be harmonic and therefore independent of each other.

We propose instead to describe the entire soft-phonon branch via a double-well of form:

$$V(x)=\frac{1}{2}m\omega _0^2x^2+ϵ(e^{x^2/2\sigma ^2}1).$$
(12)
where $`ϵ`$, $`\omega _0`$ and $`\sigma `$ are wavevector dependent. Provided that $`ϵ>m\omega _0^2\sigma ^2`$, there are minima at:

$$x=x_\pm \pm \sqrt{2\sigma ^2\mathrm{log}(ϵ/m\omega _0^2\sigma ^2)},$$
(13)
separated by a barrier of height:

$`\mathrm{\Delta }V`$ $``$ $`V(0)V(x_\pm )`$ (14)
$`=`$ $`ϵm\omega _0^2\sigma ^2\left(1+\mathrm{log}\left(ϵ/m\omega _0^2\sigma ^2\right)\right).`$ (15)
This form of potential has the advantage of being approximately quadratic in both the low-energy and high-energy limits. Specifically, for the low-energy case at $`x=x_\pm `$, we have:

$$\frac{d^2V}{dx^2}=2m\omega _0^2\mathrm{log}\left(ϵ/m\omega _0^2\sigma ^2\right),$$
(16)
which is equivalent to an harmonic oscillator of frequency:

$$\omega _0^{}=\sqrt{2\omega _0^2\mathrm{log}(ϵ/m\omega _0^2\sigma ^2)}.$$
(17)
On the other hand, for the high-energy case, the potential approximates that of an harmonic oscillator of frequency $`\omega _0`$.

Soft modes do not usually show abnormal dependence on temperature—except in the vicinity of the phase transition. Therefore we expect our model to be more widely applicable than models where soft modes are treated as quartic and other modes as harmonic.

The (imaginary) harmonic frequency $`\omega _c`$ about the (unstable) center of the well is given by:

$$\omega _c^2=\omega _0^2\frac{ϵ}{m\sigma ^2}.$$
(18)
### B Isolated double-well oscillators

We now consider the problem of motion in an isolated potential double-well.

#### 1 Classical solution

To evaluate the mechanical energy, we assume the mode is in thermal contact with a heat bath at the appropriate temperature. The mean energy is given by:

$`E`$ $`=`$ $`{\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}H(p,x)e^{\beta H(p,x)}𝑑p𝑑x}{_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}e^{\beta H(p,x)}𝑑p𝑑x}}`$ (19)
$`=`$ $`{\displaystyle \frac{k_BT}{2}}+{\displaystyle \frac{_0^{\mathrm{}}V(\sigma z)e^{\beta V(\sigma z)}𝑑z}{_0^{\mathrm{}}e^{\beta V(\sigma z)}𝑑z}},`$ (20)
where $`H(p,x)=p^2/2m+V(x)`$ is the Hamiltonian of the isolated mode as a function of $`x`$ and $`p`$, the canonical momentum conjugate to $`x`$. $`\beta =1/k_BT`$ where $`k_B`$ is Boltzmann’s constant and $`T`$ is the temperature.

For a given energy $`E`$, the frequency of our isolated mode can be evaluated using the action-angle method. The action variable is:

$$j\frac{1}{2\pi }p𝑑x.$$
(21)
If the mode has energy $`Eϵ`$, so that there is sufficient energy to cross the barrier each libration, we find that:

$$j=\frac{2\sigma }{\pi }_0^{x_M(E)/\sigma }\sqrt{2mEm^2\omega _0^2\sigma ^2z^22mϵe^{z^2/2}}𝑑z,$$
(22)
where $`x_M`$ is the positive solution of $`V(x)=E`$. On the other hand, if $`E<ϵ`$, so that the motion is confined to one side of the double-well, we find that:

$$j=\frac{\sigma }{\pi }_{x_m(E)/\sigma }^{x_M(E)/\sigma }\sqrt{2mEm^2\omega _0^2\sigma ^2z^22mϵe^{z^2/2}}𝑑z,$$
(23)
where $`x_M`$ is the greater of the two positive solutions and $`x_m`$ is the lesser.

In either case, the corresponding frequency $`\omega `$ may be evaluated using Equation 24:

$$\frac{1}{\omega }=\frac{j}{E}.$$
(24)
Taken together, the results of Equations 20, 22, 23 and 24 allow us to calculate numerically the frequency of an isolated oscillator moving with the mean thermal energy as a function of temperature. Example results are shown in Figure 1. Typical “soft-mode” behavior is observed, with the frequency dropping to zero in a cusp at the transition temperature. This simple approach was used in early studies of MgSiO3 .

#### 2 Quantum solution

The energy eigenfunctions of a particle moving in a symmetric potential must be either symmetric or antisymmetric. Furthermore, by considering building up the Gaussian barrier adiabatically, it is clear that the symmetry of the eigenfunctions must be the same as for those of the harmonic oscillator.

The definite symmetry of the wavefunctions leads to a “paradox” for wells of finite separation. If we know the energy of the oscillator then it is in an energy eigenstate and the wavefunction is either symmetric or antisymmetric. Hence the probability distribution is symmetric about the center of the double-well and we cannot meaningfully say which side the oscillator is confined in, even if its energy is much less than the barrier height. Thus it is not conceptually clear that equating the mean thermal energy with the barrier height gives the correct transition temperature. We discuss this further in Section III D.

The Hamiltonian operator for a particle moving in the quadratic potential without the additional Gaussian potential is:

$$\widehat{H}^0=\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+\frac{1}{2}m\omega _0^2x^2.$$
(25)
The well-known energy eigenvalues and eigenfunctions for the time-independent Schrödinger equation $`\widehat{H}^0\varphi _n=E_n^0\varphi _n`$ are:

$$E_n^0=(n+1/2)\mathrm{}\omega _0$$
(26)
and

$$\varphi _n=\frac{2^{n/2}}{\sqrt{n!}}\left(\frac{m\omega _0}{\mathrm{}\pi }\right)^{1/4}e^{m\omega _0x^2/2\mathrm{}}H_n\left(\sqrt{\frac{m\omega _0}{\mathrm{}}}x\right)$$
(27)
for $`n𝒩_0`$, where $`H_n(x)`$ is the $`n`$th Hermite polynomial.

The Hamiltonian operator for a particle moving in the double-well is $`\widehat{H}=\widehat{H}^0+V_1`$ where $`V_1=ϵ(e^{x^2/2\sigma ^2}1)`$ is the extra Gaussian term. Let the eigenfunctions and eigenenergies of the full Hamiltonian be $`\psi _n`$ and $`E_n`$.

The eigenfunctions of the simple harmonic oscillator (SHO) are chosen as the basis of wavefunction space. This choice makes the computation particularly simple, as will be seen below.

The matrix elements of the Hamiltonian with respect to our chosen basis are:

$`\varphi _i|\widehat{H}\varphi _j`$ $`=`$ $`\varphi _i|\widehat{H}^0\varphi _j+\varphi _i|V_1\varphi _j`$ (28)
$`=`$ $`(i+1/2)\mathrm{}\omega _0\delta _{i,j}`$ (30)
$`+{\displaystyle \frac{ϵ}{2^{(i+j)/2}\sqrt{i!j!\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{Kz^2}H_i(z)H_j(z)𝑑z,`$
where $`K\mathrm{}/2m\omega _0\sigma ^2+1`$. Note that the matrix is real and symmetric.

The Hermite polynomials satisfy $`H_i(x)=(1)^iH_i(x)`$, so that $`\varphi _i|V_1\varphi _j=0`$ if $`i+j`$ is odd. If $`i+j`$ is even then the integrand is an even function. Hence, in this case:

$$\varphi _i|V_1\varphi _j=\frac{2ϵ}{2^{(i+j)/2}\sqrt{i!j!\pi }}_0^{\mathrm{}}e^{Kz^2}H_i(z)H_j(z)𝑑z.$$
(32)
The eigenvalues of the matrix of the Hamiltonian are the allowed energy levels. For energies that are large compared with the barrier height the particle will spend most of its time away from the center of the well; hence we expect that the system will behave as a SHO in this limit. Indeed, it is clear from Equation 32 that the elements of the matrix $`\varphi _i|V_1\varphi _j`$ fall off rapidly as $`i`$ and $`j`$ increase. Hence, for large $`i`$ or $`j`$, the eigenvalues of the Hamiltonian tend to to those of the SHO. Thus we only need to diagonalize the upper left-hand corner (say, the $`(n_c+1)\times (n_c+1)`$ submatrix) in order to obtain the first $`n_c+1`$ energy levels $`E_0`$ to $`E_{n_c}`$. Beyond $`n_c`$ the eigenvalues may be taken to be those of the SHO. The comparative ease with which the Hamiltonian matrix can be diagonalized is one of the advantages of the quadratic-plus-gaussian double-well over the quartic double-well potential, although there is no analytic form equivalent to equation 9

Assuming that $`n_c`$ is sufficiently large, the canonical partition function can be written as:

$`Z`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{n_c}{}}}e^{\beta E_n}+{\displaystyle \underset{n=n_c+1}{\overset{\mathrm{}}{}}}e^{\beta (n+1/2)\mathrm{}\omega _0}`$ (33)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{n_c}{}}}e^{\beta E_n}+{\displaystyle \frac{e^{\beta \mathrm{}\omega _0(n_c+1)}}{e^{\beta \mathrm{}\omega _0/2}e^{\beta \mathrm{}\omega _0/2}}}.`$ (34)
The free energy of the double-well oscillator can then be evaluated as:

$$F_1=k_BT\mathrm{log}\left(Z\right).$$
(35)
### C Practical implementation in the quasiharmonic method

Having proposed that each soft mode at a given wavevector be described by a double-well of the form given in Equation 12, we now describe how the parameters $`ϵ`$, $`\sigma `$ and $`\omega _0`$ can be determined. Note that if the $`\{x_i\}`$ are mass-reduced phonon coefficients then we may, without loss of generality, set $`m=1`$.

Consider the phonon dispersion curve of a crystal structure in which imaginary frequencies are present. Those branches that remain real throughout the whole of the Brillouin zone are treated as harmonic and, for each mode at each wavevector, Equation 9 may be used to find the corresponding free energy. For those branches that are imaginary in some region of the Brillouin zone, however, we propose the following treatment:

1. At each symmetry point of the Brillouin zone the eigenvector corresponding to the relevant mode should be evaluated and the displacement pattern frozen into the crystal. Ab initio techniques can then be used to find the corresponding low-symmetry structure, if desired. Using three total energy calculations with different amplitudes of the soft phonon frozen into the structure***In fact, if we have calculated the total energy of the high-symmetry phase (corresponding to zero amplitude of the phonon) then we only need to carry out a further two frozen phonon calculations., we may evaluate the parameters $`ϵ`$, $`\sigma `$ and $`\omega _0`$ of the double-well. Note that it is possible to fit equation 12 to every branch even if the mode is not imaginary since that Equation 12 does not necessarily describe a double-well. In practice the harmonic approximation is used for all-real branches, it is equivalent to setting $`ϵ=0`$.
2. For each branch we use our results for the double-well parameters at the symmetry points to construct interpolating polynomials over the whole of the Brillouin zone for the $`ϵ`$ and $`\sigma `$ parameters.
3. For any wavevector in the Brillouin zone, we may find the spectrum of corresponding (possibly imaginary) frequencies. Provided we know to which branch these modes belong, we have sufficient information to determine the parameters of the appropriate double-well for each mode. $`ϵ`$ and $`\sigma `$ are found by interpolating to our wavevector and the unstable frequency gives $`\omega _c`$ (see Equation 18), from which we may find $`\omega _0^2=\omega _c^2+ϵ/m\sigma ^2`$.
4. Hence, for any given wavevector, the free energy of each mode, whether harmonic or soft, can be evaluated. These free energies can be summed to give the free energy contribution from all modes at the given wavevector.
5. The free energy can then be integrated over all wavevectors in the Brillouin zone to give the total lattice thermal free energy. By using a grid-based scheme to integrate over an irreducible wedge of the zone and by making use of the continuity of the gradient of each branch, it is possible to keep track of which branch is which—necessary if the appropriate values of $`ϵ`$ and $`\sigma `$ are to be interpolated in the pressence of imaginary branch crossings. The problem of interpolation over an irreducible wedge of the Brillouin zone has been studied extensively in the context of electronic eigenvalues: see, for example, Reference .

The high-symmetry dynamically stabilized phase and the low-symmetry “frozen-phonon” phase can now be treated as being distinct. Hence the methodology of Section II B can be applied to find the phase diagram.

### D Interpretation of the soft mode transition

Consider an isolated symmetric double-well oscillator. For the probability density to be asymmetric—necessary if we are to meaningfully say that the particle is in one well or the other—we must have a mixture of symmetric and antisymmetric energy eigenstates. Therefore we cannot simultaneously know the energy of our particle and which well it is in unless we break the symmetry (e.g. by allowing the crystal to distort under phonon-strain coupling).

If the oscillating particle’s wavefunction is a superposition of different energy eigenfunctions then the expansion coefficients will evolve in time (provided the system remains both undisturbed and unobserved) according to the Schrödinger equation. Hence the quantum mechanical expectation value of the particle’s position, $`x`$ changes in time.

The time-dependent wavefunction can be written as:

$$\mathrm{\Psi }(x,t)=\underset{n=0}{\overset{\mathrm{}}{}}c_ne^{iE_nt/\mathrm{}}\psi _n(x).$$
(36)
Hence the expectation of $`x`$ can be written as:

$`x`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}c_n^{}c_me^{i(E_nE_m)t/\mathrm{}}\psi _n|x\psi _m`$ (37)
$`=`$ $`2{\displaystyle \underset{n\mathrm{even}}{}}{\displaystyle \underset{m\mathrm{odd}}{}}\left|c_n\right|\left|c_m\right|`$ (39)
$`\times \mathrm{cos}(\omega _{nm}t+\eta _{nm})\psi _n|x\psi _m`$
where $`\omega _{nm}=|E_nE_m|/\mathrm{}`$ and $`\eta _{nm}=(\mathrm{arg}(c_m)\mathrm{arg}(c_n))\mathrm{sgn}(\omega _{nm})`$. Note that we use the fact that $`\psi _n|x\psi _m=0`$ if $`\psi _n`$ and $`\psi _m`$ are either both odd or both even since $`x`$ itself is odd. We also use the fact that $`\psi _n|x\psi _m=\psi _m|x\psi _n`$.

We assume that the energy levels are initially populated according to Boltzmann statistics; thus $`|c_n|^2=Z^1e^{\beta E_n}`$, where $`\beta =1/k_BT`$ and $`Z=_{n=0}^{\mathrm{}}e^{\beta E_n}`$ is the canonical partition function.

So we have:

$$x=\underset{n\mathrm{even}}{}\underset{m\mathrm{odd}}{}\mathrm{\Gamma }_{nm}\mathrm{cos}(\omega _{nm}t+\eta _{nm}),$$
(40)
where $`\mathrm{\Gamma }_{nm}=2Z^1e^{\beta (E_n+E_m)/2}\psi _n|x\psi _m`$ is the amplitude of the sinusoidal component in the expansion of $`x`$ with frequency $`\omega _{nm}`$.

For the harmonic oscillator potential, the frequencies of the oscillations in $`x`$ are of the form $`\omega _{nm}=|E_nE_m|/\mathrm{}=|nm|\omega _0`$. Thus the lowest oscillation frequency is $`\omega _0`$. For the symmetric double-well, however, we end up with a set of pairs of energy levels that are very close to each other (becoming degenerate in the limit that the barrier height goes to infinity). These give rise to oscillation frequencies very much lower than $`\omega _0`$.

As the temperature is increased, higher frequency components have larger $`\mathrm{\Gamma }`$ coefficients. We suggest that the soft mode phase transition be judged to occur when the frequency with the highest coefficient exceeds the frequency of the experimental probe. When this has happened, the predominant sinusoidal component of $`x`$ has a frequency higher than can be measured by the experimental probe, and so it appears to the experimenter that $`x=0`$. Below this temperature, measurements of $`x`$ will tend to find it in one well or the other. This definition is different from the polymorphic one (Section II B) because of the contribution to the free energy from the symmetry-breaking distortion of the lattice that inevitably accompanies the transition. In particular, the quasiharmonic transition is first-order while this measurement dependent “soft mode” one is second-order .

Thus the first-order transition is determined by comparing:

1. The free energy of the soft-mode phase, calculated as above, expanded about an unstable frozen-ion structure without strain-phonon coupling.

with
2. The free energy of the low symmetry phase, calculated by expanding about the minimum of total energy.

Typically, the former will have higher entropy (sampling from both wells) while the latter has lower energy.

### E Absorption of low-frequency photons

Figure 2 shows the energy difference between neighboring energy levels $`(E_nE_{n1})/\mathrm{}`$ plotted against the mean of the two energies $`(E_n+E_{n1})/2`$ for the quantum double-well oscillator. Absorption of photons at frequency $`(E_nE_{n1})/\mathrm{}`$ is symmetry-allowed.

For energies in excess of the barrier height the frequency $`(E_nE_{n1})/\mathrm{}`$ is virtually identical to the classical frequency for energy $`(E_n+E_{n1})/2`$, obtained using the method of Section III B 1. In the very high-energy limit the frequency behaves as that of the quadratic potential well without the Gaussian barrier.

For energies less than the barrier height the energy levels tend to degenerate pairs of levels. The frequencies given by the difference between the energy levels of neighboring pairs again correspond to the classical frequencies. However, the pairs of almost-degenerate eigenstates imply the existence of very low-frequency absorption peaks. (These are the very low frequencies that alternate with the classical frequencies below the transition energy in Figure 2.) It should be noted that these frequencies are not associated with the normal modes of the low-symmetry phase and do not, therefore, contribute to the quasiharmonic thermal energy. They are a feature of the quantum double-well oscillator without classical analog.

## IV Application of quasiharmonic methods to periclase

### A Computational details

#### 1 The cold-curve

For each lattice parameter the total energy is evaluated using the CASTEP software package , which utilizes density-functional theory in the generalized gradient approximation (GGA) . The ionic cores are accounted for using ultrasoft pseudopotentials  in Kleinman-Bylander form . The wavefunctions of the valence electrons are expanded in a plane wave basis set up to an energy cutoff of 540 eV.

For the B1 phase, the simulation cell consists of a single cubic unit cell. The Brillouin zone is sampled at 20 special points generated from an $`8\times 8\times 8`$ mesh using the Monkhorst-Pack scheme . For the B2 phase, the simulation cell is a single cubic primitive cell. The Brillouin zone is sampled at 35 special points from a $`9\times 9\times 9`$ mesh. In each case the point symmetries of the crystal are enforced .

The equilibrium lattice parameter of the B1 phase at zero external pressure (which corresponds to the minimum of the cold-curve) calculated using CASTEP is $`a=4.259`$ Å, which may be compared with an experimentally determined parameter $`a=4.2115(1)`$ Å . The difference between the theoretical and experimental values is about 1%.

#### 2 Determination of the matrix of force constants

The supercells simulated to determine the matrices of force constants for the B1 phase consist of $`2\times 2\times 2`$ cubic unit cells (64 atoms). Thus the interactions between a given ion and its third-closest shell of neighbors are included in our calculations. For the B2 phase, the cells used consist of $`2\times 2\times 2`$ cubic primitive unit cells (16 atoms). In these supercells the crystal symmetry is such that only two ionic displacements are required to complete the entire matrix of force constants. The plane-wave cutoff energy is 540 eV and the Brillouin zone is sampled at 6 special points from a $`4\times 4\times 4`$ mesh. In each simulation the ion displaced from equilibrium is moved by $`0.4\%`$ of the lattice parameter.

As demonstrated by Parlinski  it is possible to calculate the Born effective charge tensors from first principles using simulations of elongated supercells. Note that because of the symmetry of the B1 and B2 phases, the Born effective charge tensors are isotropic (so $`Z_{n,\alpha ,\beta }=Z_n\delta _{\alpha ,\beta }`$). Furthermore, the sum of the Born effective charges over the ions in a unit cell must be zero  (so $`Z_{\mathrm{Mg}}=Z_\mathrm{O}`$). Hence there is effectively only one undetermined parameter in the non-analytic term: $`Z_{\mathrm{Mg}}(\mathrm{𝟎})/\sqrt{ϵ_0(\mathrm{𝟎})}`$.

We choose $`Z_{\mathrm{Mg}}(\mathrm{𝟎})/\sqrt{ϵ_0(\mathrm{𝟎})}=4.4`$ to give the LO branch in the dispersion curve of Figure 3. Our values for the effective charge tensors of the B1 phase are such that the calculated LO branch for lattice parameter 4.2 Å are in reasonable agreement the experimental results of Peckham . We neglect the variation of the effective charge tensors with lattice parameter.

### B Approximations and errors

#### 1 Errors in ab initio total energy calculations

Total energy differences between structures calculated using density-functional theory in the generalized gradient approximation are thought to be reliable to within a few percent .

The cutoff energy of the plane-wave basis at 540 eV is sufficient for convergence of the total energies of the crystals to within 10 $`\mu `$eV per ion, several orders of magnitude less than the likely error due to the use of the GGA. The Hellmann-Feynman forces are converged to within 1 meV Å-1, at least two orders of magnitude less than the dominant forces arising when an ion is displaced.

#### 2 The harmonic approximation

We investigate the range of validity of the harmonic approximation. This is done for the B1 phase with lattice parameter $`a=4.2`$ Å.

We evaluate the force constant of the restoring force on a magnesium ion as it is displaced in the $`x`$-direction. The results are shown in Figure 4. It can be seen that the force constant starts to increase when the ionic displacement reaches $`a_{\mathrm{max}}0.084`$ Å, about 2% of the lattice parameter, at which point the potential energy is about 0.4 eV. Other displacements are similar; hence it is reasonable to assume that the forces remain linear (and the quasiharmonic assumption is valid) for temperatures up to several thousand Kelvin.

#### 3 Other sources of error

Another potential source of error is the limited size of the simulation supercells. However, the results shown in Table I for the B1 phase make it clear that interactions beyond the third-closest shell of neighbors may be safely neglected .

Contributions to the free energy from the thermal excitation of the valence electrons, from coupled electron-phonon excitations and from the equilibrium population of defects are thought to be negligible in comparison with the frozen-ion and lattice thermal energies .

Poor convergence of the Hellmann-Feynman forces can result in the violation of Newton’s third law for the matrix of force constants. Typically this results in the acoustic branches of the dispersion curve failing to pass through zero at the center of the Brillouin zone . As discussed in Section II A 2, Newton’s third law is imposed on the matrix of force constants. However, even without this, the calculated acoustic branches pass very close to zero at the zone center.

The method by which long-range polarization effects are accounted for is also approximate. The effects of this are discussed below.

### C Ab initio phonons

#### 1 The B1 phase

Inelastic neutron scattering experiments were carried out by Peckham  and used to generate dispersion curves for the B1 phase . We compare our theoretical dispersion curve with these results in Figure 3. (Note that our dispersion curve was generated for lattice parameter 4.2 Å, whereas Peckham’s results were obtained under ambient conditions where the lattice parameter is 4.212 Å.) Our theoretical results are in reasonable agreement with experiment. (The lattice parameter of $`4.2`$ Å corresponds to a pressure of about 7 GPa at zero temperature.)

The specific frequency density-of-states function is shown in Figure 5. Without the addition of the non-analytic term to the dynamical matrix, the longitudinal optic branch is degenerate with the transverse optic branch at the $`\mathrm{\Gamma }`$-point, and this is also shown. Although only the LO branch is altered substantially, it can be seen that the inclusion of the non-analytic term has a significant effect on the density of states.

We compare sound velocities calculated from our dispersion curves with the experimental results of Reichmann et al  obtained using ultrasonic interferometry. Reichmann obtains a P-wave sound speed of 9119 ms-1 in the $`[100]`$-direction whereas our longitudinal-acoustic mode has $`\left(d\omega _{\mathrm{LA}}/dk_{[100]}\right)_{𝐤=0}=11367`$ ms-1.

On the other hand, in the $`[111]`$-direction, the experimental P-wave velocity is 10125 ms-1 which may be compared with our theoretical value of 10818 ms-1.

Thus Reichmann’s experimental results show a higher degree of anisotropy than do our theoretical results. The discrepancy would appear to be (at least partly) caused by the imposition of symmetry and Newton’s third law on the matrix of force constants : if this procedure is not carried out then we find that $`\left(d\omega _{\mathrm{LA}}/dk_{[100]}\right)_{𝐤=0}=10770`$ ms-1 and $`\left(d\omega _{\mathrm{LA}}/dk_{[111]}\right)_{𝐤=0}=12651`$ ms-1. Although Reichmann’s P-wave velocities are still somewhat less than these theoretical velocities, both are now in agreement that the velocity in the $`[100]`$-direction is less than the velocity in the $`[111]`$-direction.

#### 2 The B2 phase

Typical dispersion curves for the B2 phase at lattice parameters 2.0 Å and 2.7 Å are shown in Figures 6 and 7. At zero temperature these lattice parameters correspond to pressures of $`653`$ GPa and $`6`$ GPa respectively.

Note the presence of unstable modes at low pressures in the dispersion curve of the B2 phase (Figure 7). We find that the B2 phase is structurally unstable for pressures below about 82 GPa. In our calculations for the phase boundary in periclase we do not require our novel method for dealing with soft modes because the imaginary frequencies are only found in the B2 phase at pressures for which the B1 phase is clearly favored.

### D Ab initio equation of state

We plot the pressure against specific volume for a range of temperatures in Figure 8. This is the desired thermodynamic equation of state for periclase. Also shown is a third-order Birch-Murnaghan equation of state generated from the isothermal bulk modulus and its first derivative with respect to volume, which were obtained by means of ultrasonic sound velocity measurement .

As a test of the validity of our results, the isothermal bulk modulus at zero temperature and external pressure (where periclase is entirely in the B1 phase) can be compared with experimental results. We calculate the bulk modulus to be $`v\left(p/v\right)_{T,v}=155`$ GPa, whereas an experimentally determined value is $`160.3\pm 0.3`$ GPa . The theoretical and experimental results differ by about 3%.

We may also compare the pressure derivative of the bulk modulus at zero pressure and temperature with experimental results. The bulk modulus is found to be almost, but not quite, a linear function of pressure. Fitting a straight line to the data from $`8.3`$ GPa to $`21.12`$ GPa we find that the the gradient is $`4.11`$, in excellent agreement with the experimentally determined value of $`4.2\pm 0.2`$ . It is found that the pressure derivative of the bulk modulus decreases slightly as the pressure is increased.

### E Ab initio phase diagram

The theoretical phase diagram of solid periclase is shown in Figure 9. At pressures and temperatures below and to the left of the phase boundary shown in the diagram, periclase exists in the B1 phase; above and to the right of the boundary it exists in the B2 phase. Duffy  has shown experimentally that, at room temperature, the B1 phase is stable to pressures of at least $`227`$ GPa. This lies well within the B1 region of our theoretical phase diagram. It can also be noted that the B2 phase is not favored at the pressures at which it is structurally unstable.

Also shown in Figure 9 is the theoretical phase diagram obtained by Strachan et al  using molecular dynamics simulation. Although qualitatively similar, the calculated position and orientation of the B1–B2 phase boundary is very different. Karki  obtained a B1–B2 transition pressure of 460 GPa at 0 K, which also differs substantially from that of Strachan. Figure 10 shows the Gibbs free energy plotted against pressure for the two phases at two different temperatures. The difficulty in ascertaining the transition pressures at which the curves cross is apparent. This consideration will affect all theoretical calculations of the phase diagram of periclase.

We confirm the difficulty in locating the transition by attempting to reproduce Karki’s zero-temperature results in which zero-point lattice vibrational energy is neglected. We simply calculate the enthalpy against pressure for the two phases using CASTEP. The results are shown in Figure 11. We find a transition pressure of GPa, different from that of Karki (460 GPa). The possibility of a metallization transition lowering the energy of the B2 phase was investigated but found not to occur at relevant pressures.

We consider the fractional error in unit cell volume for the B2 phase to which the difference between the specific enthalpies of the B1 and B2 phases at 451 GPa corresponds. For given unit cell enthalpy for the B2 phase, we find the fractional volume difference as:

$$\left|\frac{\mathrm{\Delta }V}{V}\right|=2\left|\frac{h_{B2}h_{B1}}{h_{B2}+h_{B1}}\right|,$$
(41)
where $`h_{B1}=34.4043`$ eVÅ-3 and $`h_{B2}=34.3978`$ eVÅ-3 are the specific enthalpies of the two phases at 451 GPa. Thus a discrepancy of $`\mathrm{\Delta }V/V\times 10^40.02`$% in the volume of either phase would correspond to a 200GPa change in the transition pressure. This illustrates the sensitivity of the transition point to the details of how the total energy calculations are carried out.

We use the generalized gradient approximation and ultrasoft pseudopotentials, whereas Karki used the local density approximation and $`Q_c`$ tuned pseudopotentials . The difference between calculated volume using these two methods is typically of the order of 5%; hence this is likely to be responsible for the difference between our results and Karki’s.

The Clausius-Clapeyron equation for the coexistence line between two phases in $`(p,T)`$-space is:

$$\frac{dp}{dT}=\frac{\mathrm{\Delta }s}{\mathrm{\Delta }v},$$
(42)
where $`p(T)`$ is the coexistence line and $`\mathrm{\Delta }s`$ and $`\mathrm{\Delta }v`$ are the specific entropy and volume differences between the two phases across the line. At zero temperature the entropy of the two phases should be zero by the third law of thermodynamics which is valid within our method (though not, e.g. in classical MD); therefore, provided the phases have different densities, the coexistence line should satisfy $`dp/dT=0`$ at $`T=0`$. Figure 9.) does not appear to satisfy this requirement.

The explanation for this apparent discrepancy lies with the fact that the densities of the two phases are very similar (and converging) at their predicted zero-temperature transition point: the specific volumes of the B1 and B2 phases are 0.202 Å3amu-1 and 0.198 Å3amu-1 respectively.

## V Conclusions

We have described an extension to the quasiharmonic method that allows the free energy contribution from “soft” phonons in dynamically stabilized crystals to be evaluated. Our approach is based on a form of potential double-well different to that used in previous work on soft phonons: a parabola-plus-Gaussian form that has the advantage of being harmonic in both the low- and high-temperature limits.

We argue that the first-order nature of the phase transitions found using our extended quasiharmonic method arises because of the coupling of the relevant phonon to strain in the crystal. Without this coupling, the transition would be second-order. We have suggested a criterion for judging when a second-order soft-mode phase transition has occurred, taking into account the quantum mechanical nature of the problem.

At energies less than the height of the central barrier in our symmetric potential double-well, the allowed energy levels consist of near-degenerate pairs. Hence we suggest there must exist extremely low-frequency photon absorption peaks for soft-mode materials in their low-temperature phase, corresponding to photon-induced transitions between such pairs.

We have evaluated the equation of state and the phase boundary for the B1–B2 transition in periclase using ab initio calculations in the quasiharmonic approximation. We predict that this transition will occur, but that it is well outside the ranges encountered inside the Earth. Locating the B1–B2 phase boundary with precision is difficult, however, because the Gibbs free energy curves of the two phases are very similar when they cross.

## VI Acknowledgments

We would like to thank D. C. Swift and M. C. Warren for helpful discussion. This work was carried out as part of the UKCP-MSI collaboration, supported by EPRSC GR/N02337. Copies of the codes used for the free energy calculations are available from the authors on request.


format: text

text:
# Untitled Document

Revision of upper estimate of percolation threshold on square lattice

Yu.P. Virchenko, Yu.A. Tolmacheva

The more exact upper estimate of the percolation threshold for the site problem on the quadratic lattice $`^2`$ have been found on the basis of the cluster decomposition. It is done by the number estimate of cycles on $`^2`$ which maybe external boundaries of finite clusters.

1. Introduction

Generally, percolation theory studies the connectedness relation for random sets in topological spaces. Problems representing greatest interest in the theory appear in noncompact topological spaces.

The existence of a noncompact connected component in random set realizations with nonzero probability is the crucial question in the theory . In such a general setting of the problem, however, it is impossible to obtain any strong results. That is why the basic objects of the percolation theory are random sets in $`^d`$ and $`^d`$ at $`d=2,3`$ generated by stationary random fields . The first case corresponds to so-called continuous percolation theory and the second one – to discrete percolation theory. This limitation is a consequence of the fact that it is the problem set in such a way is mainly needed in physical applications. Therefore, it may be considered as the object of mathematical physics. However, despite the above restriction, the basic problem is very complicated and hardly lends itself to exact mathematical study. Rigorous results referring to discrete percolation theory which is its mostly developed direction have been summarized in the monograph . Later results and also results referring to continuous theory have been summarized in the excellent review . Physical literature dealing with different nonrigorous heuristic approaches to the investigation of the percolation phenomenon, computer experiments and special applications of the theory in theoretical physics is giant and we shall not pay attention to it in this short introduction to the problem.

Despite of the great flow of publications at 80-th in physical literature and simultaneously appearance of some reviews in which mathematical results are summarized, main questions of the percolation theory considered as the object of mathematical physics remain without answers up to now. It is the situation even for simplest case being studied in the theory when random set is generated by Bernuolli’s field on $`^2`$. Even for this case any algorithm of the percolation threshold calculation (see §2) and, more generally, any algorithm of percolation probability evaluation was not found for an arbitrary periodic graph . Our communication is dealing with the revision of upper estimate of percolation threshold as compared to what may be found in literature  for the periodic graph on $`^2`$ which is called the square lattice.

2. Percolation theory problem on $`^2`$

Let us consider an infinite graph with the vertex set $`^2`$. For simplicity of statements and arguments, we shall study it as one immersed into $`^2`$. The adjacency relation on the graph is defined by the set of pairs $`\{(x,y)^2:x\phi y\}`$ where $`x\phi y`$ if and only if $`y=x\pm 𝐞_1`$ or $`y=x\pm 𝐞_2,𝐞_1=(1,0),𝐞_2=(0,1)`$. We shall follow the terminology of statistical physics and that is why we shall call such a graph the square lattice and denote it by the same symbol $`^2`$.

Let $`\{\stackrel{~}{c}(x)\}`$ is the Bernoulli random field with the parameter

$$c=\mathrm{Pr}\left\{\stackrel{~}{c}\left(x\right)=1\right\}$$

which is said to be the concentration. Below, the tilde marks the randomness of used objects. The field $`\{\stackrel{~}{c}(x)\}`$ induces the random set with realizations $`\{\stackrel{~}{M}:\stackrel{~}{M}^2\}`$ where $`\stackrel{~}{M}=\{x:\stackrel{~}{c}(x)=1\}`$. Sometimes we shall call these realizations configurations of completed vertices or, simply, configurations. This probability distribution for possible realizations $`\stackrel{~}{M}`$ is clearly determined completely by the following probability collection

$$\mathrm{Pr}\{\stackrel{~}{M}:A\stackrel{~}{M}\}=c^{\left|A\right|},A^2.$$

Here and further $`||=\mathrm{Card}\{\}`$.

Naturally, the adjacency relation $`\phi `$ induces the connectedness relation for vertices having included into any configuration $`\stackrel{~}{M}`$. We shall call two vertices $`x`$ and $`y`$ linked in $`\stackrel{~}{M}`$ if there exists a path $`(x_i;i=0,1,2,\mathrm{},n)`$, $`x_i\stackrel{~}{M}`$, $`x_0=x,x_n=y`$ and $`x_i\phi x_{i+1}`$, $`i=0,1,2,\mathrm{},n1`$. The connectedness of vertices of $`\stackrel{~}{M}`$ is the equivalence relation. Therefore, each configuration $`\stackrel{~}{M}`$ is uniquely decomposed on some disjoint equivalence classes $`\stackrel{~}{M}={\displaystyle \underset{j}{}}\stackrel{~}{W}_j`$ generated by the connectedness relation of vertices. These classes we shall name clusters. The collection of clusters related to configurations $`\stackrel{~}{M}`$ will be denoted by $`𝒲[\stackrel{~}{M}]=\{\stackrel{~}{W}_j;j\}`$. If $`x\stackrel{~}{W}_j`$ for a $`j`$ in configuration $`\stackrel{~}{M}`$, then we shall denote this cluster $`\stackrel{~}{W}_j`$ by $`\stackrel{~}{W}(x)`$ .

Let us introduce the random field $`\{\stackrel{~}{a}(x);x^2\}`$ on the base of the set $`\stackrel{~}{M}`$,

$$\stackrel{~}{a}\left(x\right)=\{\begin{array}{cc}1;\hfill & x\stackrel{~}{W}\stackrel{~}{M},\left|\stackrel{~}{W}\left(x\right)\right|<\mathrm{}\text{ ,}\hfill \\ 0;\hfill & \text{otherwise\hspace{0.17em}.}\hfill \end{array}$$

In view that $`\{\stackrel{~}{c}(x)\}`$ is the uniform field (i.e. the probability measure is invariant relative to translations on vectors $`(n_1𝐞_1+n_2𝐞_2)`$, $`n_1,n_2`$), the random field $`\{\stackrel{~}{a}(x)\}`$ is uniform too. Therefore, the probability

$$Q\left(c\right)=\mathrm{Pr}\left\{\stackrel{~}{a}\left(x\right)=1\right\}$$

does not depend on $`x^2`$. This probability is a nondecreasing function of $`c`$ . If $`Q(c)>0`$, then the percolation on the $`\{\stackrel{~}{c}(x)\}`$ is said to exist. In connection with this fact, the following characteristic value is introduced

$$c^{}=sup\{c:Q\left(c\right)=0\}.$$

$`c^{}`$ is called the percolation threshold.

3. Finite clusters on $`^2`$

At construction of some upper estimates of the probability $`Q(c)`$, it is necessary to enumerate all finite clusters $`W(x)`$ including the fixed vertex $`x`$. We perform the enumeration by using the concept of the external boundary for each finite cluster. To this end, let us introduce on $`^2`$ the following new concept of the adjacency relation $`\overline{\phi }`$ side by side with the adjacency relation $`\phi `$ . We denote $`\overline{}^2`$ the corresponding graph having the same vertex set but with the adjacency relation $`\overline{\phi }`$.

Vertices $`x`$ and $`y`$ are named the $`\overline{\phi }`$-adjacent ones if one of two cases takes place: 1) $`x\phi y`$, 2) either $`y=x+𝐞_1\pm 𝐞_2`$ or $`y=x𝐞_1\pm 𝐞_2`$.

The relation $`\overline{\phi }`$ induces the new relation of vertex connectedness on each configuration $`\stackrel{~}{M}`$. This connectedness relation is the equivalence relation too and it leads to a decomposition of $`\stackrel{~}{M}`$ on some connected sets of vertices relative to $`\overline{\phi }`$.

Let us introduce the external boundary concept of a finite cluster on $`^2`$.

D e f i n i t i o n 1. The set $`W`$ is named the boundary of the cluster $`W`$ if it consists of those vertices $`y`$ which are $`\phi `$-adjacent to vertex $`x`$ in the cluster $`W`$ but do not belong to it.

D e f i n i t i o n 2. The external boundary $`\overline{}W`$ of the cluster $`W`$ is the set of vertices $`uW`$ such that for each of them there exists an infinite $`\phi `$-path $`\alpha (u)`$ on completed vertices of $`^2`$ and, moreover, $`u`$ is the unique vertex in $`\alpha (u)`$ belonging to the union $`WW`$.

Further classification of finite clusters $`W(x)`$ is performed by the enumeration of all possible external boundaries $`\overline{}W(x)`$. Following statement is the key for such enumeration. It repeats the corresponding statement in the monograph  with the exception of the last item.

Theorem 1Let $`W(x)`$ be a finite cluster for a fixed vertex $`x`$, $`|W(x)|<\mathrm{}`$. Then $`W(x)`$ has a nonempty finite external boundary $`\overline{}W(x)`$ having following properties.

1$`\overline{}W(x)`$ is the $`\overline{\phi }`$-connected vertex set in $`\overline{}^2`$ which presents the cycle i.e. $`\overline{}W(x)=(x_1,x_2,\mathrm{},x_n)`$ where $`x_ix_j,ij,n=|\overline{}W(x)|,x_i\overline{\phi }x_{i+1},i=1,2,\mathrm{},n`$, $`x_{n+1}=x_0`$ and each vertex has only two $`\overline{\phi }`$-adjacent vertices in the $`\gamma `$. It is possible to introduce the definite orientation on $`\overline{}W(x)`$. (it will be used the counter-clockwise orientation in further arguments.)

2. The vertex $`x`$ is contained in the finite set $`\mathrm{Int}[\overline{}W(x)]`$ defined by

$$\mathrm{Int}\left[\overline{}W\left(x\right)\right]\{u\overline{}W\left(x\right):(\alpha \left(u\right):|\alpha \left(u\right)|=\mathrm{})(\alpha \left(u\right)\overline{}W\left(x\right)\mathrm{})\};$$

3. Let $`u,v,w\overline{}W(x)`$ be three $`\overline{\phi }`$-adjacent vertices following one after another according to the introduced orientation. Then if $`u\phi v`$, the vertex $`w`$ belongs necessarily to one of the following collection:

a) in the case $`u\phi v`$ the set of three elements $`\{2vu,\mathrm{\hspace{0.17em}2}vu\pm 𝐞\}`$ where $`𝐞`$ is a unit basis vector that is orthogonal to the vector $`(vu)`$;

b) in the case $`v=u+\epsilon (𝐞_1𝐞_2)`$ the set of five elements

$$\{2vu,\mathrm{\hspace{0.17em}2}vu+\epsilon 𝐞_1,\mathrm{\hspace{0.17em}2}vu+\epsilon 𝐞_2,\mathrm{\hspace{0.17em}2}vu\pm 𝐞^{}/2\}$$

in the case $`v=u+\epsilon (𝐞_1+𝐞_2)`$ ;

$$\{2vu,\mathrm{\hspace{0.17em}2}vu+\epsilon 𝐞_1,\mathrm{\hspace{0.17em}2}vu\epsilon 𝐞_2,\mathrm{\hspace{0.17em}2}vu\pm 𝐞^{}/2\}$$

where $`𝐞^{}`$ be an orthogonal vector to $`(vu)`$ and $`|𝐞^{}|=\sqrt{2}`$ and $`\epsilon =\pm 1`$.

The proof of first and second statements of Theorem 1 is obvious (see ). What about the last statement, it is easy to verify its justification from the explaining figure 1(a, b). The full rigorous proof is found by us; it is quite tedious and goes beyond the scope of this communication. Unlike the proof of the analogous theorem in , our arguments are not used an application of Jordan’s theorem. Further we omits the proof and pass to the proof of the main statement in next paragraphs.

4. The cluster decomposition on $`^2`$

We shall consider the probability $`1Q(c)=\mathrm{Pr}\{\stackrel{~}{a}(0)=0\}`$. Let $`𝒜=\{W:W=W(0)\text{ is cluster},|W|<\mathrm{}\}`$ be the collection of all finite clusters including the vertex 0. We define the event

$$A\left(W\right)=\{\stackrel{~}{M}:\stackrel{~}{W}=W,\mathrm{\hspace{0.17em}0}\stackrel{~}{M},\stackrel{~}{W}\left\{\stackrel{~}{W}_j;j\right\}\}$$

for any cluster $`W𝒜`$. This event has the definite probability

$$\mathrm{Pr}\left\{A\left(W\right)\right\}=c^{\left|W\right|}\left(1c\right)^{\left|W\right|}.$$
$`\left(1\right)`$
According to statements in the previous section, any cluster of this collection $`𝒜`$ is corresponded to a cycle $`\gamma `$. Vertices of this cycle are $`\overline{\phi }`$-adjacent and such that $`0\mathrm{Int}[\gamma ]`$. In this connection, let us introduce in consideration the collection $``$ of all $`\overline{\phi }`$-cycles having the last property. Let $`B(\gamma )`$ be the event

$$B\left(\gamma \right)=\{\stackrel{~}{M}:\mathrm{\hspace{0.17em}0}\stackrel{~}{M},\stackrel{~}{W}\left(0\right)𝒲\left[\stackrel{~}{M}\right],\overline{}\stackrel{~}{W}\left(0\right)=\gamma \}$$
$`\left(2\right)`$
defined for any $`\overline{\phi }`$-cycle $`\gamma `$. It is represented by the finite union of mutually disjoint events

$$B\left(\gamma \right)=\underset{W𝒜:\overline{}W=\gamma }{}A\left(W\right).$$
$`\left(3\right)`$
Hence, according to (1) and (3), such event has the definite probability

$$P\left(\gamma \right)=\mathrm{Pr}\left\{B\left(\gamma \right)\right\}.$$

It is equal

$$P\left(\gamma \right)=\underset{W𝒜:\overline{}W=\gamma }{}\mathrm{Pr}\left\{A\left(W\right)\right\}=\underset{W𝒜:\overline{}W=\gamma }{}c^{\left|W\right|}\left(1c\right)^{\left|W\right|}.$$

Let us note that

$$\left\{\stackrel{~}{a}\left(0\right)=0\right\}=\underset{W𝒜}{}A\left(W\right).$$
$`\left(4\right)`$
The collection $`𝒜`$ is decomposed on some disjoint classes of clusters. Each class consists of those clusters $`W𝒜`$ which have the same external boundary. It is realized $`\overline{}W=\gamma `$ for them. Therefore, the following representation is correct

$$\underset{W𝒜}{}\mathrm{}=\underset{\gamma }{}\underset{W𝒜:\overline{}W=\gamma }{}\mathrm{}.$$

Then we obtain

$$\left\{\stackrel{~}{a}\left(x\right)=0\right\}=\underset{\gamma }{}B\left(\gamma \right)$$

using (3). Finally, we come to the statement

Theorem 2The probability $`1Q(c)`$ is represented by the decomposition

$$1Q\left(c\right)=\underset{\gamma }{}P\left(\gamma \right).$$
$`\left(5\right)`$
Usually, such decomposition is said to be the cluster one in the percolation theory .

5. The main theorem

The cluster decomposition (5) is represented by the sum of probabilities of some disjoint events. Therefore, definitely, the cluster decomposition is convergent. The function $`Q(c)`$ is not equal to zero only if $`c>c^{}>0`$. So, it is not an analytic function on the concentration $`c`$ and $`c=c^{}`$ is its the singular point. Earlier , it was obtained that $`6/7>^{}>1/3`$. We give the improvement of the upper estimate at the below-formulated statement.

Theorem 3The inequality $`c^{}c_0=3\sqrt{5}`$ is correct for Bernoulli’s random field on $`^2`$.

P r o o f . We use the elementary estimate

$$P\left(\gamma \right)\left(1c\right)^{\left|\gamma \right|},$$

that follows from the Definition 2 and Eq.(2). Using it and Eq.(5), we come to the upper boundary

$$\underset{\gamma }{}P\left(\gamma \right)\underset{\gamma }{}\left(1c\right)^{\left|\gamma \right|}=\underset{\gamma :\left|\gamma \right|=2k}{}\left(1c\right)^{2k}r_k$$
$`\left(6\right)`$
where $`r_k=\mathrm{Card}\{\gamma :|\gamma |=2k\}`$, $`k2`$.

Further, we shall find the upper estimate for the value $`r_k`$. Consider the infinite path $`\alpha (0)=(j𝐞_1;j=0,1,2\mathrm{})`$ with initial vertex $`0`$. Then, according to Theorem 1 (2), each cycle $`\gamma `$ necessarily crosses the path in some vertices. Let us select the vertex in this set of all intersection vertices which is nearest to the vertex $`0`$. Denote it by $`z_\gamma `$. All cycles of $``$ is decomposed on disjoint classes $`𝒞_\gamma `$, i.e. the cycles having the same vertex $`z_\gamma `$ belong to the same class. This class decomposition induces the decomposition of the cycle set $`\{\gamma :|\gamma |=2k\}`$ on corresponding classes $`𝒞_l^{(k)}`$ where $`l`$ is the distance from 0 to $`z_\gamma `$. In addition, $`l<k`$ and therefore,

$$\{\gamma :\left|\gamma \right|=2k\}=\underset{l=1}{\overset{k1}{}}𝒞_l^{\left(k\right)},r_k=\underset{l=1}{\overset{k1}{}}\left|𝒞_l^{\left(k\right)}\right|.$$

Let $`\gamma =(z_\gamma =x_0,x_1,\mathrm{},x_{2k1},x_{2k}=z_\gamma )`$. According to the earlier introduced orientation, the vertex $`x_1`$ that follows after $`x_0`$ in the cycle $`\gamma `$ may be one of the set only

$$(z_\gamma +𝐞_1,z_\gamma +𝐞_1+𝐞_2,z_\gamma +𝐞_2,z_\gamma 𝐞_1+𝐞_2).$$
$`\left(7\right)`$
Using this ordering, cycles belonging to the same class $`𝒞_l^{(k)}`$ are distributed on disjoint collections $`𝒞_l^{(k,i)}`$, $`i=1,2,3,4`$ independence of the selection of the vertex $`x_1`$ in the cycle $`\gamma `$. Hence,

$$𝒞_l^{\left(k\right)}=\underset{i=1,2,3,4}{}𝒞_l^{(k,i)},\left|𝒞_l^{\left(k\right)}\right|=\underset{i=1,2,3,4}{}\left|𝒞_l^{(k,i)}\right|.$$

We obtain as the result

$$r_k=\underset{l=1}{\overset{k1}{}}\underset{i=1,2,3,4}{}\left|𝒞_l^{(k,i)}\right|.$$

We must find upper estimate of the value $`\left|𝒞_l^{(k,i)}\right|`$. It is easy to see that $`𝒞_l^{(k,i)}𝒫_l^{(2k1,i)}`$ where $`𝒫_l^{(2k1,i)}`$ is the set of paths beginning at the vertex $`z_\gamma `$ and the distance from $`z_\gamma `$ to the vertex 0 along the path $`\alpha (0)`$ is equal to $`l`$. Besides, they have the length $`2k1`$ and the vertex $`x_1`$ is the $`i`$-th in the finite sequence (7). We sort out paths $`\gamma `$ with the length $`2k1`$ but not with $`2k`$. It is connected with the fact that the last edge $`x_{2k1}\overline{\phi }x_0`$ of the cycle $`(x_0,x_1,\mathrm{},x_{2k1},x_0)`$ is fixed by the collection of vertices $`x_0,x_1,\mathrm{},x_{2k1}`$. Therefore,

$$\left|𝒞_l^{(k,i)}\right|\left|𝒫_l^{(2k1,i)}\right|=s_{2k1}.$$

It is easy to see that $`s_n`$ does not depend on $`l`$ and $`i`$ due to our construction. Let us represent

$$s_n=s_n^++s_n^\times $$
$`\left(8\right)`$
where

$$s_n^+=\mathrm{Card}\{\gamma 𝒫_l^{(n,i)}:x_{n1}\phi x_n\},$$

$$s_n^\times =\mathrm{Card}\{\gamma 𝒫_l^{(n,i)}:\overline{x_{n1}\phi x_n},x_{n1}\overline{\phi }x_n\}.$$

We introduce the two-component vector $`[s_n^+,s_n^\times ]`$ (not distinguishing column-vectors and tuple-vectors). Then the following equation takes place

$$\left[\begin{array}{c}s_n^+\\ s_n^\times \end{array}\right]=𝖳\left[\begin{array}{c}s_{n1}^+\\ s_{n1}^\times \end{array}\right]$$

with the transfer matrix

$$𝖳=\left[\begin{array}{cc}1& 2\\ 2& 3\end{array}\right].$$

according to Theorem 1 (3). Hence,

$$\left[\begin{array}{c}s_n^+\\ s_n^\times \end{array}\right]=𝖳^{𝗇\mathrm{𝟣}}\left[\begin{array}{c}s_1^+\\ s_1^\times \end{array}\right].$$
$`\left(9\right)`$
Eigenvalues of the matrix $`𝖳`$ is equal to $`\lambda _+=2+\sqrt{5},\lambda _{}=2\sqrt{5}`$ and corresponding eigenvectors (they are orthogonal but are not normalized) have the form

$$𝗍_+=\left[\begin{array}{c}1\\ \frac{1+\sqrt{5}}{2}\end{array}\right],𝗍_{}=\left[\begin{array}{c}1\\ \frac{1\sqrt{5}}{2}\end{array}\right].$$

Decomposing the vector $`[s_1^+,s_1^\times ]=[2,2]`$ on eigenvectors $`𝗍_+,𝗍_{}`$,

$$[s_1^+,s_1^\times ]=g_+𝗍_++g_{}𝗍_{},$$

we obtain $`g_+=1+\sqrt{5}/5`$, $`g_{}=1\sqrt{5}/5`$. Therefore, $`[s_n^+,s_n^\times ]=g_+\lambda _+^{n1}𝗍_++g_{}\lambda _{}𝗍_{}`$ and, according to (8), (9), we find

$$s_n=\left(g_+\lambda _+^{n1}+g_{}\lambda _{}^{n1}\right)+\frac{1}{2}\left(g_+\lambda _+^{n1}\left(1+\sqrt{5}\right)+g_{}\lambda _{}^{n1}\left(1\sqrt{5}\right)\right)$$

$$4\left(2+\sqrt{5}\right)^{n1}.$$
$`\left(10\right)`$
Since $`r_k<4(k1)s_{2k1}`$, hence, putting $`n=2k1`$ at the estimate having found, we obtain the majorant series from (6). It is summable at $`c>3\sqrt{5}`$. Then, if this condition is satisfied, we shall obtain

$$\underset{\gamma }{}P\left(\gamma \right)\mathrm{}.$$
$`\left(11\right)`$
Completion of the proof is performed on the base of arguments which are standard in the percolation theory (see, for example, ). The inequality (11) permits to apply the Borel-Cantelli lemma (see, for example, ) to the event family $``$. According to this statement, the probability of the event that consists in simultaneous realization of an infinite set of events belonging to the family, is equal to zero. Then, with the probability one, there exists a maximal cycle $`\gamma `$. For this cycle, there is a vertex $`z`$ out of it and, besides, there is the infinite path $`\alpha (z)`$ without intersections and with the initial vertex $`z`$. The path does not intersect the cycle. This means that the event $`\{\stackrel{~}{M}:\left(\stackrel{~}{W}𝒲[\stackrel{~}{M}]\right)\left(|\stackrel{~}{W}|=\mathrm{}\right)\}`$ of Bernoulli’s field $`\{\stackrel{~}{c}(x)\}`$ has the probability 1. On the other hand, the countable decomposition

$$\{\stackrel{~}{M}:\left(\stackrel{~}{W}𝒲\left[\stackrel{~}{M}\right]\right)\left(\left|\stackrel{~}{W}\right|=\mathrm{}\right)\}=\underset{v^2}{}\{\stackrel{~}{M}:\stackrel{~}{W}\left(v\right)𝒲\left[\stackrel{~}{M}\right],\left|\stackrel{~}{W}\left(v\right)\right|=\mathrm{}\}$$
$`\left(12\right)`$
takes place. In view of the random field $`\{\stackrel{~}{a}(x)\}`$ is uniform, the probability

$$\mathrm{Pr}\{\stackrel{~}{M}:\stackrel{~}{W}\left(v\right)𝒲\left[\stackrel{~}{M}\right],\left|\stackrel{~}{W}\left(v\right)\right|=\mathrm{}\}=Q\left(c\right)$$

does not depend on $`v`$. Therefore, it cannot be equal to zero, since the following inequality is correct according to (12)

$$1\underset{v^2}{}\mathrm{Pr}\{\stackrel{~}{M}:\stackrel{~}{W}\left(v\right)𝒲\left[\stackrel{~}{M}\right],|\stackrel{~}{W}\left(v\right)|=\mathrm{}\}.\mathrm{}$$

ConsequenceAt $`>c_0`$, the probability $`Q(c)`$ may be represented with any preassigned accuracy by finite sum of series (5). It is defined by summands which correspond to cycles $`\gamma `$ with the length not exceeding $`2m`$ with fixed $`m`$. In this case, the following error estimate takes place

$$\underset{\gamma :\left|\gamma \right|>2m}{}P\left(\gamma \right)\underset{\gamma :\left|\gamma \right|>2m}{}\left(1c\right)^{\left|\gamma \right|}\underset{k>m}{}\left(1c\right)^{2k1}s_{2k1}=$$

$$=16\left(1c\right)^2\underset{k=m+1}{\overset{\mathrm{}}{}}\left(k1\right)\left[\left(2+\sqrt{5}\right)\left(1c\right)\right]^{2k2}=$$

$$=16\left(1c\right)^2\xi \left[\frac{d}{d\xi }\left(\frac{\xi ^m}{1\xi }\right)\right]_{\xi =\left[\left(2+\sqrt{5}\right)\left(1c\right)\right]^2}<\mathrm{}\mathrm{}$$

Thus, the decomposition (5) is the key for solving of the main percolation problem on $`^2`$ that is to calculate the probability $`Q(c)`$ with any guaranteed accuracy. In addition, due to presence of the singular point $`=c^{}`$, the majorant series for the remainder must have a singularity in a point $`{}_{0}{}^{}>c^{}`$ for any initial finite sum of the cluster decomposition. Such a situation will take place in any way of upper estimation of the percolation threshold. Then, it follows that the solution of the main problem is impossible without building of a calculation algorithm for the percolation threshold $`c^{}`$ with any preassigned accuracy.

6. Discussion

One can see on the base of the above proof that the singularity  is possibly connected with the singularity of the generating function for the nonintersecting contours number. If it is true, so the regular method of sequential approximative calculation of this point may be based on the sequential exclusion of intersections in paths $`(z_\gamma ,x_1,\mathrm{},x_{2k1})`$ corresponding to cycles $`\gamma `$. Such an exclusion must be fulfilled with the increasing of approximation order. In this case the estimate given by Theorem 2 may be considered as the zero approximation.

Generally, the algorithm of approximations of the point $`c^{}`$ should be consist in such a construction when estimates $`{}_{n}{}^{}<c^{}<c_n^{}`$, $`c_{n+1}^{}>c_n^{}`$, $`c_{n+1}^+<c_n^+`$ are build sequentially step by step. They have the property $`\underset{n\mathrm{}}{lim}c_n^\pm =c^{}`$ and besides, for sure, $`Q(c)=0`$ at $`c<c_n^{}`$; $`Q(c)>0`$ at $`c>c_n^+`$.

References

1. Yu.P. Virchenko, Percolation. Mathematical Physics. Great Russian Encyclopedia. Moscow (1998), p.346-347 (in Russian)

2. Yu.P. Virchenko, Percolation of random field. Mathematical Physics. Great Russian Encyclopedia. Moscow (1998), p.363-364 (in Russian)

3.  H. Kesten, Percolation Theory for Mathematicians. Boston. Birkhäuser (1982)

4. M.V. Men’shikov, S.A. Molchanov, A.F. Sidorenko, Percolation Theory and Some Its Applications. – In: Advances in Sci. and Engin., Ser. Propability Theory, Mathematical Statistics and Theoretical Cybernetics , Vol. 24. Moscow (1986), VINITI Publ. 53p. (in Russian)

5. V.A. Malyshev, M.V. Men’shikov, E.V. Petrova, Introduction in Probability Theory. Moscow (1997), MGU Publ, p.119 (in Russian)

6. J. Lamperti, Probability. New York - Amsterdam (1966), Dartmuoth college