text:
# Solving large nonlinear generalized eigenvalue problems from Density Functional Theory calculations in parallel

## 1 Introduction

Within recent years it has become possible to perform quantum mechanical calculations from first principles using the Density Functional Theory (DFT) of Hohenberg, Kohn and Sham (see e.g. ). Realistic calculations of the electrons’ ground-state can be carried out for large systems consisting of tens to hundreds of atoms owing to improvements in physical methods, faster computers, and more efficient algorithms. In this way modeling of a vast number of physical and chemical properties can be carried out, often taking large computer systems to their limits.

The Density Functional Theory specifies that the ground-state total energy $`E_0`$ can be obtained as the global minimum of the energy functional $`E[\rho ]`$,

$$E[\rho ]=T[\rho ]+V_{\text{ee}}[\rho ]+\rho (𝐫)v_{\text{ext}}(𝐫)𝑑𝐫+\gamma _{\text{Ewald}}$$
(1)
Here the position vector in space is denoted as $`𝐫`$, $`\rho (𝐫)`$ (or $`\rho `$ for short) is the charge density of electrons satisfying the constraint $`\rho (𝐫)𝑑𝐫=N`$, where $`N`$ is the integral number of electrons in the system. The $`v_{\text{ext}}(𝐫)`$ is the “external” potential function acting on the electrons and originating from the atomic nuclei at fixed positions in the system. The kinetic energy $`T[\rho ]`$ as well as the electron-electron interaction energy $`V_{\text{ee}}[\rho ]`$ are functionals of the density $`\rho `$, only. Finally, an electrostatic repulsion energy between the nuclei, denoted as $`\gamma _{\text{Ewald}}`$, is added. The temperature is taken to be at the absolute zero so that entropy terms are left out in the equations below.

Kohn and Sham proved that the charge density $`\rho (𝐫)`$ can be represented by a system of $`N`$ non-interacting electrons, whose complex-valued quantum-mechanical wavefunctions $`\psi _i(𝐫)`$ for $`i=1,\mathrm{},N`$ yield the charge density as

$$\rho (𝐫)=\underset{i=1}{\overset{N}{}}|\psi _i(𝐫)|^2$$
(2)
(for simplicity, we only consider integral occupation numbers of 0 or 1) and the non-interacting kinetic energy functional $`T_s`$ as

$$T_s[\rho ]=\frac{1}{2}\underset{i=1}{\overset{N}{}}\overline{\psi }_i(𝐫)^2\psi _i(𝐫)𝑑𝐫$$
(3)
($`\overline{\psi }_i(𝐫)`$ denotes the complex conjugate of $`\psi _i(𝐫)`$). The energy functional Eq. 1 can now be rewritten as

$$E[\rho ]=T_s[\rho ]+E_{\text{Hartree}}[\rho ]+E_{\text{xc}}[\rho ]+\rho (𝐫)v_{\text{ext}}(𝐫)𝑑𝐫+\gamma _{\text{Ewald}}$$
(4)
where the $`E_{\text{Hartree}}[\rho ]`$ is the classical electrostatic energy of the electrons, and $`E_{\text{xc}}[\rho ]`$ is an (unknown) exchange-correlation energy functional.

The wavefunctions $`\psi _i(𝐫)`$ are the $`i=1,\mathrm{},N`$ lowest eigensolutions of the Kohn-Sham equation

$$\widehat{H}\psi _i(𝐫)[\frac{1}{2}^2+v_{\text{eff}}(𝐫)]\psi _i(𝐫)=ϵ_i\psi _i(𝐫)$$
(5)
where $`\widehat{H}`$ is known as the Hamiltonian operator. The effective potential $`v_{\text{eff}}(𝐫)`$ is the sum of the external, the Hartree (electrostatic), and the exchange-correlation potentials:

$$v_{\text{eff}}(𝐫)=v_{\text{ext}}(𝐫)+\frac{\delta E_{\text{Hartree}}[\rho ]}{\delta \rho (𝐫)}+\frac{\delta E_{\text{xc}}[\rho ]}{\delta \rho (𝐫)}=v_{\text{ext}}(𝐫)+\frac{\rho (𝐫^{})}{|𝐫𝐫^{}|}𝑑𝐫^{}+v_{\text{xc}}(𝐫)$$
(6)
It is seen from Eq. 6 that the potential $`v_{\text{eff}}(𝐫)`$ depends on the charge density $`\rho (𝐫)`$, which itself is given by Eq. 2 as a sum of squared wavefunctions that are determined as eigensolutions of Eq. 5. Hence a self-consistent solution of these equations must be found starting from an initial guess of $`\rho (𝐫)`$, from which the $`v_{\text{eff}}(𝐫)`$ is constructed; solving the Kohn-Sham eigenvalue equation yields through Eq. 2 a renewed charge density. A minimization algorithm for $`E[\rho ]`$ must be used to obtain the ground-state total energy $`E_0`$ from such iterations.

An additional complication in some of the current approaches is the use of ultrasoft pseudopotentials that puts a lower requirement on the basis set size, but at the price of a more complicated and computationally more demanding Hamiltonian $`\widehat{H}`$. This leads to replacing the Kohn-Sham equation Eq. 5 by a generalized eigenvalue equation,

$$\widehat{H}\psi _i(𝐫)=ϵ_i\widehat{S}\psi _i(𝐫)$$
(7)
where the $`\widehat{S}`$ denotes the overlap projection operator.

In most numerical calculations the wavefunctions $`\psi _i(𝐫)`$ are expanded on a suitable, finite set of basis functions. The present work employs the widely used pseudopotential approximation together with plane-wave (or Fourier expansion) basis sets. In the plane-wave basis the Hamiltonian Eq. 5 is diagonally dominant owing to the kinetic energy term. In this method full, complex Hermitian matrices of sizes $`10^310^5`$ are usually encountered, whereas the number of electrons $`N`$ may typically be two orders of magnitude smaller, $`10^110^3`$. However, the algorithms described below remain valid for very general classes of problems using other basis function sets.

The translational symmetry of crystals is dealt with by introducing a summation over the so-called k-points (see e.g. ). This leads to a set of independent eigenvalue problems for each k-point, which are coupled only through the addition of charge densities in Eq. 2. Similarly, electron spin can double the number of independent eigenvalue problems to be solved. Such almost-decoupled eigenvalue problems should always be solved independently in parallel, since this will be optimally efficient and can easily be combined with the parallelization described below, leading to the possibility of utilizing large numbers of processors efficiently.

The discretized problem originating from Eq. 7, which essentially has to be solved for the smallest $`N`$ eigenvalues, is the non-linear matrix eigenvalue problem,

$`𝐇(\underset{¯}{\rho })\underset{¯}{\psi _i}`$ $`=ϵ_i𝐒\underset{¯}{\psi _i},i=1,\mathrm{},N`$ (8)
$`\underset{¯}{\rho }`$ $`={\displaystyle \underset{i=1}{\overset{N}{}}}\underset{¯}{\psi _i}^H𝐒\underset{¯}{\psi _i},`$ (9)
where the overlap matrix $`𝐒`$ is a symmetric positive definite matrix when ultrasoft pseudopotentials are used (and the identity matrix otherwise). This implies that the eigenvalues, $`ϵ_i`$ are real and the eigenvectors, $`\underset{¯}{\psi _i}`$ mutually $`𝐒`$-orthogonal.

Given an initial charge density vector $`\underset{¯}{\rho }`$ (reasonably good initial values for $`\rho `$ can usually be constructed using the densities of the constituent free atoms), Eq. 8 can be solved for the wave functions and the resulting charge density can subsequently be computed using Eq. 9. A new input charge density can finally be formed, and the cycle can be iterated until convergence. The overall procedure is commonly termed Self Consistent Field (SCF).

Alternative approaches to the SCF method exist, notably direct minimization of the energy as a functional of the wavefunctions $`\psi _i(𝐫)`$ instead of the charge density $`\rho `$, see for example . It has been shown that this method has about the same efficiency as the traditional SCF method. Another active area of research tries to achieve order-N scaling by, e.g., carefully selecting alternative basis functions. However, a discussion of the numerical methods proposed in various order-N methods are beyond the scope of this paper but a recent review can be found in .

## 2 Solving the Non-linear Eigenvalue Problem

Using a quasi newton method for the iterations of Eq. 8 and Eq. 9 generally results in convergence within a small number of iterations. The iterations are usually carried out by the method Pulay, see e.g.  for a review of current methods. We propose a more robust and slightly more efficient approach, a quasi newton method as described in Appendix A. Experiments with the limited memory BFGS method indicated that it has an almost similar efficiency but is less robust.

Solving the eigenvalue problem Eq. 8 is computationally by far the most expensive part of the algorithm because of the cost of the $`𝐇(\underset{¯}{\rho })\underset{¯}{\psi _i}`$ product. The method of Davidson has been especially successful for electronic structure computations, and has been improved by iterating on several eigenpairs simultaneously (see e.g. ) and using better preconditioners. Lately, better theoretical understanding of the methods has contributed to further generalizations such as the Jacobi-Davidson method and restarting techniques along with other iterative methods for large eigenvalue problems based on Rayleigh quotient iterations, inverse iteration or the Lanczos method. A currently advocated method in this field is the residual minimization method, RMM-DIIS but as this method essentially just computes eigenpairs closest to the initial ones and therefore easily can result in eigenpairs being missed, we cannot recommend this approach.

A major drawback—in the present context—of the above methods is that they all (except for RMM-DIIS) focus on solving the eigenvalue problem to a fairly high accuracy. However, this turns out not to be necessary during every step of the selfconsistency cycle; in fact, our calculations have shown that solving the eigenvalue problem to a high accuracy generally does not decrease the number of overall steps in the quasi newton/fixpoint iteration noticeably, and in certain cases even increases the number of steps.

Thus we seek an eigenproblem solver which, given an initial estimate of $`N`$ eigenpairs, is able to improve this estimate (i.e. perform a relatively small number of iterations) for each of the eigenpairs. In addition, an important aspect of large-scale applications is that the memory requirement must be minimized. In order to achieve these objectives, we investigated a number of algorithms for improving eigenpairs so that residuals are improved only until a certain, adjustable limit.

Using deflation techniques we experienced convergence problems when deflating eigenpairs which were not well converged, and therefore the overall eigenvalue problem had to be solved to a much higher accuracy (and at a higher computational cost) than actually required. Similar experiences have been reported in , where the convergence tolerance was loosened as more and more eigenpairs converged.

We instead propose an algorithm for solving the generalized eigenvalue problem Eq. 8, which essentially is a blocked Davidson-like algorithm using a non-orthogonal basis, but where a maximum number of expansions on each eigenpair can be imposed.

The reason for choosing a non-orthogonal basis, $`𝐁`$ (and thus obtaining a generalized eigenvalue problem in the projected space, step 1 and 3iC) is, besides from reducing the computational cost compared to using an $`𝐒`$-orthonormal basis, that it reduces the number of synchronization points in the parallel implementation (even when compared to an orthogonalization algorithm implemented with delayed summation as in e.g. ). The drawback is that projected eigenvalue problem, step 1 and 3iC could be poorly conditioned and thus slow down the convergence\[17, Sec. 11.10\]. Our initial basis (step 1) is however well conditioned for every quasi newton iteration (as the estimate of eigenvectors $`\underset{¯}{\psi }_1,\mathrm{},\underset{¯}{\psi }_N`$ is almost $`𝐒`$-orthogonal and $`𝐒`$ is close to the identity) and as only a few expansions are performed for each eigenpair (typically 2 to 3) a poorly conditioned (projected) eigenvalue problem is unlikely to build up to the point where it becomes a problem. This is due to the fact that the eigenvectors do not have to be converged to a very high accuracy. In fact, introducing orthogonalization in the algorithm has not led to fewer iterations for a wide variety of tested physical problems.

The algorithm is as follows:

1. Given an estimate $`\underset{¯}{\psi }_1,\mathrm{},\underset{¯}{\psi }_N`$ for the $`N`$ eigenvectors, initialize the subspace basis $`𝐁[\underset{¯}{b}_1,\mathrm{},\underset{¯}{b}_N]=[\underset{¯}{\psi }_1,\mathrm{},\underset{¯}{\psi }_N]`$ and solve the small, projected subspace eigenvalue problem of dimension $`N`$:

$$𝐁^H\mathrm{𝐇𝐁}\underset{¯}{\varphi }_i=\lambda _i𝐁^H\mathrm{𝐒𝐁}\underset{¯}{\varphi }_i,i=1,\mathrm{},N,\text{ where }\lambda _1\lambda _2\mathrm{}\lambda _N.$$

Select an iteration block-size $`n_b`$ and a maximum number of basis vectors $`n_{\text{max}}`$ and a maximum number of iterations on each eigenpair $`k_{\text{max}}`$.
2. Initialize the set $`𝒦=\mathrm{}`$.
3. Loop over all eigenpairs, $`i=1,\mathrm{},N`$

1. If the $`i`$’th eigenpair $`\lambda _i,\underset{¯}{\varphi }_i`$ has less than $`k_{\text{max}}`$ iterations then compute the corresponding residual $`\underset{¯}{r}_i=\mathrm{𝐇𝐁}\underset{¯}{\varphi }_i\lambda _i\mathrm{𝐒𝐁}\underset{¯}{\varphi }_i`$. If $`\underset{¯}{r}_i`$ is larger than some tolerance, add $`i`$ to $`𝒦`$.
    2. If $`𝒦`$ has $`n_b`$ elements or $`i=N`$:

    1. Do an iteration on eigenpairs given by indexes in $`𝒦`$:

        1. “Precondition” the residuals111 The usual physical approximation in  is used, where the Hamiltonian is approximated by the diagonal (and diagonally dominant) kinetic energy term in Eq. 5. A smooth function is multiplied onto the kinetic term in order to damp the high-frequency errors.  corresponding to the unconverged eigenpairs, i.e. solve approximately for $`i𝒦:(𝐇\lambda _i𝐒)\underset{¯}{x}_i=\underset{¯}{r}_i`$.
            2. Extend the basis, $`i𝒦:𝐁=[𝐁,\underset{¯}{x}_i/\underset{¯}{x}_i]`$.
            3. Solve the updated subspace eigenproblem for the $`N`$ smallest eigenvalues,

            $$\begin{array}{c}𝐁^H\mathrm{𝐇𝐁}\underset{¯}{\varphi }_i=\lambda _i𝐁^H\mathrm{𝐒𝐁}\underset{¯}{\varphi }_i,i=1,\mathrm{},N,\hfill \\ \hfill \text{ where }\lambda _1\mathrm{}\lambda _N.\end{array}$$
        2. If the number of basis vectors in $`𝐁`$ is too large (i.e. larger than $`n_{\text{max}}n_b`$) then collapse the subspace (restart):

        $$𝐁=𝐁[\underset{¯}{\varphi }_1,\mathrm{},\underset{¯}{\varphi }_N].$$
        3. Set $`𝒦=\mathrm{}`$.
4. Update $`[\underset{¯}{\psi }_1,\mathrm{},\underset{¯}{\psi }_N]=𝐁[\underset{¯}{\varphi }_1,\mathrm{},\underset{¯}{\varphi }_N]`$

The tolerance of step 3a is set to an order of magnitude smaller than the current residual norm of the quasi newton process (i.e. $`\underset{¯}{r}^k_2`$, page A), but no larger than $`0.1`$. The tolerance will therefore decrease with the convergence of the quasi newton iteration.

Since the most time consuming part of the eigenvalue solver is computing the products of $`𝐇`$ and $`𝐒`$ by vectors in the steps 1, 3a and 3iC these vectors are stored and reused whenever possible. Thus the dominant memory requirement of the algorithm is $`𝒪(3nn_{\text{max}})`$, where $`n`$ is the vector length given by the dimension of Eq. 8.

Typical problem sizes might be $`N10^110^3`$, and as it will be shown in Sec. 4 the parameters should typically be chosen so $`k_{\text{max}}1\text{ to }2`$, $`n_b10`$ and $`n_{\text{max}}N+5n_b`$, so the memory requirement is very close to the absolute minimum, which is the $`3N`$ vectors $`\underset{¯}{\psi _i}`$, $`𝐇\underset{¯}{\psi _i}`$ and $`𝐒\underset{¯}{\psi _i}`$ for $`i=1,\mathrm{},N`$. Lower memory consumption could be achieved by eliminating the storage of $`𝐇\underset{¯}{\psi _i}`$ and $`𝐒\underset{¯}{\psi _i}`$ and recomputing these products as required in steps 3a and 3iC, however with the tradeoff of requiring about twice as many multiplications of $`𝐇`$ and $`𝐒`$ by vectors.

Most floating-point operations stem from the very complicated products of $`𝐇`$ and $`𝐒`$ by vectors (the number of such operations should be minimized by any algorithm). In addition, the present algorithm contains significant amounts of linear algebra (approximately $`8nN`$ flops for each $`𝐇`$ product) which however can be carried out efficiently using BLAS level 3 operations.

The asymptotic scaling of the total time spent on the matrix-products for each quasi newton iteration is $`Nnlog(n)`$ due to the FFT-technique used, whereas the linear algebra scales as $`N^2n`$. Since the matrix size $`n`$ increases roughly linearly with $`N`$, it is obvious that the linear algebra would dominate asymptotically. For the sizes of our problems, however, the matrix products generally take longer than the linear algebra.

## 3 Parallel Implementation

The products of $`𝐇`$ times vector parallelize particularly poorly, since they involve 3D FFTs of fairly small size (typically $`50`$) in each dimension. Therefore the parallelization for a distributed memory computer is not viable solely across the plane waves (vectors of length $`n`$). By using a block size $`n_b`$ equal to the number of available processors, it is however possible to perform $`n_b`$ $`𝐇`$-products in parallel using local FFTs exclusively – each on a different processor (the $`n_b`$ may also be a multiple of the number of processors, requiring additional memory, however). This approach is used in the implementation for the initial calculation of the $`𝐇`$ and $`𝐒`$ products in step 1 and for step 1, where the $`n_b`$ residual vectors each are gathered on different processors before preconditioning them. After the preconditioning, $`𝐇`$ and $`𝐒`$ products are formed for the preconditioned vector, and the results are finally redistributed across the processors. This requirement on the block size $`n_b`$ can be lifted on computers that perform small-size parallel 3D FFTs with high efficiency.

For the remaining part of the code the parallelization is performed across the plane waves (i.e. each basis-vector is partitioned among the processors). This requires a fairly minimal amount of communication (essentially only a small number of reduction operations).

The parallel implementation was done using MPI, and some typical results indicating the scaling for runs carried at on the IBM SP222Located at UNI$``$C in Lyngby, Denmark. are given in Table 1333 This Platinum (Pt) problem is a an infinite slab with 35 atoms in the unit cell. The parameters used for the run were: vector length $`n=16409`$ (25 Ry plane-wave cutoff energy), $`N=210`$, $`n_b=\text{\# of processors}`$, $`n_{\text{max}}=N+8n_b`$, number of selfconsistent iterations = 4. The execution time is given for the eigensolver but with time for the complete runs including initialization (not yet fully parallelized) given in parenthesis..

It is seen that the speed-up of the eigenvalue solver is very satisfactory for these problem sizes, and larger problems are expected to scale up to a correspondingly higher number of processors. Notice that the superlinear speed-up is partly due to the fact that the blocksize and maximum subspace dimension varies with the number of processors which implies that the 1 processor run did not benefit from BLAS level 3 operations. A sequential run with the same parameters as the 4 processor parallel run led to a speed-up of 3.0 on 4 processors.

These runs were made using only a single “k-point” (see Introduction). In production runs the additional trivial parallelism of a number of k-points and spins are utilized. The parallelization of the eigenvalue solver is duplicated across mutually exclusive sets of processors for each k-point and spin in an optimally efficient way. Our work uses MPI communicators to implement this in a simple way. However, memory consumption is not reduced by exploiting this extra parallel dimension.

## 4 Numerical Examples

In order to investigate the stability of the proposed algorithm and its sensitivity to choice of parameters, a large number of different examples have been tested. Some representative examples444The Gold (Au) problem shown in the figures is an infinite slab with 18 atoms in the unit cell. The parameters used for the run were: vector length $`n=13665`$ (30 Ry plane-wave cutoff energy), $`N=198`$, number of selfconsistent iterations = 15.  are shown in Fig. 1, 2 and 3.

The experiments indicate that the size of $`n_b`$ is not critical as long as it is chosen large enough for the solution of the projected eigensystem not to be unnecessarily large. The maximum subspace dimension, $`n_{\text{max}}`$ is however important – if $`n_{\text{max}}`$ is too small the overhead in restarting becomes dominant, on the other hand if $`n_{\text{max}}`$ becomes unnecessarily large the memory consumption and the cost of solving the projected eigenvalue problem becomes large.

From Fig. 3 it is seen that the overall work grows linearly when increasing the precision to which we solve the eigenvalue problem, indicating that it is not worthwhile to solve the eigenvalue problem to a high accuracy. This is furthermore supported by Fig. 4 which shows that the convergence of the quasi newton process is largely unaffected by increasing the number of iterations on each eigenpair, $`k_{\text{max}}`$.

## 5 Conclusion

The algorithm proposed in this paper presents an efficient way to solve the nonlinear generalized eigenvalue problems occuring in DFT calculations on multiple processors. Part of the algorithm is a variant of the Davidson algorithm which focuses on using only small amounts of work for each eigenpair in each of the quasi newton iterations, and thereby being able to use a non-orthogonal basis.

In contrast, our best effort using traditional deflation techniques typically required 3 to 5 times the amount of overall work, as each of the eigenpairs had to be converged to a much higher precision in order for the deflation to be robust. Attempts were made to improve on this and an algorithm using deflation which was almost as efficient as the algorithm proposed on page 1 was devised. However, it involved complicated heuristics especially for treating eigenvalues with higher multiplicity (which often occur due to symmetries in the systems) and this additional effort did not seem worthwhile.

In summary, the algorithm proposed in the present paper has been designed for nonlinear eigenvalue problems such as those occurring in selfconsistent DFT calculations. The algorithm has shown robustness and efficiency, is flexible in accomodating memory limitations, and allows constraints in the number of iterations and the tolerance of eigenpairs. Several parameters can be used to optimize the convergence, and we propose recommended values for DFT iterations. Good parallel performance is achieved on distributed memory computers, but since the eigenvalue problem in the projected subspace is not parallelized, the parallel scalability may be limited, independent of the efficiency of the interconnecting network.

## Appendix A The Projected Newton Method

Given a series of input charge densities, $`\underset{¯}{\rho }_{\text{in}}^{km},\mathrm{},\underset{¯}{\rho }_{\text{in}}^k`$ ($`km1`$) for Eq. 8, the corresponding output charge densities defined by Eq. 9 will be denoted as $`\underset{¯}{\rho }_{\text{out}}^{km},\mathrm{},\underset{¯}{\rho }_{\text{out}}^k`$. The residual of the charge density iteration is then defined by

$$\underset{¯}{r}^k=\underset{¯}{\rho }_{\text{out}}^k\underset{¯}{\rho }_{\text{in}}^k$$

of which the Jacobian can be approximated by $`𝐓_k𝐆_k^{}`$, where $`𝐆^{}`$ denotes a pseudo inverse of $`𝐆`$, and we have defined

$`𝐆_k`$ $`=(\underset{¯}{\mathrm{\Delta }\rho }^{km+1},\mathrm{},\underset{¯}{\mathrm{\Delta }\rho }^k),`$ $`\underset{¯}{\mathrm{\Delta }\rho }^k`$ $`=\underset{¯}{\rho }_{\text{in}}^k\underset{¯}{\rho }_{\text{in}}^{k1}`$
$`𝐓_k`$ $`=(\underset{¯}{\mathrm{\Delta }r}^{km+1},\mathrm{},\underset{¯}{\mathrm{\Delta }r}^k),`$ $`\underset{¯}{\mathrm{\Delta }r}^k`$ $`=\underset{¯}{r}^k\underset{¯}{r}^{k1}.`$
Applying Newton’s method gives

$$\underset{¯}{\rho }_{\text{in}}^{k+1}=\underset{¯}{\rho }_{\text{in}}^k𝐆_k𝐓_k^{}\underset{¯}{r}^k$$
(10)
where $`\underset{¯}{\alpha }_k𝐓_k^{}\underset{¯}{r}^k`$ is obtained as the least squares solution of

$$𝐓_k\underset{¯}{\alpha }_k=\underset{¯}{r}^k.$$
(11)
It is seen that the fixpoint solution $`\underset{¯}{\rho }^{}\underset{¯}{\rho }^k`$ corresponding to $`\underset{¯}{r^k}=\underset{¯}{0}`$ will never be found unless $`\underset{¯}{\rho }^{}\underset{¯}{\rho }^i`$ (where $`i`$ is equal to the initial $`k`$ used) belongs to the span of $`𝐆_i`$ (i.e., the span of $`𝐆_k`$ remains constant). Therefore, a term $`\beta (\underset{¯}{r}^k𝐓_k\underset{¯}{\alpha }_k)`$ is added to the right hand side of Eq. 10 for some arbitrarily chosen full rank, diagonal matrix $`\beta `$ as proposed in  (notice that $`\underset{¯}{r}^k=𝐓_k\underset{¯}{\alpha }_k`$ if $`𝐓_k`$ has full rank).

The least squares system Eq. 11 is solved via the thin QR factorization on the augmented matrix $`[𝐓_k\underset{¯}{\alpha }_k]`$ as described in e.g. , and the factorization is updated rather than recomputed in the following iteration using the routines from .

As the above algorithm can not be used for the first iteration, the value $`\underset{¯}{\rho }_{\text{in}}^1=\underset{¯}{\rho }_{\text{in}}^0+\beta \underset{¯}{r}^0`$ is used instead. For the subsequent iterations the above procedure is applied with $`m`$ incremented by one for each iteration (starting with $`m=1`$) until a maximum value $`m_{\text{max}}`$ is reached. Choosing $`m_{\text{max}}`$ in the range of $`5`$ to $`10`$ generally works well for the present DFT problems.


format: text

text:
# Calorons in Weyl Gauge

## Abstract

We demonstrate by explicit construction that while the untwisted Harrington-Shepard caloron $`A_\mu `$ is manifestly periodic in Euclidean time, with period $`\beta =\frac{1}{T}`$, when transformed to the Weyl ($`A_0=0`$) gauge, the caloron gauge field $`A_i`$ is periodic only up to a large gauge transformation, with winding number equal to the caloron’s topological charge. This helps clarify the tunneling interpretation of these solutions, and their relation to Chern-Simons numbers and winding numbers.

PACS: 11.15.-q, 11.10.Wx, 11.27.+d

Instantons have many important applications in particle physics . They may be defined as classical solutions to the Euclidean Yang-Mills equations that minimize the Yang-Mills action within a given topological charge sector. Explicit instanton solutions  can be found by solving the first-order self-duality equation

$$F_{\mu \nu }=\frac{1}{2}ϵ_{\mu \nu \rho \sigma }F_{\rho \sigma }$$
(1)
At zero temperature, instantons provide a semiclassical description of tunneling between classical minima of the Yang-Mills potential that are connected by large gauge transformations . At finite temperature, the simplest instanton solutions are the untwisted ’calorons’ of Harrington and Shepard . These calorons are solutions to the self-duality equation (1) that are periodic in Euclidean time, with period $`\beta =\frac{1}{T}`$ equal to the inverse temperature. At first sight, it may seem puzzling that a periodic solution, with $`A_\mu (\stackrel{}{x},t=\frac{\beta }{2})=A_\mu (\stackrel{}{x},t=+\frac{\beta }{2})`$, could describe tunneling from one sector to another since the gauge field $`A_\mu `$ is the same at either end. The answer is very simple but illustrative, and does not appear to have been spelled out in detail before. In the Weyl ($`A_0=0`$) gauge, which is most convenient for describing Hamiltonian processes such as tunneling, the coordinate gauge field $`A_i`$ is not periodic; rather, $`A_i(\stackrel{}{x},t=+\frac{\beta }{2})`$ is related to $`A_i(\stackrel{}{x},t=\frac{\beta }{2})`$ by a fixed-time large gauge transformation whose winding number is equal to the integer topological charge of the caloron solution. Here we present a simple explicit construction of this gauge transformation, yielding the Weyl gauge form of the caloron. We conclude with some comments concerning the distinction between the caloron solutions and the so-called “periodic instantons” studied in . For simplicity, we concentrate in this paper on untwisted calorons , which are related to tunneling between sectors of different Chern-Simons charge; we will not address twisted calorons , which involve tunneling between different flux sectors. Twisted calorons are important for Yang-Mills-Higgs theory with symmetry breaking, and especially for understanding the monopole content of the theory. The twisted calorons have nontrivial holonomy, characterized by a nontrivial Polyakov loop at spatial infinity: $`𝒫\mathrm{exp}\left(_0^\beta 𝑑tA_0\right)`$, and the expectation value of $`A_0`$ at infinity plays the role of a Higgs expectation value . Thus, it is simpler to discuss first the Weyl gauge for untwisted calorons which have a trivial holonomy.

In discussing SU(2) instantons, it is convenient to use a radially symmetric ansatz form . Here, we adopt the conventions of  for the symmetric ansatz for the SU(2) instanton field $`A_\mu =\frac{\sigma ^a}{2i}A_\mu ^a`$ :

$`A_i^a`$ $`=`$ $`{\displaystyle \frac{\varphi _21}{r^2}}ϵ_{iak}x_k+{\displaystyle \frac{\varphi _1}{r}}\left(\delta _{ia}{\displaystyle \frac{x_ix_a}{r^2}}\right)+A_1{\displaystyle \frac{x_ix_a}{r^2}}`$ (2)
$`A_0^a`$ $`=`$ $`A_0{\displaystyle \frac{x_a}{r}}`$ (3)
where $`A_0`$, $`A_1`$, $`\varphi _1`$, and $`\varphi _2`$ are functions of $`r=\sqrt{\stackrel{}{x}^2}`$ and $`t`$ only. The charge 1 instanton for SU(2) can then be expressed in terms of a single function $`\rho (r,t)`$, with the ansatz functions in (3) given by:

$`A_0`$ $`=`$ $`_r\mathrm{ln}\rho `$ (4)
$`A_1`$ $`=`$ $`_0\mathrm{ln}\rho `$ (5)
$`\varphi _1`$ $`=`$ $`r_0\mathrm{ln}\rho `$ (6)
$`\varphi _2`$ $`=`$ $`1+r_r\mathrm{ln}\rho `$ (7)
Then the gauge field described by (3,7) satisfies the instanton equation (1) if $`\rho `$ satisfies the linear equation: $`\mathrm{}\rho =0`$. For the zero temperature, charge 1, instanton in SU(2), the function $`\rho (r,t)`$ can be chosen :

$`\rho =1+{\displaystyle \frac{\lambda ^2}{r^2+t^2}}`$ (8)
This form corresponds to a single instanton located at the origin in $`R^4`$, and with scale parameter $`\lambda `$. This choice is known as the ‘singular gauge’.

To facilitate the transformation of the gauge field in (3) into the Weyl ($`A_0=0`$) gauge, we use the fact  that the ansatz fields $`A_0`$, $`A_1`$, $`\varphi _1`$ and $`\varphi _2`$, which appear in (3) and (7), transform in an abelian-like manner under a particular type of SU(2) gauge transformation on the non-abelian gauge field $`A_\mu =\frac{\sigma ^a}{2i}A_\mu ^a`$. Namely, if the gauge transformation matrix $`U`$ has the special form

$`U(\stackrel{}{x},t)=e^{\frac{i}{2}f(r,t)\widehat{x}\stackrel{}{\sigma }}=\mathrm{cos}({\displaystyle \frac{f}{2}})i(\widehat{x}\stackrel{}{\sigma })\mathrm{sin}({\displaystyle \frac{f}{2}})`$ (9)
then the SU(2) gauge transformation, $`A_\mu \stackrel{~}{A}_\mu =U^1A_\mu U+U^1_\mu U`$, has the following effect on the ansatz functions $`A_0`$, $`A_1`$, $`\varphi _1`$ and $`\varphi _2`$ :

$`A_0`$ $``$ $`\stackrel{~}{A_0}=A_0+_0f`$ (10)
$`A_1`$ $``$ $`\stackrel{~}{A_1}=A_1+_rf`$ (11)
$`\varphi _1`$ $``$ $`\stackrel{~}{\varphi }_1=\mathrm{cos}f\varphi _1+\mathrm{sin}f\varphi _2`$ (12)
$`\varphi _2`$ $``$ $`\stackrel{~}{\varphi }_2=\mathrm{sin}f\varphi _1+\mathrm{cos}f\varphi _2`$ (13)
Thus, to achieve the Weyl gauge, $`\stackrel{~}{A}_0^a=0`$, we simply need a gauge transformation $`U(\stackrel{}{x},t)`$ of the form in (9) where the function $`f(r,t)`$ satisfies

$$_0f=_r\mathrm{ln}\rho $$
(14)
For the zero temperature, charge 1, instanton with $`\rho `$ given by (8), this leads to

$`f(r,t)=\left(\pi +2\mathrm{arctan}{\displaystyle \frac{t}{r}}\right){\displaystyle \frac{r}{\sqrt{r^2+\lambda ^2}}}\left(\pi +2\mathrm{arctan}{\displaystyle \frac{t}{\sqrt{r^2+\lambda ^2}}}\right)`$ (15)
where we have chosen the constant of integration so that $`f(r,t=\mathrm{})=0`$. The arctan function is understood to take values in the $`[\frac{\pi }{2},\frac{\pi }{2}]`$ branch. The gauge field, $`\stackrel{~}{A}_i^a`$, in Weyl gauge has the ansatz form (3), with ansatz functions $`\stackrel{~}{\varphi }_1`$, $`\stackrel{~}{\varphi }_2`$, $`\stackrel{~}{A}_1`$, (and $`\stackrel{~}{A}_00`$) given by (13).

The tunneling interpretation of the instanton solution becomes evident by noting that $`\stackrel{~}{A}_i`$ is pure gauge at $`t=\pm \mathrm{}`$. At $`t=\mathrm{}`$,

$`\stackrel{~}{A}_i^a(\stackrel{}{x},t=\mathrm{})=0`$ (16)
while at $`t=+\mathrm{}`$,

$`\stackrel{~}{A}_i^a(\stackrel{}{x},t=+\mathrm{})`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}\left(\frac{2\pi r}{\sqrt{r^2+\lambda ^2}}\right)1}{r^2}}ϵ_{iak}x_k{\displaystyle \frac{\mathrm{sin}\left(\frac{2\pi r}{\sqrt{r^2+\lambda ^2}}\right)}{r}}\left(\delta _{ia}{\displaystyle \frac{x_ix_a}{r^2}}\right){\displaystyle \frac{2\pi \lambda ^2}{(r^2+\lambda ^2)^{3/2}}}{\displaystyle \frac{x_ix_a}{r^2}}`$ (17)
$`=`$ $`\left(W^1_iW\right)^a`$ (18)
where the static gauge function $`W(\stackrel{}{x})`$ is

$`W(\stackrel{}{x})=e^{\frac{i}{2}f(r,t=+\mathrm{})\widehat{x}\stackrel{}{\sigma }}=\mathrm{exp}\left[i{\displaystyle \frac{\pi r}{\sqrt{r^2+\lambda ^2}}}\widehat{x}\stackrel{}{\sigma }\right]`$ (19)
Thus, as $`t`$ goes from $`t=\mathrm{}`$ to $`t=+\mathrm{}`$, the instanton gauge field $`\stackrel{~}{A}_i`$ in Weyl gauge interpolates between the pure gauge $`\stackrel{~}{A}_i=0`$ in (16), and the pure gauge $`\stackrel{~}{A}_i=W^1_iW`$ in (18). These fields $`\stackrel{~}{A}_i`$ are neighboring minima (in fact, zeros) of the classical Yang-Mills potential

$`V={\displaystyle \frac{1}{2}}{\displaystyle d^3xB_i^aB_i^a}`$ (20)
¿From (16) and (18) we see that these two minima are related by the large gauge transformation $`W(\stackrel{}{x})`$ in (19). Furthermore, the fixed-time gauge transformation in (19) is $`W(\stackrel{}{x})=U(\stackrel{}{x},t=+\mathrm{})`$, where $`U(\stackrel{}{x},t)`$ is the (time-dependent) gauge transformation used to transform to the Weyl gauge. To evaluate the winding number of the fixed-time gauge transformation $`W`$, we note that if a gauge transformation has the form $`W(\stackrel{}{x})=e^{\frac{i}{2}g(r)\widehat{x}\stackrel{}{\sigma }}`$, then the winding number $`N`$ is determined  by the values of $`g(r)`$ at $`r=0`$ and $`r=\mathrm{}`$ :

$`N{\displaystyle \frac{1}{24\pi ^2}}{\displaystyle d^3xϵ_{ijk}\mathrm{tr}\left(W^1_iWW^1_jWW^1_kW\right)}={\displaystyle \frac{1}{2\pi }}\left[g(r)\mathrm{sin}g(r)\right]_{r=0}^{r=\mathrm{}}`$ (21)
Thus, for the gauge transformation in (19), we find $`N=1`$.

Having reviewed this familiar tunneling interpretation of zero temperature instantons, we now turn to finite temperature instantons, or ‘calorons’. At finite T, the Euclidean time direction is compactified onto the finite interval $`t[\frac{\beta }{2},+\frac{\beta }{2}]`$, where $`\beta =\frac{1}{T}`$ is the inverse temperature. The manifold is now $`R^3\times S^1`$, rather than $`R^4`$. Instantons are solutions to the same self-duality equation (1), but the gauge fields, being bosonic, satisfy periodic boundary conditions  on the Euclidean time interval $`[\frac{\beta }{2},+\frac{\beta }{2}]`$.

The untwisted Harrington-Shepard (HS) caloron solution is consistent with the ansatz form described in (3) and (7). The only change is that the function $`\rho (r,t)`$ becomes

$`\rho (r,t)`$ $`=`$ $`1+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{2\pi \lambda }{\beta }}\right)^2\left({\displaystyle \frac{\beta }{2\pi r}}\right){\displaystyle \frac{\mathrm{sinh}\left(\frac{2\pi r}{\beta }\right)}{\mathrm{cosh}\left(\frac{2\pi r}{\beta }\right)\mathrm{cos}\left(\frac{2\pi t}{\beta }\right)}}`$ (22)
$``$ $`1+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\overline{\lambda }^2}{\overline{r}}}{\displaystyle \frac{\mathrm{sinh}\overline{r}}{\mathrm{cosh}\overline{r}\mathrm{cos}\overline{t}}}`$ (23)
where we have defined the convenient short-hand: $`\overline{r}=\frac{2\pi r}{\beta }`$, etc. Note that $`\rho `$ in (23) satisfies $`\mathrm{}\rho =0`$. The caloron solution given by (3), (7) and (23) corresponds to an infinite line of zero temperature instantons, each of scale $`\lambda `$ and each at spatial position $`\stackrel{}{x}=0`$, lined up along the t axis, with period $`\beta `$. It is straightforward to check that in the zero temperature limit ($`\beta \mathrm{}`$), the expression for $`\rho (r,t)`$ in (23) reduces smoothly to that in (8), and we recover the zero temperature instanton. At finite T, the caloron solution (3)-(7) is manifestly periodic in $`t`$, $`A_\mu (\stackrel{}{x},t+\beta )=A_\mu (\stackrel{}{x},t)`$, since $`\rho (r,t)`$ in (23) has period $`\beta `$.

To explore the tunneling interpretation of this caloron solution, we transform it to the Weyl gauge. Once again, the Weyl gauge can be attained most easily by using the abelian-like transformation properties (13) of the ansatz functions in (3). The required SU(2) transformation matrix has the form in (9), with $`f(r,t)`$ satisfying (14), and $`\rho `$ is now given by (23). The solution to (14) in this case is

$`f(r,t)`$ $`=`$ $`\left(\pi +2\mathrm{arctan}\left[\mathrm{coth}({\displaystyle \frac{\overline{r}}{2}})\mathrm{tan}({\displaystyle \frac{\overline{t}}{2}})\right]\right)+`$ (25)
$`{\displaystyle \frac{(\overline{\lambda }^22\overline{r}^2)\mathrm{sinh}\overline{r}\overline{\lambda }^2\overline{r}\mathrm{cosh}\overline{r}}{\overline{r}\sqrt{\mathrm{sinh}\overline{r}(4\overline{\lambda }^2\overline{r}\mathrm{cosh}\overline{r}+(\overline{\lambda }^4+4\overline{r}^2)\mathrm{sinh}\overline{r})}}}\left(\pi +2\mathrm{arctan}\left[\mathrm{tan}({\displaystyle \frac{\overline{t}}{2}}){\displaystyle \frac{\overline{\lambda }^2\mathrm{sinh}\overline{r}+2\overline{r}(1+\mathrm{cosh}\overline{r})}{\sqrt{\mathrm{sinh}\overline{r}(4\overline{\lambda }^2\overline{r}\mathrm{cosh}\overline{r}+(\overline{\lambda }^4+4\overline{r}^2)\mathrm{sinh}\overline{r})}}}\right]\right)`$
where the constant of integration has been chosen so that $`f(r,t=\frac{\beta }{2})=0`$. Note that this function $`f(r,t)`$, and hence the gauge transformation $`U(\stackrel{}{x},t)=\mathrm{exp}(\frac{i}{2}f(r,t)\widehat{x}\stackrel{}{\sigma })`$ that transforms the HS caloron to the Weyl gauge, is manifestly periodic in t, with period $`\beta `$. It is straightforward to check that $`f(r,t)`$ in (25) reduces smoothly, as $`\beta \mathrm{}`$, to the zero temperature function in (15).

Given the function $`f(r,t)`$ in (25), the caloron fields $`\stackrel{~}{A}_i^a`$ in the Weyl gauge are obtained simply from (3) and (13). To investigate the relation of this caloron solution to quantum tunneling, we evaluate $`\stackrel{~}{A}_i^a`$ at $`t=\pm \frac{\beta }{2}`$. At $`t=\frac{\beta }{2}`$, the ansatz fields satisfy $`\stackrel{~}{A}_1=\stackrel{~}{\varphi }_1=0`$, while $`\stackrel{~}{\varphi }_2=[1+r_r\mathrm{ln}\rho (r,t=\frac{\beta }{2})]`$, where $`\rho (r,t=\frac{\beta }{2})=1+\frac{\overline{\lambda }^2}{2\overline{r}}\mathrm{tanh}(\frac{\overline{r}}{2})`$. Thus,

$`\stackrel{~}{A}_i^a(\stackrel{}{x},t={\displaystyle \frac{\beta }{2}})={\displaystyle \frac{\lambda ^2\pi \left(\mathrm{tanh}(\frac{\pi r}{\beta })\frac{\pi r}{\beta }\mathrm{sech}^2(\frac{\pi r}{\beta })\right)}{\beta r^3\left(1+\frac{\lambda ^2\pi }{r\beta }\mathrm{tanh}(\frac{\pi r}{\beta })\right)}}ϵ_{iak}x_k`$ (26)
For low temperatures,

$`\stackrel{~}{A}_i^a(\stackrel{}{x},t={\displaystyle \frac{\beta }{2}}){\displaystyle \frac{2}{3}}\lambda ^2\pi ^4T^4ϵ_{iak}x_k`$ (27)
which is nonzero, but which vanishes at zero temperature, as in (16). At $`t=+\frac{\beta }{2}`$, the expression for $`\stackrel{~}{A}_i^a(\stackrel{}{x},t=+\frac{\beta }{2})`$ is more cumbersome to write out, but it is completely specified by the ansatz form (3) with

$`\stackrel{~}{\varphi }_1(r,t=+{\displaystyle \frac{\beta }{2}})`$ $`=`$ $`\mathrm{sin}f(r,t=+{\displaystyle \frac{\beta }{2}})\varphi _2(r,t=+{\displaystyle \frac{\beta }{2}})`$ (28)
$`\stackrel{~}{\varphi }_2(r,t=+{\displaystyle \frac{\beta }{2}})`$ $`=`$ $`\mathrm{cos}f(r,t=+{\displaystyle \frac{\beta }{2}})\varphi _2(r,t=+{\displaystyle \frac{\beta }{2}})`$ (29)
$`\stackrel{~}{A}_1(r,t=+{\displaystyle \frac{\beta }{2}})`$ $`=`$ $`_rf(r,t=+{\displaystyle \frac{\beta }{2}})`$ (30)
and where $`f`$ is given by (25), and $`\varphi _2`$ is as in (7).

We now ask how $`\stackrel{~}{A}_i^a(\stackrel{}{x},t=+\frac{\beta }{2})`$ is related to $`\stackrel{~}{A}_i^a(\stackrel{}{x},t=\frac{\beta }{2})`$. Since the ansatz functions satisfy $`\stackrel{~}{A}_1(r,t=\frac{\beta }{2})=0`$, and $`\varphi _2(r,t=+\frac{\beta }{2})=\stackrel{~}{\varphi }_2(r,t=\frac{\beta }{2})`$, this shows that $`\stackrel{~}{A}_i^a(\stackrel{}{x},t=+\frac{\beta }{2})`$ and $`\stackrel{~}{A}_i^a(\stackrel{}{x},t=\frac{\beta }{2})`$ are related by a gauge transformation

$`\stackrel{~}{A}_i(\stackrel{}{x},t=+{\displaystyle \frac{\beta }{2}})=W^1\stackrel{~}{A}_i(\stackrel{}{x},t={\displaystyle \frac{\beta }{2}})W+W^1_iW`$ (31)
where the static large gauge transformation $`W(\stackrel{}{x})`$ is

$`W(\stackrel{}{x})`$ $`=`$ $`e^{\frac{i}{2}f(r,t=+\frac{\beta }{2})\widehat{x}\stackrel{}{\sigma }}`$ (32)
$`=`$ $`\mathrm{exp}\left[i\pi {\displaystyle \frac{(\overline{\lambda }^22\overline{r}^2)\mathrm{sinh}\overline{r}\overline{\lambda }^2\overline{r}\mathrm{cosh}\overline{r}}{\overline{r}\sqrt{\mathrm{sinh}\overline{r}(4\overline{\lambda }^2\overline{r}\mathrm{cosh}\overline{r}+(\overline{\lambda }^4+4\overline{r}^2)\mathrm{sinh}\overline{r})}}}(\widehat{x}\stackrel{}{\sigma })\right]`$ (33)
Using (21), the winding number of this gauge transformation is

$`N`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\left[f(r,t=+{\displaystyle \frac{\beta }{2}})\mathrm{sin}f(r,t=+{\displaystyle \frac{\beta }{2}})\right]_{r=0}^{r=\mathrm{}}`$ (34)
$`=`$ $`1`$ (35)
for all inverse temperatures $`\beta `$ and scales $`\lambda `$.

Notice that (31) shows that in the Weyl gauge the gauge fields $`\stackrel{~}{A}_i`$ are not strictly periodic. Rather, $`\stackrel{~}{A}_i`$ is periodic up to a large gauge transformation. This fixed-time gauge transformation is $`W(\stackrel{}{x})=U(\stackrel{}{x},t=+\frac{\beta }{2})`$, where $`U(\stackrel{}{x},t)`$ is the time-dependent gauge transformation used to transform to the Weyl gauge. So, just as in the $`T=0`$ case, the initial and final gauge fields, $`\stackrel{~}{A}_i(\stackrel{}{x},t=\frac{\beta }{2})`$, are related by a static large gauge transformation of unit winding number. But at finite T, these initial and final gauge fields $`\stackrel{~}{A}_i(\stackrel{}{x},t=\frac{\beta }{2})`$ are not pure gauge. This means that the tunneling at finite T is not between minima of the classical Yang-Mills potential (20), but between fields of higher potential energy. With the ansatz (3), the corresponding magnetic field strength is (and since the instanton solution is self-dual, this also gives the electric field strength)

$`B_i^a`$ $`=`$ $`{\displaystyle \frac{\left(_r\varphi _1A_1\varphi _2\right)}{r^2}}ϵ_{iak}x_k+{\displaystyle \frac{\left(_r\varphi _2+A_1\varphi _1\right)}{r}}\left(\delta _{ia}{\displaystyle \frac{x_ix_a}{r^2}}\right)+{\displaystyle \frac{\left(1\varphi _1^2\varphi _2^2\right)}{r^4}}x_ix_a`$ (36)
It is straightforward to verify that $`\stackrel{~}{B}_i^a(\stackrel{}{x},t)`$ is nonvanishing at $`t=\frac{\beta }{2}`$. The magnetic potential energy (20) can be expressed in terms of the ansatz functions in (3) and (7) \[since $`V`$ is gauge invariant, it doesn’t matter whether we use the original gauge or the Weyl gauge\] :

$`V=2\pi {\displaystyle _0^{\mathrm{}}}𝑑r\left\{2(_r\varphi _1A_1\varphi _2)^2+2(_r\varphi _2+A_1\varphi _1)^2+{\displaystyle \frac{1}{r^2}}(1\varphi _1^2\varphi _2^2)^2\right\}`$ (37)
At zero temperature, $`V(t=\mathrm{})=0`$, since the zero temperature instanton is a pure gauge at $`t=\mathrm{}`$. But at finite temperature, $`V(t=\frac{\beta }{2})`$ is nonzero, even though it approaches zero as $`\beta \mathrm{}`$. It is straightforward to evaluate numerically $`V`$ in (37) for the zero and finite temperature instantons discussed above. In Figure 1 we plot the potential energy as a function of Euclidean time for the zero temperature instanton for which $`\rho (r,t)`$ is given by (8), and in Figure 2 for the finite temperature caloron for which $`\rho (r,t)`$ is given by (23).

We now consider the instanton’s topological charge and its relation with Chern-Simons numbers and winding numbers. Recall that the instanton charge may be expressed as

$`Q`$ $`=`$ $`{\displaystyle \frac{1}{32\pi ^2}}{\displaystyle d^4xϵ_{\mu \nu \rho \sigma }\mathrm{tr}\left(F_{\mu \nu }F_{\rho \sigma }\right)}`$ (38)
$``$ $`{\displaystyle d^4x_\mu K_\mu },\mathrm{where}K_\mu ={\displaystyle \frac{1}{8\pi ^2}}ϵ_{\mu \nu \rho \sigma }\mathrm{tr}\left(A_\nu _\rho A_\sigma +{\displaystyle \frac{2}{3}}A_\nu A_\rho A_\sigma \right)`$ (39)
This leads to a natural separation of Q as

$`Q={\displaystyle 𝑑t\frac{d}{dt}CS(t)}+{\displaystyle _0^{\mathrm{}}}𝑑r{\displaystyle \frac{d}{dr}}K(r)`$ (40)
In terms of the ansatz fields in (3), the Chern-Simons term, $`CSd^3xK_0`$, is

$`CS(t)={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑r\left[\varphi _1\varphi _2^{}\varphi _2\varphi _1^{}A_1(1\varphi _1^2\varphi _2^2)+\varphi _1^{}\right]`$ (41)
while

$`K(r)={\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑t\left[\varphi _1\dot{\varphi _2}+\varphi _2\dot{\varphi _1}+A_0(1\varphi _1^2\varphi _2^2)\dot{\varphi _1}\right]}`$ (42)
Here $`\varphi ^{}`$ and $`\dot{\varphi }`$ denote derivative with respect to $`r`$ and $`t`$ respectively. At zero temperature, the $`t`$ integration in (42) is from $`\mathrm{}`$ to $`+\mathrm{}`$, while at finite temperature it is over $`[\frac{\beta }{2},+\frac{\beta }{2}]`$. Note that both $`CS(t)`$ and $`K(r)`$ are gauge dependent, while $`Q`$ itself is gauge invariant. Also note that the last terms in each of $`CS(t)`$ and $`K(r)`$ cancel one another in the decomposition (40), but are required for the correct transformation properties of the respective terms. For example, under a gauge transformation of the form (9), the Chern-Simons term changes as

$`\stackrel{~}{CS}`$ $`=`$ $`CS+{\displaystyle \frac{1}{8\pi ^2}}{\displaystyle d^3xϵ_{ijk}_i\mathrm{tr}\left(_jUU^1A_k\right)}+{\displaystyle \frac{1}{24\pi ^2}}{\displaystyle d^3xϵ_{ijk}\mathrm{tr}\left(U^1_iUU^1_jUU^1_kU\right)}`$ (43)
$`=`$ $`CS+{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑r{\displaystyle \frac{d}{dr}}\left[(\varphi _21)\mathrm{sin}f+\varphi _1(\mathrm{cos}f1)\right]{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑r{\displaystyle \frac{d}{dr}}\left[f\mathrm{sin}f\right]`$ (44)
It is a simple but instructive exercise to check that this nonabelian transformation law agrees with the transformation of (41) under the abelian-like transformations (13) of the ansatz functions, provided one includes the $`\varphi _1^{}`$ term in $`CS(t)`$ .

For the instanton solutions studied in this paper, one can evaluate numerically the functions $`CS(t)`$ and $`K(r)`$ in (41) and (42), in the original gauge (7), and in the Weyl gauge $`\stackrel{~}{A}_i`$ where $`\stackrel{~}{A}_0=0`$. At zero temperature, with the scalar function $`\rho (r,t)`$ given by (8), the functions $`CS(t)`$ and $`K(r)`$ are plotted in Figures 3 and 4. In this gauge, $`CS(t=+\mathrm{})CS(t=\mathrm{})=0`$, and $`K(r=\mathrm{})K(r=0)=1`$. Thus, in this gauge the entire contribution to the instanton charge $`Q`$ comes from the $`K(r)`$ term in the decomposition (40). On the other hand, in the Weyl gauge, the functions $`\stackrel{~}{CS}(t)`$ and $`\stackrel{~}{K}(r)`$ are plotted in Figures 5 and 6. In the Weyl gauge, the entire contribution to $`Q`$ comes from the change in the Chern-Simons number. Also note that in the Weyl gauge the Chern-Simons number is an integer at $`t=\pm \mathrm{}`$; this is because at $`t=\pm \mathrm{}`$ the Weyl gauge field is pure gauge \[see (16) and (18)\], and for a pure gauge field the Chern-Simons term is equal to the winding number of the corresponding group element. The time-dependent gauge transformation $`U=\mathrm{exp}(\frac{i}{2}f(r,t)\widehat{x}\stackrel{}{\sigma })`$, with $`f(r,t)`$ in (15), that transforms from the original singular gauge to the Weyl gauge can be considered to have a t-dependent winding number $`N(t)`$ that, from (21) and (15), is a step function in time:

$`N(t)={\displaystyle \frac{1}{2}}\left(1+\mathrm{sign}(t)\right)`$ (45)
This is consistent with the plots in Figures 3 and 5, since the middle term in (44) vanishes in this case. Thus, in the Weyl gauge, the Chern-Simons number provides a meaningful label for the distinct classical wells of the Yang-Mills potential. At $`t=0`$, the Weyl gauge Chern-Simons number is $`\frac{1}{2}`$, corresponding to the sphaleron at the peak of the barrier between neighboring wells. The winding number of the gauge transformation $`U(\stackrel{}{x},t)`$ is also $`\frac{1}{2}`$ when $`t=0`$.

At finite temperature, the situation is similar. The only change is that the scalar function $`\rho (r,t)`$ is now given by (23), and the t integration in (42) is over the finite interval $`[\frac{\beta }{2},+\frac{\beta }{2}]`$. In the Harrington-Shepard gauge (7,23), the functions $`CS(t)`$ and $`K(r)`$ are plotted in Figures 7 and 8. (These plots, and those in Figs. 9 and 10, are for scale parameter $`\overline{\lambda }=1`$. For different $`\overline{\lambda }`$, the plots have the same form.) Once again, in this gauge, the entire contribution to the instanton charge comes from the change in $`K(r)`$. The corresponding plots for the Weyl gauge are in Figures 9 and 10, and we see that the entire contribution to $`Q`$ comes from the change $`CS(t)`$. From the fact (31) that in the Weyl gauge $`\stackrel{~}{A}_i(t=+\frac{\beta }{2})`$ is related to $`\stackrel{~}{A}_i(t=\frac{\beta }{2})`$ by a static gauge transformation of winding number 1, we know that the difference $`\stackrel{~}{CS}(t=+\frac{\beta }{2})\stackrel{~}{CS}(t=\frac{\beta }{2})=1`$. But it is interesting to note further that in the Weyl gauge the Chern-Simons number is an integer at $`t=\pm \frac{\beta }{2}`$, even though the gauge field is not a pure gauge \[see (26) and (31)\]. At $`t=0`$, the Weyl gauge Chern-Simons number is $`\frac{1}{2}`$, as is the winding number of the gauge transformation $`U(\stackrel{}{x},t)`$ that connects the HS gauge to the Weyl gauge.

In conclusion, we have demonstrated explicitly that while the HS caloron solution is manifestly periodic, so that $`A_\mu (\stackrel{}{x},t=+\frac{\beta }{2})=A_\mu (\stackrel{}{x},t=\frac{\beta }{2})`$, when transformed to the Weyl ($`A_0=0`$) gauge, the gauge field is periodic up to a large gauge transformation, $`\stackrel{~}{A}_i(\stackrel{}{x},t=+\frac{\beta }{2})=W^1\stackrel{~}{A}_i(\stackrel{}{x},t=\frac{\beta }{2})W+W^1_iW`$, where $`W`$ is a fixed-time gauge transformation with integer winding number equal to the caloron’s topological charge. This gauge transformation $`W`$ connects different Chern-Simons-number sectors of the classical Yang-Mills potential energy, just as is familiar from zero temperature. However, the caloron solution in Weyl gauge interpolates, as $`t`$ goes from $`\frac{\beta }{2}`$ to $`+\frac{\beta }{2}`$, not between classical minima but between classical gauge fields of nonzero potential energy.

Having understood the tunneling interpretation of the HS calorons in Weyl gauge, we conclude with some brief comments concerning the difference between these calorons and the ”periodic instantons” studied in . These are very different objects, but the names can be somewhat confusing since calorons are themselves instantons that are periodic in Euclidean time. This distinction, and some of its physical consequences, have been stated clearly in , but we re-iterate these distinctions here in the light of our Weyl gauge calorons. The ”periodic instantons” are most easily introduced in the simple context of a quantum mechanical model; for example, a double well potential $`V=(q^21)^2/4`$. This captures most of the essential physics. This potential model has well-known instanton solutions  that are periodic in Euclidean time t:

$`q(t)=\sqrt{{\displaystyle \frac{2\nu }{1+\nu }}}\mathrm{sn}({\displaystyle \frac{t}{\sqrt{1+\nu }}}K;\nu ),`$ (46)
where $`\mathrm{sn}(x;\nu )`$ is a Jacobi elliptic function , and the modulus parameter $`0\nu 1`$ is related to the energy by $`E=(\nu 1)^2/(4(\nu +1)^2)`$. These solutions are manifestly periodic: $`q(t=\frac{\beta }{2})=q(t=+\frac{\beta }{2})`$, where the period $`\beta =4\sqrt{1+\nu }K(\nu )`$, and $`K(\nu )`$ is the elliptic quarter-period . Physically, these periodic instantons correspond to tunneling from one well to the next, and then back again to the original starting position (since they are periodic!). By contrast, the Weyl gauge calorons describe tunneling from one classical well (characterized by Chern-Simons number 0) to the next well (with CS=1); they do not tunnel back again, since the Weyl gauge field $`\stackrel{~}{A}_i(t=+\frac{\beta }{2})`$ differs from $`\stackrel{~}{A}_i(t=\frac{\beta }{2})`$ by a large gauge transformation, as we have demonstrated explicitly in (31,33). This difference between calorons and periodic instantons also manifests itself through the fact that the calorons have nonzero topological charge per period (the explicit ones studied here have Q=1), while the ”periodic instantons” have zero net topological charge per period. For low temperatures, the periodic instantons can be approximated by an infinite chain of alternating zero temperature instantons and anti-instantons. Indeed, this is particularly easy to see in the quantum mechanical double-well potential case (46), owing to the remarkable elliptic function identity:

$`\mathrm{sn}({\displaystyle \frac{t}{\sqrt{1+\nu }}};\nu )={\displaystyle \frac{\pi }{2\sqrt{\nu }K^{}}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}(1)^n\mathrm{tanh}\left({\displaystyle \frac{\pi }{2\sqrt{1+\nu }K^{}}}(tn{\displaystyle \frac{\beta }{2}})\right)`$ (47)
Here $`K^{}(\nu )K(1\nu )`$. This leads to the simple interpretation of the periodic instanton in (46) as a series of alternating instantons and anti-instantons. The zero temperature limit corresponds to $`\nu 1`$, in which case the period $`\beta \mathrm{}`$, and $`E0`$, and we approach the usual zero temperature kink-like instantons. But, on the full period $`[\frac{\beta }{2},\frac{\beta }{2}]`$, the $`sn`$ function reduces to a pair of an anti-kink and a kink. This explains why the net topological charge per period is zero. It is only on the half-period $`[0,\frac{\beta }{2}]`$ that $`q(t)`$ approaches the familiar zero temperature instanton: $`q=\mathrm{tanh}(t/\sqrt{2})`$. By contrast, the HS caloron is a periodic sum of instantons (and an HS anti-caloron is a periodic sum of anti-instantons), not a periodic sum of alternating instantons and anti-instantons. Indeed, in Yang-Mills theory there is no known exact ”periodic instanton” that is an exact alternating chain of instantons and anti-instantons, although approximate chains of this form have been extensively studied . This difference between calorons and periodic instantons can be reconciled because even though the caloron solutions are manifestly periodic, $`A_\mu (t+\beta )=A_\mu (t)`$, since this is a gauge theory, the $`A_\mu `$ are not all independent physical fields; when we transform the caloron to the Weyl gauge we see that the coordinate fields $`\stackrel{~}{A}_i`$ are periodic only up to a large gauge tranformation. We hope that our explicit construction of the Weyl gauge caloron fields helps to clarify this important distinction between calorons and ”periodic instantons”.

Acknowledgments:

GD is supported by the U.S. DOE grant DE-FG02-92ER40716.00. GD also thanks Balliol College, Oxford, the Theoretical Physics Department at Oxford, and the CSSM at Adelaide University, where this work was begun. BT is supported by PPARC.


text:
# Quasi-local Update Algorithms for Numerical Simulations of d=3 SU(2) LGT in the Dual Formulation

## 1 THE DUAL FORMULATION

A brief introduction to the dual formalism for $`d=3`$ $`SU(2)`$ lattice gauge theory as well as to the techniques that have been developed for its numerical simulations was presented at LATTICE’99 . As was stated there, the partition function of the conventional LGT can be converted, upon using the techniques of character expansion and group integrations, into the dual form

$$Z_d=\underset{\{j\}}{}(2j+1)C_{j_a}(\beta )\underset{i=1}{\overset{5}{}}\left\{\begin{array}{ccc}a_i& b_i& c_i\\ d_i& e_i& f_i\end{array}\right\}$$
(1)
This partition sum is defined over the dual lattice where $`\{j_a\}`$ live over the links and $`\{j_b\}`$ over the diagonals to the plaquettes. The convention for the diagonals is that they connect the vertices of the odd sublattice. We have collectively designated $`\{j_a\},\{j_b\}`$ by $`\{j\}`$. Each cube of the dual lattice is spanned by 5 tetrahedra of which one is spanned entirely by $`\{j_b\}`$ while four are spanned by three $`\{j_a\}`$ and three $`\{j_b\}`$. In what follows, the b-links will be shown by dashed lines while the a-links will be shown by solid lines. Each tetrahedron carries a weight factor which is the $`SU(2)`$ 6-j symbol $`\left\{\begin{array}{ccc}a& b& c\\ d& e& f\end{array}\right\}`$. Periodic b.c for the original lattice is crucial for this construction.

## 2 QUASI-LOCAL UPDATES

An important ingredient in the techniques for numerical simulation of such classes of models is an updating technique that respects all the triangle inequalities. In  we introduced two distinct classes of updating algorithms: local and quasi-local. The local algorithms work very much like the usual updating techniques for Monte Carlo simulations: all links but one are held fixed and the free link is updated by either Metropolis or heat bath methods. In the present context, one can see immediately that such local updates are not ergodic. To see this, let us recall that the representations of $`SU(2)`$ can be grouped into two classes i.e half-integral which can be called fermionic and integral which can be called bosonic. Let us call a move that takes fermionic links to bosonic, and vice versa, a $`Z_2`$-flip. It is obvious that in every triangle exactly two links must be $`Z_2`$-flipped simulateneously if we are to maintain the triangle inequalities. As the local moves change only one link at a time, they can not accommadate $`Z_2`$-flips and the moves therefore never take configurations out of their $`Z_2`$-classes.

To circumvent this the quasi-local moves were introduced. Since in each triangle exactly two $`Z_2`$-flips have to be carried out, these flips can proliferate throughout the lattice unless a way out is found. Therefore one can look for those clusters of $`Z_2`$-flips that have the smallest volume possible. In this way we identified the globeven moves at even-sublattice points which flip the six interior a-links as shown in fig.1. All the ”exterior” b-links are held fixed during a globeven move and this contains the $`Z_2`$-flips to the interior. It should be noted that eight tetrahedra get updated during each globeven move.

While the globeven moves are capable of $`Z_2`$-flipping the a-links, they leave the b-links unchanged. For this reason, these moves are also not ergodic. Since b-links connect only odd sites, any quasi-local move affecting them must take place at the odd sites. It then follows that the smallest volume of $`Z_2`$-flippings involves all the 6 a-links and 12 b-links at the odd site while keeping all the other a-links and b-links fixed with 32 tetrahedra getting updated. This is illustrated in fig.2.

## 3 KAGOME VARIABLES

Though we have identified the geometry of the cluster of links to be updated during quasi-local moves, we have to address the issue of how the new spin values are to be chosen such that all triangle inequalities are respected. Let us first recall the form of triangle inequalities. Considering the triangle $`(j_1,j_2,j_{12})`$ of fig.3, these inequalities take the form

$`|j_1j_2|`$ $`j_{12}`$ $`j_1+j_2`$
$`|j_1j_{12}|`$ $`j_2`$ $`j_1+j_{12}`$
$`|j_2j_{12}|`$ $`j_1`$ $`j_2+j_{12}`$ (2)
We illustrate the K-variable strategy with the example of a single tetrahedra. The idea behind these variables is to replace inequalities by equalities to some extent. They were first introduced by Bargman  and were later developed by Anishetty et al . Again, let us first introduce them for the triangle considered above. Their construction involves a pair $`(n_1,n_2)`$ for each link or $`j`$-value such that

$`n_1+n_2`$ $`=j_1`$
$`n_1+n_3`$ $`=j_2`$
$`n_3+n_2`$ $`=j_{12}`$ (3)
These relations are invertible:

$`2n_1`$ $`=j_1+j_2j_{12}`$
$`2n_2`$ $`=j_{12}+j_1j_2`$
$`2n_3`$ $`=j_2+j_{12}j_1`$ (4)
From eqn (4) it is clear that $`n_i`$ are either integers or half-integers and also positive semi-definite. The moral of the story is that we can choose any triplet $`(n_1,n_2,n_3)`$ and the $`j_i`$ constructed from them according to eqn (3) will automatically satisfy the triangle inequalities of eqn (2).

Now let us consider the tetrahedron of fig.3 and ask how we can choose $`(j_1,j_2,j_3,j_{12},j_{23},j_{13})`$ such that the four sets of triangle inequalities are all satisfied. At the level of $`j_i`$ this can be done as follows: freely choose 3 $`j`$’s originating at a vertex e.g $`(j_1,j_2,j_3)`$. Next choose $`j_{12}`$ in the range given in eqn (2) and choose $`j_{23}`$ likewise. Finally choose $`j_{13}`$ such that it simultaneously satisfies the triangle inequalities relevant for the two triangles $`(j_1,j_3,j_{13})`$ and $`(j_{12},j_{23},j_{13})`$. This procedure gets increasingly tedious. Not only is this procedure tedious, it does not even guarantee a solution always even though in the case of the single tetrahedron a solution is always possible.

That a solution using this strategy may not always be possible can be seen by applying it to the gobeven geometry of fig.1. Here one could start by first selecting a $`Z_2`$-flipped value for, say, OE. The link OA could then be chosen to satisfy the inequalities for OAE, the link OB chosen to satisfy OBE and OAB. So far the situation is like that for a single tetrahedron and all would be fine. Next, OC will have to be chosen to satisfy OCE and OBC; already at this stage there is no guarantee that these two sets can be simultaneously realised. OD will have to be chosen to satisfy ODE, AOD and COD; OF will have to satisfy FOA, FOB, FOC and FOD! The situation only gets worse with these additional requirements.

Now we show how the same problem can be efficiently addressed by using the K-variables.Again we illustrate the method by first applying it to the simple case of a single tetrahedron. For this purpose we introduce the variables $`(n_1,n_2,n_3)`$ for the set $`(j_1,j_2,j_{12})`$, $`(n_4,n_5,n_6)`$ for $`(j_2,j_3,j_{23})`$, $`(n_7,n_8,n_9)`$ for $`(j_1,j_3,j_{13})`$, and $`(l_1,l_2,l_3)`$ for $`(j_{12},j_{23},j_{13})`$. Writing down the various equalities we get

$`j_1`$ $`=n_2+n_3`$ $`=n_8+n_9`$
$`j_2`$ $`=n_1+n_3`$ $`=n_5+n_6`$
$`j_3`$ $`=n_4+n_6`$ $`=n_7+n_9`$ (5)
and

$`j_{12}`$ $`=n_1+n_2`$ $`=l_2+l_3`$
$`j_{23}`$ $`=n_4+n_5`$ $`=l_1+l_3`$
$`j_{13}`$ $`=n_7+n_8`$ $`=l_2+l_1`$ (6)
From eqn(5) it follows that $`(n_1,n_2,n_3,n_4)`$ can be chosen to be unrestricted, one of $`(n_5,n_6)`$ can be chosen to be independent but with restricted range and one of $`(n_7,n_8,n_9)`$ to be also independent but with restricted range. Choosing $`(n_6,n_9)`$ to be the independent ones, the relevant restricted ranges are

$`0`$ $`n_6`$ $`n_1+n_3`$
$`0`$ $`n_9`$ $`min(n_2+n_3,n_4+n_6)`$ (7)
Though $`(l_1,l_2,l_3)`$ are determined by the $`n_i`$, the fact they all have to be positive semidefinite yields additional restrictions which are best seen by solving for $`l_i`$ using eqn (6):

$`l_1`$ $`=n_3+n_4n_9`$
$`l_2`$ $`=n_2+n_6n_9`$
$`l_3`$ $`=n_1+n_9n_6`$ (8)
Summarising, the solution to the single tetrahedron problem in terms of the K-variables is: $`(n_1,n_2,n_3,n_4)`$ are unrestricted except for positive semi-definiteness, $`(n_6,n_9)`$ are independent but restricted in range and $`(n_5,n_7,n_8)`$ are dependent:

$`0`$ $`n_6`$ $`n_1+n_3`$
$`n_6n_1`$ $`n_9`$ $`min(n_2,n_4)+min(n_3,n_6)`$ (9)
and

$`n_5`$ $`=n_1+n_3n_6`$
$`n_7`$ $`=n_4+n_6n_9`$
$`n_8`$ $`=n_2+n_3n_9`$ (10)
### 3.1 TWO TETRAHEDRA CELL

Our solution to the single tetrahedron given by eqn(9) will form the basis for the eventual solution of the globeven and globodd problems. The final solution will be arrived at in three stages: i) the solution of the two-tetrahedra problem, ii) the solution of the four-tetrahedra problem and iii) globeven and globodd using the results of ii). We now show how the two-tetrahedra problem is solved. For this purpose, consider two tetrahedra glued together as shown in fig.4. We take the common triangle ABC with $`AB=j_1,BC=j_2,CA=j_{12}`$ to be associated with the independent K-variables $`(n_1,n_2,n_3)`$. For the second tetrahedron ABCD we have to introduce 6 more variables but there will be 3 more relations associated with the links AB,BC,CA thus reducing the additional variables to 3.

The spin assignments are taken to be : $`BE=j_3,AE=j_{13},CE=j_{23}`$ and $`BD=k_3,AD=k_{13},CD=k_{23}`$ and the additional K-variable assignments to be : $`(n_4,n_5,n_6)`$ with $`(j_3,j_2,j_{23})`$, $`(n_7,n_8,n_9)`$ with $`(j_1,j_3,j_{13})`$, $`(m_4,m_5,m_6)`$ with $`(k_3,j_2,k_{23})`$, and $`(m_7,m_8,m_9)`$ with $`(j_1,k_3,k_{13})`$. We give only the final results: the unrestricted variables are $`(n_1,n_2,n_3,n_5,m_5)`$ and the restricted but independent variables are $`(n_6,n_9,m_6,m_9)`$ and the remaining dependent variables are $`(n_4,n_7,n_8,m_4,m_7,m_8)`$. The equation corresponding to eqn (9) are:

$`0`$ $`n_6`$ $`n_1+n_3`$
$`n_6n_1`$ $`n_9`$ $`min(n_2,n_5)+min(n_3,n_6)`$
$`0`$ $`m_6`$ $`n_1+n_3`$
$`m_6n_1`$ $`m_9`$ $`min(n_2,m_5)+min(n_3,m_6)`$
and

$`n_4`$ $`=n_1+n_3n_6`$
$`n_7`$ $`=n_5+n_6n_9`$
$`n_8`$ $`=n_2+n_3n_9`$
$`m_4`$ $`=n_1+n_3m_6`$
$`m_7`$ $`=m_5+m_6m_9`$
$`m_8`$ $`=n_2+n_3m_9`$ (12)
### 3.2 FOUR TETRAHEDRA CELL

As the final ingredient in the K-variable solution to the quasi-local updates, we now consider four tetrahedra glued together as shown in fig.5. In the two tetrahedra case we finally had 15 K-variables of which 9 were independent and 6 dependent. The K-variable assignments are summarised below:

$`(OA,AE,OE)`$ $`(n1,n2,n3)`$
$`(DA,AE,DE)`$ $`(n4,n5,n6)`$
$`(OA,AD,OD)`$ $`(n7,n8,n9)`$
$`(OA,AC,OC)`$ $`(p1,p2,p3)`$
$`(DA,AC,DC)`$ $`(p4,p5,p6)`$
$`(BA,AC,BC)`$ $`(q4,q5,q6)`$
$`(BA,AE,BE)`$ $`(m4,m5,m6)`$
$`(OA,AB,OB)`$ $`(m7,m8,m9)`$ (13)
Of the 24 K-variables 6 are unrestricted, 7 are restricted but independent and 11 are dependent. These are summarised in the following:

$`n_1+n_3`$ $`n_6`$ $`0`$
$`n_1+n_3`$ $`m_6`$ $`0`$
$`n_2+n_3`$ $`p_2`$ $`0`$ (14)
$`min(n_2,n_5)+min(n_3,n_6)`$ $`n_9`$ $`n_6n_1`$
$`min(n_2,m_5)+min(n_3,m_6)`$ $`m_9`$ $`m_6n_1`$
$`min(n_5+n_6,n_2+n_3+p_1p_2)`$ $`p_6`$ $`0`$
$`min(m_5+m_6,n_2+n_3+p_1p_2)`$ $`q_6`$ $`0`$
The new feature at this level is that range restrictions are now of two types: those that depend only on the unrestricted set as in the case of $`(n_6,m_6,p_2)`$ and those that depend on the already restricted set as in the case of $`(n_9,m_9,p_6,q_6)`$. The dependent variables are given by:

$`n_4`$ $`=n_1+n_3n_6`$
$`n_7`$ $`=n_5+n_6n_9`$
$`n_8`$ $`=n_2+n_3n_9`$
$`m_4`$ $`=n_1+n_3m_6`$
$`m_7`$ $`=m_5+m_6m_9`$
$`m_8`$ $`=n_2+n_3m_9`$
$`p_3`$ $`=n_2+n_3p_2`$
$`p_4`$ $`=p_1+p_3p_6`$
$`p_5`$ $`=n_5+n_6p_6`$
$`q_4`$ $`=p_1+p_3q_6`$
$`q_5`$ $`=m_5+m_6q_6`$ (16)
## 4 GLOBEVEN AND GLOBODD

Finally we are in a position to solve the globeven and globodd problems. For globeven, one takes two copies of the 4-tetrahedra cell and glues them such that the rectangular bases of the two cells coincide. As before the number of K-variables required will be twice the number for each 4-tetrahedra cell minus the number associated with the rectangular base. But because of the 8 links that are shared between the two cells, the independent variables will be reduced by 8. The results are similar to eqns(15,16) but more lengthy and will not be displayed here.

For the globodd case we have to take six copies of the 4-tetrahedra cells and glue them together such that the six rectangular bases now form the skin of the 8-cube cluster shown in fig.2. Now the total number of independent K-variables will be six times those of the 4-tetrahedra cell reduced by the 12 new relations that arise due to the sharing of the 12 b-links.

This completes the solution to the problem of finding spin-configurations for the globoeven and globodd updates that will respect all the triangle inequalities. Yet, the resulting equations are quite involved and may not result in efficient algorithms for numerical simulations. A practical approach that can be followed to accomplish a $`Z_2`$-flip is to simply increase all the interior spins by $`\frac{1}{2}`$ except those that are already at the maximum value permitted by the code; for these, the spin value can be reduced by $`\frac{1}{2}`$. It is easily seen that this ”brute force” $`Z_2`$-flip will also maintain the triangle inequalities. But such moves do not necessarily respect the balance conditions required of Monte Carlo simulations. As a way out of this, we swept through the relevant clusters many times randomly setting the links to new values allowed by triangle inequalities. This was achieved through the local update methods mentioned before. New values were accepted according to a metropolis algorithm or a heat bath method carried over a small chosen subset. The random numbers were generated uniformly in the range $`(0,1)`$..Whether this method is picking new link samples in an unbiased way is yet to be understood satisfactorily .


text:
# The relationship between fragility, configurational entropy and the potential energy landscape of glass forming liquids

Glass is a microscopically disordered, solid form of matter that results when a fluid is cooled or compressed in such a fashion that it does not crystallise. Almost all types of materials are capable of glass formation – polymers, metal alloys, and molten salts, to name a few. Given such diversity, organising principles which systematise data concerning glass formation are invaluable. One such principle is the classification of glass formers according to their fragility. Fragility measures the rapidity with which a liquid’s properties such as viscosity change as the glassy state is approached. Although the relationship between features of the energy landscape of a glass former, its configurational entropy and fragility have been analysed previously (e. g.,), an understanding of the origins of fragility in these features is far from being well established. Results for a model liquid, whose fragility depends on its bulk density, are presented in this letter. Analysis of the relationship between fragility and quantitative measures of the energy landscape (the complicated dependence of energy on configuration) reveal that the fragility depends on changes in the vibrational properties of individual energy basins, in addition to the total number of such basins present, and their spread in energy. A thermodynamic expression for fragility is derived, which is in quantitative agreement with kinetic fragilities obtained from the liquid’s diffusivity.

Glass forming liquids grow increasingly viscous upon cooling, till the viscosity becomes so large that they fail to flow on experimental time scales. While remaining microscopically disordered like a liquid, they manifest mechanical properties of a solid; i. e., they transform to a glass. By convention, the glass transition temperature $`T_g`$ is that where the viscosity reaches a value of $`10^{12}`$ Pa sec. The approach to this large viscosity, however, differs from one liquid to another. When displayed in an Arrhenius plot \[log(viscosity) vs. inverse temperature $`1/T`$\], some liquids (e. g., silica) show a steady, linear increase, while others display a much steeper dependence on $`1/T`$, as illustrated in Fig. 1(a) (inset). The former are strong liquids, and the latter, fragile. This range of behaviours is implicit in the Vogel-Fulcher-Tammann-Hesse (VFT) form, observed to describe the $`T`$ dependence of viscosity (as well as diffusivity and relaxation times) in many glass formers. The VFT relation may be written as

$$\eta =\eta _0\mathrm{exp}\left[\frac{1}{K_{VFT}(T/T_01)}\right],$$
(1)
where $`T_0`$ is the temperature of apparent divergence of viscosity, and $`K_{VFT}`$ is a material specific parameter quantifying the kinetic fragility; more fragile liquids have larger $`K_{VFT}`$ values. This behaviour of viscosity correlates with the jump in heat capacity at the glass transition – more fragile the liquid, sharper and bigger the jump. Rationalisation of this correlation comes from the Adam-Gibbs relationship which predicts a dependence of viscosity on the configurational entropy $`S_c`$ (described later) of the liquid:

$$\eta =\eta _0\mathrm{exp}\left[\frac{A}{TS_c}\right],$$
(2)
which results in the VFT relation if $`S_c`$ has the form $`TS_c=K_{AG}(T/T_K1)`$, where $`T_K`$ is the ideal glass transition or Kauzmann temperature below which $`S_c=0`$. The jump in heat capacity is proportional to $`K_{AG}`$, and is larger for more fragile glass formers.

Fragility is analysed here for a binary mixture of particles interacting via the Lennard-Jones potential, widely studied as a model glass former. The $`T`$ dependence of diffusivity $`D`$ (to characterise dynamics) and configurational entropy are obtained for a range of bulk densities $`\rho `$. The diffusivities $`D(\rho ,T)`$, are shown in Fig. 1(a) in a ‘fragility plot’: $`\mathrm{log}[D(\rho ,T)]`$ vs. $`T_r/T`$, where $`T_r`$ is a reference temperature (akin to $`T_g`$ in the usual fragility plot) at which $`D`$ at each density reaches a fixed, small value. By comparison with the inset, the liquid is seen to become more fragile as its density increases. Fig. 1(b) shows the kinetic fragility index, $`K_{VFT}`$, obtained from VFT (Eq.(1)) fits to diffusivity data at each density. The density dependence of $`K_{VFT}`$ quantifies the increase of fragility with density, apparent from data in Fig. 1(a).

A thermodynamic explanation of fragility, as discussed above (also ), lies in explaining the rapidity of change with $`T`$ of the configurational entropy, $`S_c`$. In the inherent structure (IS) formalism employed here, the configuration space of a liquid is divided into basins of local potential energy minima (‘inherent structures’). $`S_c`$ at a given $`T`$ derives from the multiplicity of local potential energy minima sampled by the liquid at that $`T`$, and equals the difference of the total entropy and the vibrational entropy of typical basins sampled, $`S_{vib}`$:

$$S_c(\rho ,T)=S_{total}(\rho ,T)S_{vib}(\rho ,T).$$
(3)
For all the results presented here, the basin entropy $`S_{vib}`$ is calculated by approximating each basin as a harmonic well (valid at sufficiently low $`T`$). In this approximation, the entropy of a given basin, arising from vibrational motion within the basin (the suffix vib emphasizes this point) is given by

$$S_{vib}=k_B\underset{i=1}{\overset{3N}{}}1\mathrm{log}(\frac{h\nu _i}{k_BT}),$$
(4)
where $`\nu _i`$ are the vibrational frequencies of the given basin, $`k_B`$ is Boltzmann’s constant and $`h`$ is Plank’s constant. The basin free energy is $`F_{vib}=U_{vib}TS_{vib}`$, where $`U_{vib}=3Nk_BT`$ is the internal energy. Apart from the explicit temperature dependence of $`S_{vib}`$ above, an implicit dependence also exists because vibrational frequencies change from one basin to another, and, different basins are sampled at different temperatures. From the expression above, it is clear that the entropy difference between any two basins arises solely from the difference in their vibrational frequencies, and is independent of temperature.

The configurational entropy is evaluated in this letter using two methods, referred to as the thermodynamic integration (TI) method, and the potential energy landscape (PEL) method, respectively. In both cases, inherent structures are obtained by performing local potential energy minimisation for a subset of configurations sampled by a liquid at each $`\rho `$ and $`T`$. Vibrational frequencies are calculated for each inherent structure, and the basin entropy $`S_{vib}`$ is obtained from Eq. (4). In the TI method (described in detail elsewhere; see also caption of Fig. 2), the total entropy $`S_{total}`$ is obtained using pressure and internal energy values from simulations. Eq. (3) is used to obtain $`S_c`$. In contrast, the PEL method is based on constructing the distribution of inherent structures in order to calculate $`S_{total}`$ and $`S_c`$, as described later.

Fig. 2(a) shows the Adam-Gibbs plot, where $`\mathrm{log}D`$ is plotted against $`(TS_c)^1`$. The Adam-Gibbs relation (Eq. (2)) is seen to apply well at each density, as also observed previously for the hard sphere liquid and water. Fig. 2(b) shows (as lines) the $`T`$ dependence of $`TS_c`$ (from TI). All curves are nearly linear, but with slight increases of slope at lower $`T`$. The slopes of $`TS_c`$ (in the $`T`$ range where diffusivities are measured), $`K_{AG}^{TI}`$, are used as the quantitative index of thermodynamic fragility, and are plotted in Fig. 1(b). $`K_{AG}^{TI}`$ are nearly the same as the kinetic fragility index $`K_{VFT}`$, thereby validating the idea that the rate of increase of $`TS_c`$ does indeed determine the fragility.

However, this analysis provides no insight into how the fragility relates to features of the energy landscape. To this end, the distribution of energy minima, and properties of individual basins, are now examined, at fixed density $`\rho `$. The probability distribution $`P(\mathrm{\Phi },T)`$ that an inherent structure of energy $`\mathrm{\Phi }`$ is sampled at temperature $`T`$ depends on (a) the energy $`\mathrm{\Phi }`$, (b) the structure of the basin, contained in the basin free energy $`F_{vib}(\mathrm{\Phi },T)`$ and (c) the number density of inherent structures with energy $`\mathrm{\Phi }`$, $`\mathrm{\Omega }(\mathrm{\Phi })`$:

$$P(\mathrm{\Phi },T)=\mathrm{\Omega }(\mathrm{\Phi })exp\left[\beta (\mathrm{\Phi }+F_{vib}(\mathrm{\Phi },T))\right]/Q_N(\rho ,T),$$
(5)
where $`Q_N`$, the partition function, normalises $`P`$ and is given by $`Q_N=𝑑\mathrm{\Phi }P(\mathrm{\Phi },T)`$ over all possible $`\mathrm{\Phi }`$. The configurational entropy density is defined by $`𝒮_c(\mathrm{\Phi })k_B\mathrm{ln}\mathrm{\Omega }(\mathrm{\Phi })`$ and is related to the $`T`$ dependent configurational entropy by $`S_c(T)=𝑑\mathrm{\Phi }𝒮_c(\mathrm{\Phi })P(\mathrm{\Phi },T)`$. The quantities $`P(\mathrm{\Phi },T)`$, $`F_{vib}`$ and $`Q_N`$ are obtained from simulations, and are used to invert the relation in Eq. (5) to obtain the configurational entropy density $`𝒮_c(\mathrm{\Phi })`$. The same $`𝒮_c(\mathrm{\Phi })`$ estimates must result (though not for the same range of $`\mathrm{\Phi }`$) for different $`T`$, when the harmonic assumption is valid, which is indeed seen to be the case from data shown in Fig. 3(a). With increasing density, the maximum $`𝒮_c(\mathrm{\Phi })`$ value reached decreases, and the distribution becomes broader. It has been argued that the distribution $`\mathrm{\Omega }(\mathrm{\Phi })`$ should be Gaussian (equivalently, $`𝒮_c(\mathrm{\Phi })`$ an inverted parabola). Data in Fig. 3(a) are fully consistent with this expectation, and are fitted to the form,

$$\frac{𝒮_c(\mathrm{\Phi })}{Nk_B}=\alpha \frac{(\mathrm{\Phi }\mathrm{\Phi }_o)^2}{\sigma ^2}$$
(6)
where $`\alpha `$ is the height of the parabola and determines the total number of configurational states, i. e. energy minima (the total number is proportional to $`\mathrm{exp}(\alpha N)`$ ), $`\mathrm{\Phi }_0`$ and $`\sigma ^2`$ respectively define the mean and the variance of the distribution. The fits are displayed in Fig. 3(a) and the fit parameters listed in the figure caption. Speedy has considered a model of the thermodynamics of a glass former with a Gaussian distribution of configurational states, and the basin specific heat that does not depend on the internal parameter (here the IS energy). He predicts that the fragility of a liquid should be proportional to the logarithm of the number of states, i. e. on $`\alpha `$. From Fig. 3(a) it is clear that $`\alpha `$ values, and hence the number of states, decrease with increasing density $`\rho `$. But the fragilities (Fig. 1(b)) show the opposite trend. However, fragility has an additional dependence (implicit in ref. ) on the variance of the distribution $`\sigma `$ (see equation (6)).

Further, the variation of vibration frequencies, and hence the basin entropy, from basin to basin, also contributes to the fragility of the liquid. Fig. 3(b) displays the $`T`$ independent (see above) basin entropy difference at a given energy $`\mathrm{\Phi }`$ with respect to the value at $`\mathrm{\Phi }_0`$, $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })S_{vib}(\mathrm{\Phi },T)S_{vib}(\mathrm{\Phi }_0,T)`$. At all densities, $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })`$ displays a clear $`\mathrm{\Phi }`$ dependence, which is inconsistent with assumptions in . Although such basin dependence of vibrational frequencies and hence the basin entropy have been discussed previously, with reference to both experimental data and model calculations, its implications on a liquid’s fragility have not clearly been established. Fig. 3(b) shows that basins at higher energy (which are occupied at higher temperature) have lower basin entropies, or, higher frequencies (see Eq. (4)). This trend, contrary to what one expects for systems studied experimentally, arises because the constant density model system here cannot expand when $`T`$ is raised.

The dependence of $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })`$ on $`\mathrm{\Phi }`$ is to a very good extent linear. Hence, for each density, $`\mathrm{\Delta }S_{vib}`$ is fit to the form $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })=\delta S(\mathrm{\Phi }\mathrm{\Phi }_0)`$ with fit parameters listed in the caption of Fig. 3. The basin free energy, which can be written as $`F_{vib}(\mathrm{\Phi },T)=F_{vib}(\mathrm{\Phi }_0,T)T\delta S(\mathrm{\Phi }\mathrm{\Phi }_0)`$, is also linear in $`\mathrm{\Phi }`$. With a Gaussian form for $`\mathrm{\Omega }`$, the partition function $`Q_N=𝑑\mathrm{\Phi }P(\mathrm{\Phi },T)`$ ($`P`$ as in Eq. (5)) is easily evaluated, and results in the following predictions for the $`T`$ dependence of the average IS energy, and the configurational entropy:

$$<\mathrm{\Phi }>(T)=\mathrm{\Phi }_0^{eff}\frac{\sigma ^2}{2Nk_BT},$$
(7)
where $`\mathrm{\Phi }_0^{eff}=\mathrm{\Phi }_0+\frac{\sigma ^2\delta S}{2Nk_B}`$, and

$$TS_c(T)=K_{AG}^{PEL}(T)(T/T_K1);K_{AG}^{PEL}(T)=\left(\frac{\sigma \sqrt{\alpha }}{2}+\frac{\sigma ^2\delta S}{4Nk_B}\right)\left(1+\frac{T_K}{T}\right)\frac{\sigma ^2\delta S}{2Nk_B},$$
(8)
where $`T_K=\sigma (2Nk_B\sqrt{\alpha }+\sigma \delta S)^1`$ is the ideal glass transition temperature where $`S_c(T_K)=0`$.

The form above for $`S_c`$ results in the VFT relation, if $`K_{AG}^{PEL}(T)`$ is constant. This is not quite the case, but the $`1/T`$ term becomes rapidly irrelevant for $`T>T_K`$ and $`K_{AG}^{PEL}(T)`$ approaches a constant. Further, $`TS_c`$ values plotted (as points) in Fig. 2(b) are in extremely good agreement with the ones from thermodynamic integration, including deviations from linearity for small $`T/T_K`$. A check of the overall consistency of the model described here is made in Fig. 3(c), which shows a scaled plot of $`<\mathrm{\Phi }>(T)\mathrm{\Phi }_0^{eff}`$ vs. $`\frac{\sigma ^2}{2Nk_BT}`$. From Eq. (7), one expects $`\mathrm{\Phi }(T)`$ data for all densities to collapse onto a straight line of (negative) unit slope, with expected deviations at high $`T`$ due to basin anharmonicity. Data in Fig. 3(c) clearly meet this expectation, although with noticeable deviations for $`\rho =1.35\rho _0`$, as also in Fig. 2(b).

The expression for $`K_{AG}^{PEL}`$ in Eq. (8), the thermodynamic fragility index obtained from energy landscape parameters, shows that fragility depends not only on the multiplicity of states ($`\alpha `$) but also on their spread ($`\sigma `$), and how much the basin entropy changes from the lowest to the highest energy basins sampled by the liquid. The extent of this change is quantified by parameters $`\sigma `$ and $`\delta S`$ (listed in the caption of Fig. 3). Fig. 1(b) shows the average slopes of $`TS_c`$ (from Eq. (8) above, in the $`T`$ range where diffusivities are measured), $`K_{AG}^{PEL}`$, which are nearly the same as $`K_{AG}^{TI}`$, and in very good agreement with the kinetic fragility index $`K_{VFT}`$. This agreement establishes the quantitative relationship between fragility and the energy landscape of the liquid since $`K_{AG}^{PEL}`$ is obtained solely in terms of a quantitative description of the liquid’s energy landscape.

Acknowledgements: I thank C. A. Angell, G. P. Johari, K. J. Rao, F. Sciortino, R. Seshadri, R. J. Speedy and U. V. Waghmare for very useful discussions and/or comments on the manuscript.

Correspondence should be addressed to the authour (email: sastry@jncasr.ac.in).

Figure 1: Fragility Plot of diffusivities, kinetic and thermodynamic fragility indices. Results shown are obtained from molecular dynamics simulations of $`204`$ $`A`$ type and $`52`$ $`B`$ type particles interacting via the Lennard-Jones (LJ) potential. Argon units are used for $`A`$ type particles: LJ parameter $`ϵ_{AA}=k_B\times 119.8K`$, $`\sigma _{AA}=0.3405\times 10^9m`$, $`m_A=6.6337\times 10^{26}kg`$, and $`ϵ_{AB}/ϵ_{AA}=1.5`$, $`ϵ_{BB}/ϵ_{AA}=0.5`$, $`\sigma _{AB}/\sigma _{AA}=0.8`$, and $`\sigma _{BB}/\sigma _{AA}=0.88`$, $`m_B/mA=1.`$ Densities are reported in units of $`\rho _0=2.53\times 10^{28}m^3`$ or $`1.678gcm^3`$. The liquid-gas critical point is located roughly at $`T_c=130K`$, $`\rho _c=0.416\rho _0`$. (a) Diffusivities $`D`$ from molecular dynamics simulations (details of the simulations are as in ) are displayed in a ‘fragility plot’ for $`5`$ densities. The reference temperature $`T_r`$ is chosen such that $`D(T_r)=2.44\times 10^8cm^2/s`$, which is slightly below the lowest $`D`$ values measured. The VFT extrapolation is used to locate $`T_r`$. Fit values $`D_0`$ (preexponent) are in the range of $`0.87`$ to $`2.0\times 10^4cm^2/s`$, close to the experimentally observed values. (b) Fragility index $`K`$ obtained (i) from VFT fits to diffusivities $`D`$, and from the configurational entropy obtained (ii) from thermodynamic integration (‘TI’), and (iii) from analysis of the potential energy landscape (‘PEL’). Dimensionless thermodynamic fragility indices are obtained by dividing $`K_{AG}^{TI}`$, $`K_{AG}^{PEL}`$ by $`Nϵ_{AA}`$. Comparison with data in  shows that at the low density end, the fragility of the model liquid compares with that of considerably fragile liquids such as toluene, orthoterphenyl and salol, while at the high density end, the fragility is extremely high. Temperatures of vanishing diffusivity $`T_0`$ from VFT fits are $`18.76K`$, $`24.98K`$, $`35.97K`$, $`47.86K`$ and $`75.75K`$ for $`\rho /\rho _0=1.1`$, $`1.15`$, $`1.2`$,$`1.25`$ and $`1.35`$ respectively.

Figure 2: Adam-Gibbs plot, configurational entropy vs. temperature. (a) The ‘Adam-Gibbs’ plot displaying $`\mathrm{log}(D)`$ vs. $`(TS_c)^1`$. $`TS_c`$ is expressed in units of $`ϵ_{AA}`$. The Adam-Gibbs expectation of a linear dependence is well satisfied at each density. (b) Configurational entropy obtained from thermodynamic integration (‘TI’) (lines, not labeled), and from analysis of the potential energy landscape (points, labeled as ‘PEL’), plotted against temperature $`T`$ scaled to the Kauzmann temperature $`T_K`$ obtained from thermodynamic integration. In the TI method, the total entropy $`S_{total}`$ of the liquid is obtained from its total free energy $`F_{total}`$, from $`F_{total}=U_{total}TS_{total}`$, where $`U_{total}`$ is the internal energy. Thermodynamic identity $`P=\rho ^2\left(\frac{F_{total}/N}{\rho }\right)_T`$ where $`P`$ is the pressure, is used to integrate $`P`$ from simulations to obtain $`F_{total}`$ with respect to the known ideal gas free energy. The identity $`U_{total}=\left(\frac{F_{total}/T}{1/T}\right)_\rho `$ is used to integrate at fixed density $`U_{total}`$ from simulations to obtain $`F_{total}`$ at any desired $`T`$. Extrapolations to low $`T`$ beyond the range of simulations utilize the form for the potential energies $`E`$, $`E=E_0+E_1T^{E_2}`$, where $`E_23/5`$ (see ). Basin entropies $`S_{vib}`$ are obtained by averaging over typical inherent structures at each $`T`$, whose vibrational frequencies are obtained numerically by diagonalising the Hessian (matrix of second derivatives of the potential energy). The plot covers at each density the range of $`T`$ where simulations have been performed. Scaling temperatures $`T_K^{TI}`$, the Kauzmann or ideal glass transition temperature where $`S_c=0`$, are: $`20.44K,26.83K,34.35K,42.49K`$, and $`62.47K`$, for $`\rho /\rho _0=1.1,1.15,1.2,1.25`$ and $`1.35`$ respectively. $`T_K^{PEL}`$ obtained from the PEL method are, in the same order, $`20.24K,26.18K,34.74K,44.64K`$, and $`62.35K`$.

Figure 3: Distribution of energy minima, basin entropy, and temperature dependence of inherent structure energies. (a) The configurational entropy density $`𝒮_c(\mathrm{\Phi })`$ plotted against inherent structure energy $`\mathrm{\Phi }`$ minus the ‘Kauzmann’ energy $`\mathrm{\Phi }_{min}`$ below which $`𝒮_c(\mathrm{\Phi })=0`$. $`\mathrm{\Phi }_{min}`$ is the lower root of $`\frac{𝒮_c(\mathrm{\Phi })}{Nk_B}=\alpha \frac{(\mathrm{\Phi }\mathrm{\Phi }_o)^2}{\sigma ^2}=0.`$ For each density, data from $`9`$ temperatures (on average) are shown (smaller symbols) which collapse onto each other, validating the harmonic approximation for calculating basin entropies (see text). Also shown are the averages of these curves (larger symbols), for whose calculation a threshold of $`20`$ samples (out of a total of $`1000`$) is imposed in each bin. (b) The $`T`$-independent component of the basin entropy $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })`$ vs. $`\mathrm{\Phi }\mathrm{\Phi }_o`$, where $`\mathrm{\Phi }_o`$ is the value where $`𝒮_c(\mathrm{\Phi })`$ is maximum. $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })`$ curves are approximately linear. The $`T`$ independence of $`\mathrm{\Delta }S_{vib}(\mathrm{\Phi })`$ is an artefact of the classical mechanical treatment employed here, but a quantum mechanical calculation of the basin entropy does not alter the picture fundamentally (unpublished results). (c) Scaled plot of the average IS energies sampled in simulations, which are expected to collapse on to the straight line shown at low temperatures (right side of the graph), if the assumption of a Gaussian configurational entropy density $`𝒮_c(\mathrm{\Phi })`$ is valid (see text). The observed data collapse validates the Gaussian form used for $`𝒮_c(\mathrm{\Phi })`$. Deviations at high temperatures are due to the breakdown of harmonicity of basins. Fit parameters $`\alpha `$, $`\varphi _0/Nϵ_{AA}`$, $`\sigma /Nϵ_{AA}`$ and $`\delta Sϵ_{AA}/k_B`$ are, ($`\rho =1.1\rho _0`$): $`1.010,6.586,0.303,0.721,`$ ($`\rho =1.15\rho _0`$): $`0.957,6.684,0.374,0.652,`$ ($`\rho =1.2\rho _0`$): $`0.933,6.683,0.489,0.500,`$ ($`\rho =1.25\rho _0`$): $`0.927,6.564,0.642,0.314,`$ and ($`\rho =1.35\rho _0`$): $`0.849,6.088,0.873,0.188.`$

SASTRY FIGURE 1

SASTRY FIGURE 2

SASTRY FIGURE 3