# Bonding in Solids

Charbel Tannous

► To cite this version:

Charbel Tannous. Bonding in Solids. Master. Physique du Solide, Université de Brest, France. 2018, pp.31. hal-04191449

HAL Id: hal-04191449

https://hal.science/hal-04191449v1

Submitted on 30 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Bonding in Solids

C. Tannous

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

(Dated: August 28, 2023)

Crystal cohesion is built on several types of bonding varying from molecular to ionic, covalent... In this work, bonding is explained from a microscopic point of view in order to describe different types of solids ranging from rare gas solids, to ionic, covalent, metallic... Crystal thermodynamics is derived from lattice sums based on non-interacting and many-body points of view originating from pair (Molecular and Ionic solids) or triplet (angular dependent) interactions (Covalent and semiconducting solids). Lattice sums are performed with Ewald method numerically with Error functions and analytically with Jacobi-Theta functions for any dimension.

PACS numbers: 61.50.-f; 73.20.Hb; 31.15.-p; 07.05.Tp

Keywords: Structure of bulk crystals; Surfaces and interfaces; Computational techniques for electronic structure atoms and molecules; Computer modeling and simulation

## Contents

<table><tr><td>I. Introduction</td><td>1</td></tr><tr><td>A. Pair and Many-body potentials</td><td>3</td></tr><tr><td>B. Reduced Units for pair potentials</td><td>3</td></tr><tr><td>C. Properties of the pair potential</td><td>4</td></tr><tr><td>D. Pair potential types</td><td>4</td></tr><tr><td>E. Dipoles</td><td>6</td></tr><tr><td>1. Point dipoles</td><td>6</td></tr><tr><td>2. Induced dipoles</td><td>6</td></tr><tr><td>F. Many-body potentials for covalently bonded systems</td><td>8</td></tr><tr><td>G. Topological Solids</td><td>9</td></tr><tr><td>II. Interaction Potentials and Elasticity</td><td>9</td></tr><tr><td>A. Stretching deformation energy</td><td>10</td></tr><tr><td>B. Bending deformation energy</td><td>10</td></tr><tr><td>C. Pair potentials and elastic constants</td><td>11</td></tr><tr><td>III. Lattice sums</td><td>12</td></tr><tr><td>A. Ewald summation</td><td>14</td></tr><tr><td>B. Theta function summation</td><td>15</td></tr><tr><td>IV. Classical Thermodynamics</td><td>18</td></tr><tr><td>A. Thermodynamic potential functions</td><td>18</td></tr><tr><td>B. Non-ideal gas thermodynamics</td><td>22</td></tr><tr><td>V. Molecular solids</td><td>25</td></tr><tr><td>A. Van der Waals interaction</td><td>25</td></tr><tr><td>VI. Ionic solids</td><td>26</td></tr><tr><td>A. Madelung constant</td><td>26</td></tr><tr><td>VII. Covalent solids</td><td>27</td></tr><tr><td>A. Hybridization</td><td>27</td></tr><tr><td>VIII. Metals</td><td>30</td></tr><tr><td>References</td><td>31</td></tr></table>

## I. INTRODUCTION

Crystals are made from a variety of elements that belong from various parts of the periodic table implying a host of outer electron configurations such as inert gas, transition metals (d shell), lanthanides (f shell), semiconductors
(sp-hybridization)....

The outer electron configurations induce various interactions between neighboring atoms/molecules affecting the nature of bonding between them.

The interaction potential between atoms/molecules belong to different families such as pair potentials: hard spheres, Lennard-Jones, Morse..., angular dependent potentials in the case of covalent crystals

Relation between interatomic potentials and lattice properties: lattice energy, elastic constants, vacancy energy...

Empirically, bonding energy (BE) originating from interaction between atoms/molecules is related to at least three quantities (1):

* Melting temperature is $\propto BE$

* Vapor pressure is $\propto 1/BE$

* Packing fraction (PF) is $\propto 1/BE$.

For instance a Simple Cubic (CS) crystal has PF=0.524, a body-centered crystal (BCC has PF=0.680, face centered (FCC) and Hexagonal compact (HCP) have PF=0.740 whereas Diamond has PF=0.34. This indicates that crystals with weak BE like molecular crystals should be FCC or HCP whereas strong crystals with high BE like Diamond do not need a high PF.

The main source of experimental information about intermolecular interactions are:

* Direct measurement with Atomic Force Microscopy. This is particularly useful for thin films and nanostructures.

* Spectroscopic measurements (vibrational-rotational spectra, pre-dissociation, broadening of lines by pressure, etc.)

* Measurement of mechanical and thermodynamic crystal properties (elastic constants, phonon spectra, sublimation energy...)

* Measurement of thermodynamic properties of gas and liquid phases of the crystal (virial expansion coefficients, viscosity, heat transport coefficients...)

* Measurement of cleavage energies along different Miller crystal planes (2; 3). Some examples are given in Table 1.

* Scattering experiments in atomic-molecular beams, which, in some cases, allow the potentials to be determined directly from experimental data.

<table><caption>**Table 1:** *E*<sub>||</sub> is Young modulus parallel to cleavage plane whereas *E*<sub>⊥</sub> is Young modulus perpendicular to cleavage plane. CGS pressure unit dyne/cm<sup>2</sup> is equivalent to a Giga-Pascal (10<sup>9</sup> Pascals) in SI. Cleavage energy is inversely proportional (2; 3) to Young modulus. Note that Zn crystal is hexagonal, thus the four digit Miller indexing.</caption><thead><tr><th>Crystal</th><th>Cleavage plane</th><th><i>E</i><sub>||</sub> (10<sup>11</sup> dyne/cm<sup>2</sup>)</th><th><i>E</i><sub>⊥</sub> (10<sup>11</sup> dyne/cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>LiF</td><td>(100)</td><td>10.2</td><td>10.2</td></tr><tr><td>MgO</td><td>(100)</td><td>26.1</td><td>26.1</td></tr><tr><td>Fe(3% Si)</td><td>(100)</td><td>13.2</td><td>13.2</td></tr><tr><td>CaF<sub>2</sub></td><td>(111)</td><td>8.62</td><td>11.4</td></tr><tr><td>BaF<sub>2</sub></td><td>(111)</td><td>6.53</td><td>6.52</td></tr><tr><td>Si</td><td>(111)</td><td>16.8</td><td>18.6</td></tr><tr><td>Zn</td><td>(0001)</td><td>13.5</td><td>4.13</td></tr><tr><td>CaCO<sub>3</sub></td><td>(001)</td><td>14.0</td><td>9.45</td></tr></tbody></table>

Cleavage energy measurements are important for Metallurgists, Mineralogists, Electronicians as well as Precious Stone cutters. In the diamond industry, the term used for characterizing a stone is Four C that stand for Color, Clarity (purity), Cut and Carat weight. In the cutting process, one must preserve as much as possible the number of carats. A carat in the Diamond industry is 200 mg which is different from the Gold industry where the carat indicates Gold purity.
## A. Pair and Many-body potentials

* Lennard-Jones potential:

The archetypical pair potential is the Lennard-Jones:

$$ \phi(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \quad (1) $$

where $\epsilon$ is a typical bonding energy and $\sigma$ is a typical atom-atom distance. Other pair potentials are Coulomb interaction, Morse, Buckingham...

* Many-body angular-dependent potentials for BCC metals, covalently bonded systems such as Si, GaAs, Ge...

* Many-body bond order potential where bond strength depends on local environment as in Group IV elements (Si, Ge, C)

<figure><img src="image_2.png" /><figcaption>Fig.1: (Color on-line) Pair potential $\phi(r)$ where $r$ is the pair separation. In order to have bonding one needs an attractive part and a repulsive one with the equilibrium state between the two parts. Interstitials, Vacancies contribution parts to the potential are shown as well as the High-pressure region and the Elastic-Vibration regions around equilibrium.</figcaption></figure>

## B. Reduced Units for pair potentials

When a potential is simple and has only a few parameters like pair potentials, a set of reduced units is recommended for practical reasons.

Reduced units allow to apply results of a single case to different systems described by the same potential with different parameters. Thus evaluation of the potential can be more efficient in reduced units.

Sometimes, reduced units make interpretation of the results more difficult, physical meaning of the results is not immediately clear. For complex potential scaling of the parameters can be not trivial.

Natural choice of length and energy units is $\sigma$ and $\epsilon$. Then, the complete set of parameters of a system in reduced units can be chosen as:

* Length $r^* = r/\sigma$

* Energy $E^* = E/\epsilon$

* Time $t^* = t\sqrt{\epsilon/(m\sigma^2)}$

* Temperature $T^* = k_B T/\epsilon$
* Force $F^* = F\sigma/\epsilon$

* Frequency $\nu^* = \nu\sqrt{m\sigma^2/\epsilon}$

* Pressure $P^* = P\sigma^3/\epsilon$

* Density $\rho^* = \rho\sigma^3/m$

* Surface tension $\gamma^* = \gamma\sigma^2/\epsilon$

### C. Properties of the pair potential

Given a set of N interacting particles, we write the total potential energy:

$$ \phi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3...\mathbf{r}_N) = \sum_{i} \sum_{j>i} \phi(r_{ij}) \quad \text{with} \quad r_{ij} = |\mathbf{r}_i - \mathbf{r}_j| \qquad (2) $$

The forces on individual particles are given by the gradient of the potential energy such that the force on particle i is given by:

$$ F_i = -\frac{\partial\phi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3...\mathbf{r}_N)}{\partial\mathbf{r}_i} \qquad (3) $$

Potential functions like Lennard-Jones have an infinite interaction range. In practice, a cutoff radius $R_c$ is established and interactions between atoms separated by more than $R_c$ are ignored.

Physical quantities (cohesive energy, thermodynamic coefficients...) are impacted by $R_c$ truncation length and most potentials have a built-in cutoff radius limiting the potential range as well as its spatial derivatives.

### D. Pair potential types

Bonding requires a repulsive interaction at short distance and an attraction at long distance resulting in an equilibrium distance between the two parts.

At short-distance, repulsive interaction is chosen as:

* Pauli-Slater type of the form $\exp(-\alpha r)$ with $\alpha$ originating from Pauli principle, overlap and exchange-correlation interactions at short distance between respective electron clouds

* Mie type of the inverse power form $1/r^n$, with n a suitable integer

At large distances, the attractive inverse power form $-1/r^n$ is often chosen.

Commonly used examples of pair potentials:

* Lennard-Jones-Van der Waals interaction in inert gases and molecular systems. Often used to model general effects rather than properties of a specific material.

The Lennard-Jones potential provides a good description of Van der Waals interaction in inert gases and molecular systems (Ar, Kr, CH$_4$, O$_2$, H$_2$, C$_2$H$_4$, etc.).

The term $1/r^{12}$ is repulsive between neighboring atoms. Its physical origin is related to the Pauli principle: when the electronic clouds surrounding the atoms starts to overlap, the energy of the system increases abruptly. Exponent 12 is chosen to make Lennard-Jones interaction easy to compute.

The large distance term $-1/r^6$ is attractive as required to obtain cohesion. An $1/r^6$ attraction describes Van der Waals dipole-dipole interactions due to fluctuating dipoles as shown further below. These very weak interactions are responsible for bonding in closed-shell systems, such as noble gas solids.

Many-body effects are present in inert gases as well and potentials more accurate than Lennard-Jones have been developed for rare gas crystals. Many body effects can account for up to 10 % of total energy.
* Ionic - Coulomb interaction of charges, strong, long range repulsion or attraction. For a pair of charges q₁, q₂ distant by r, we have:

$$ \phi(r) = \frac{q_1 q_2}{4\pi\epsilon_0 r} \qquad (4) $$

* Mie potential

$$ \phi(r) = \epsilon \left[ \left(\frac{\sigma}{r}\right)^m - \left(\frac{\sigma}{r}\right)^n \right] \qquad (5) $$

with m, n appropriate integers.

* Morse (4) potential is similar to Lennard-Jones but is considered somehow more bonding due to presence of fast decaying double exponential terms:

$$ \phi(r) = \epsilon[e^{-2\alpha(r-r_e)} - 2e^{-\alpha(r-r_e)}] \qquad (6) $$

r_e is an equilibrium bonding distance.

It was a popular potential for simulation of metals that have FCC and HCP structures.

* exp-6 (Buckingham) potential: exp term (Slater or Born-Mayer) provides a better description of strong repulsion due to the overlap of the closed shell electron clouds:

$$ \phi(r) = A \exp(-\alpha r) - B \left(\frac{\sigma}{r}\right)^6 \qquad (7) $$

* n-6 potential: The (Pauli-Slater or Born-Mayer) term is replaced by a Mie term:

$$ \phi(r) = A \left(\frac{\sigma}{r}\right)^n - B \left(\frac{\sigma}{r}\right)^6 \qquad (8) $$

* Metals

Pair potentials do not have environmental dependence (e.g. atom in the bulk is too similar to the atom on the surface or near a defect site). In reality, bonding strength should decrease as the local environment becomes too dense resulting from Pauli exclusion principle. Since pair potentials do not depend on the pair surrounding environment, they cannot account for this decrease.

Pair potentials do not account for any directional nature of the bond. Thus covalent contributions (d orbitals) of transition metals that are highly directional cannot be described by pair potentials.

Pair potentials should work for metals where s,p orbitals are present. Unfortunately they overestimate vacancy formation energy, E_v. For instance, E_v = 0.25E_c for Au, E_v = 0.33E_c for Cu, where E_c is the cohesive (bonding) energy.

Another problem is the underestimation of the ratio between cohesive energy and melting temperature. It can be as much as 2-3 times different. Metals must have some other cohesion process that is less effective than a pair potential in keeping the system in the crystalline state.
## E. Dipoles

### 1. Point dipoles

Neutrality of matter imposes charges of opposite signs to reorganize in dipoles, quadrupoles,... or even in multipoles, as for instance in Water where we have formation of octupoles...

A point dipole is given by two opposite charges $\pm q$ distant by $\ell$ and the dipolar moment is $\mathbf{p} = q\ell\mathbf{u}$.

The dipole moment $q\ell$ is expressed practically in Debye units. If we take $q = e$ and $\ell = 1$ Angström, we get a value $\sim 10^{-30}$ in SI units.

Thus Debye suggested a practical unit on that order: 1 Debye is $10^{-10}$ esu $\times 1$ Angström. The result is 1 Debye= $3.3356 \times 10^{-30}$ C.m from the esu defined as follows.

In the CGS-ESU (Electrostatic Units) system, the unit of charge $e'$ is defined from the dyne obeying Coulomb law between two charges in Vacuum distant by $r= 1$ cm. Thus:

$$1 \text{ dyne} = \frac{e'^2}{r^2} \qquad (9)$$

Transforming to SI units: $10^{-5}$ Newton= $\frac{e'^2}{4\pi\epsilon_0 r^2}$, yields $e' = \sqrt{4\pi\epsilon_0 \times 10^{-9}} = 3.3356 \times 10^{-10}$ Coulomb. This means 1 esu of charge = 1 StatCoulomb= $e' = 3.3356 \times 10^{-10}$ Coulomb.

In order to determine the value of the electric charge in SI and CGS units, we consider a spherical shell of radius $R$ carrying a charge $e$. The SI potential energy of the shell is $\frac{1}{2}\frac{1}{4\pi\epsilon_0}\frac{e^2}{R}$, whereas the CGS expression is $\frac{1}{2}\frac{e'^2}{R'^2}$ with $e'$, $R'$ in CGS units. Using the transformation of m into cm, Joule into erg and Coulomb into StatCoulomb, we get $e' = 3 \times 10^9 e \approx ec/10$ StatCoulomb. Thus $e' \approx 4.8 \times 10^{-10}$ StatCoulomb.

Using the previously derived relation between StatCoulomb and the Coulomb, we verify indeed that $e' = 4.8 \times 3.3356 \times 10^{-20} = 1.6 \times 10^{-19}$ Coulomb.

For example, we have some values in Table 2.

<table><caption>**Table 2:** Dipole moments of free molecules and some belonging to ionic crystals. 1 Debye= $3.3356 \times 10^{-30}$ C.m.</caption><thead><tr><th>Substance</th><th>Dipole moment (Debye)</th><th>Dipole moment (C.m)</th></tr></thead><tbody><tr><td>CO<sub>2</sub></td><td>0</td><td>0</td></tr><tr><td>CH<sub>4</sub></td><td>0</td><td>0</td></tr><tr><td>CCl<sub>4</sub></td><td>0</td><td>0</td></tr><tr><td>NH<sub>3</sub></td><td>1.47</td><td>4.9</td></tr><tr><td>H<sub>2</sub>O</td><td>1.85</td><td>6.1</td></tr><tr><td>HF</td><td>1.92</td><td>6.4</td></tr><tr><td>C<sub>6</sub>H<sub>5</sub>NO<sub>2</sub></td><td>4.2</td><td>14</td></tr><tr><td>NaCl</td><td>9</td><td>30.02</td></tr><tr><td>KBr</td><td>10.5</td><td>35.02</td></tr></tbody></table>

### 2. Induced dipoles

**Note:** In this section, we leave out $4\pi\epsilon_0$ factors for simplicity.

The simplest approach to describe induced dipoles is to view them as harmonic oscillators. Thus they are dipolar harmonic oscillators.

Suppose some fluctuation on dipole 1 induces a dipole moment $\mathbf{p}_1$. Then dipole 1 exerts a field on dipole 2:

$$E_1 = \frac{3n(\mathbf{p}_1 \cdot \mathbf{n}) - \mathbf{p}_1}{r^3} \qquad (10)$$
which in turn induces a dipole moment $\mathbf{p}_2 \propto E_1 \propto 1/r^3$. The latter, generates a dipole field $E_2$ felt by $1 \propto \mathbf{p}_2/r^3 \propto 1/r^6$.
Thus, the energy of the interaction is very small and short ranged.

$$W = -\mathbf{p}_1 \cdot \mathbf{E}_2 \propto 1/r^6 \qquad (11)$$

Actually, this is no satisfactory proof since the nature of the fluctuation giving rise to the induced dipole is not explained.

A more appropriate treatment of the van der Waals interaction should rather be quantum and not classical as previously shown.

Quantum effects amount to view induced dipoles as dipolar harmonic oscillators (2).

Thus we consider a model where two identical linear harmonic oscillators 1 and 2 are separated by $R$. Each oscillator bears a dipolar charge $\pm e$ with oscillation amplitudes $x_1$ and $x_2$.

<figure><img src="image_3.png" /><figcaption>Fig.2: (Color on-line) Interacting Dipolar Harmonic oscillators. The spring constant $k = m\omega_0^2$.</figcaption></figure>

The dipoles oscillate along the $x$ axis with frequency $\omega_0$ and momenta $\mathbf{P}_1$ and $\mathbf{P}_2$. Ignoring Coulomb interaction between charges and the internal interaction between the dipolar charges, the Hamiltonian of the system is:

$$H_0 = \frac{P_1^2 + P_2^2}{2m} + \frac{1}{2}m\omega_0^2(x_1^2 + x_2^2) \qquad (12)$$

If we approximate each pair of charges as point dipoles, then they interact with a Hamiltonian:

$$H_1 \approx \frac{-3(\mathbf{p}_2 \cdot \mathbf{n})(\mathbf{p}_1 \cdot \mathbf{n}) + \mathbf{p}_1 \cdot \mathbf{p}_2}{|\mathbf{x}_1 + \mathbf{R} - \mathbf{x}_2|^3} = \frac{2p_1p_2}{R^3} = \frac{2e^2x_1x_2}{R^3} \qquad (13)$$

Note this yields an equivalent spring constant $k = m\omega_0^2$ (see fig. 2).

The total Hamiltonian $H_0 + H_1$ can be diagonalized with a normal mode transformation specifying the symmetric mode (where both oscillate in phase) and the antisymmetric (where they move in opposite phase):

$$x_s = (x_1 + x_2)/\sqrt{2}, \quad x_a = (x_1 - x_2)/\sqrt{2}, \quad P_s = (P_1 + P_2)/\sqrt{2}, \quad P_a = (P_1 - P_2)/\sqrt{2} \qquad (14)$$

After these substitutions, the total Hamiltonian becomes:

$$H = \frac{P_s^2 + P_a^2}{2m} + \frac{1}{2} \left( m\omega_0^2 - \frac{2e^2}{R^3} \right) x_s^2 + \frac{1}{2} \left( m\omega_0^2 + \frac{2e^2}{R^3} \right) x_a^2 \qquad (15)$$

The new eigenfrequencies of these two modes are then

$$\omega_s = \left( \omega_0^2 - \frac{2e^2}{mR^3} \right)^{1/2}, \quad \omega_a = \left( \omega_0^2 + \frac{2e^2}{mR^3} \right)^{1/2} \qquad (16)$$

The zero point energy of the system is:

$$E_0 = \frac{1}{2}\hbar(\omega_s + \omega_a) \approx \hbar\omega_0 \left( 1 - \frac{1}{4} \left( \frac{2e^2}{m\omega_0^2 R^3} \right)^2 + \dots \right) \qquad (17)$$

lowered through the dipole interaction by an amount

$$\Delta U \approx \frac{\hbar\omega_0}{4} \left( \frac{2e^2}{m\omega_0^2 R^3} \right)^2 \qquad (18)$$
typically a small fraction of an eV.

This is the Van der Waals-London interaction or the induced dipole-dipole interaction. It is the principal attractive interaction in crystals of noble gases and in organic molecular crystals.

## F. Many-body potentials for covalently bonded systems

Pair potentials favor close-packed structures whereas most of the covalently bonded crystals have open structures. Silicon forms diamond lattice with four nearest neighbors and has several polytypes under pressure.

Thus many-body potentials should be developed in order to account for many contributions to the potential energy coming from stretching, bending, torsion deformations of bonds and bond angles.

* Stillinger-Weber potential (5)

It is one of the first potentials for diamond lattices (e.g. Si, GaAs, Ge, C). Description of the bonding in Si requires that the potential predicts the diamond lattice (each atom has four neighbors in a tetrahedral arrangement) as the most stable atomic configuration.

It is given by:

$$ \phi = \sum_{i<j} \phi_2(r_{ij}) + \sum_{\substack{i \neq j, k \\ j < k}} \phi_3(r_{ij}, r_{ik}, r_{jk}) \quad (19) $$

where $\phi_2(r_{ij})$ is a two-body term such that we follow the Lennard-Jones parameterization such that:

$$ \phi_2(r_{ij}) = \epsilon f_2(r_{ij}/\sigma) \quad (20) $$

with:

$$ f_2(r) = A(Br^{-p} - r^{-q}) \exp\left(\frac{1}{r-a}\right), \quad r < a, \quad f_2(r) = 0, \quad r \ge a \quad (21) $$

The first term in this expression has a Lennard-Jones form, the second term is a cutoff function that smoothly terminates the potential at a distance a.

The three-body term is given by:

$$ \phi_3(r_i, r_j, r_k) = \epsilon f_3(r_i/\sigma, r_j/\sigma, r_k/\sigma) \quad (22) $$

Directional bonding is introduced in the Stillinger-Weber (5) potential through an explicit three-body (i, j, k) term based on a combinatorial function of angle $\theta_{i,j,k}$:

$$ f_3(\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k) = h(\mathbf{r}_{ij}, \mathbf{r}_{ik}, \theta_{j,i,k}) + h(\mathbf{r}_{ji}, \mathbf{r}_{jk}, \theta_{i,j,k}) + h(\mathbf{r}_{ki}, \mathbf{r}_{kj}, \theta_{i,k,j}) \quad (23) $$

The auxiliary potential function h is limited by the same cutoff length a:

$$ h(\mathbf{r}_{ij}, \mathbf{r}_{ik}, \theta_{j,i,k}) = \lambda \exp[\gamma(r_{ij} - a)^{-1} + \gamma(r_{ik} - a)^{-1}] \left(\cos \theta_{j,i,k} + \frac{1}{3}\right)^2 \quad (24) $$

Thus we have seven parameters: A, B, p, q, a, λ, γ.

The factor $(\cos \theta_{j,i,k} + \frac{1}{3})^2$ favors the perfect tetrahedral structure with angle $\theta_0$ such that $\cos \theta_0 = -\frac{1}{3}$. Thus $\theta_0$ value is given by $\cos^{-1}(-1/3) = 109.47^{\circ}$ in the case of diamond and tetrahedral structure semiconductors like Si, Ge, GaAs...
The main advantage of Stillinger-Weber potential is its simplicity and somehow realistic description of crystalline silicon. However, it has limitations since the three-body term defines only one equilibrium configuration. It is difficult to extend it to an element such as Carbon, that can have several equilibrium angles: 180° (sp hybridization), 120° (sp² hybridization) and 109.5° (sp³ hybridization).

Actually, is possible to rewrite the Stillinger-Weber bending energy in a more appealing way (6) as:

$$ \frac{\beta}{2} \sum_{ijk} (\cos \theta_{ijk} - \cos \theta_0)^2 \qquad (25) $$

leading immediately to the hexagonal value $\theta_0 = 120^\circ$ for graphite and graphene and the tetrahedral angle $\theta_0 = \cos^{-1}(-1/3) = 109.5^\circ$ in the case of Diamond and diamond-like semiconductors such as Si, Ge, GaAs as depicted in fig. 11. Recall that Diamond structure originates from two FCC structures with respective origins at (000) and $(\frac{1}{4}, \frac{1}{4}, \frac{1}{4})$.

### Bond-order potentials

These potentials are more realistic than Stillinger-Weber and have been developed for handling group IV elements (Si, Ge, C) using the concept of bond strength reliance on local environment (Tersoff potential (7) for Si and Ge, Brenner (8) potential for C).

## G. Topological Solids

One of the most exciting recent developments (9; 10) in condensed matter physics and solid-state chemistry is the discovery of topological behavior of matter leading to unique transport properties at lower dimensions.

It is remarkable to observe an important link between topological quantum materials and lithium-ion battery research. The high rate performance of a battery requires the anode to be conductive not just ionically but also electronically. This criterion has significantly stimulated the study of novel conducting anode materials by using topological quantum materials that can intrinsically exhibit high electrical conductivity with high stability protected by topology.

Current progress amounts to using 3D porous topological semimetal anodes composed of light elements such as boron, carbon, and silicon for Li, Na, and K ion batteries. Topological insulators are suggested as battery anodes.

Searching for specific materials (characterized by Atomic identities, Composition and Structure) that have required some target property (here, topological properties, but generally also novel photovoltaic semiconductors, transparent conductors, thermoelectric compounds or battery materials) generally follows an Inverse Design paradigm.

Given a compilation of materials with their characteristics, one first establishes the requisite theoretical conditions or design principles that would enable some topological property P to exist in such crystal structures. Second, one examines if the selected compound and its computed band structure have the stipulated topological property in the assumed crystal structure. Finding a positive answer establishes potential candidates for topological properties. This approach has already provided thousands of specific predictions of topological insulators, or topological crystalline insulators, Dirac and Weyl semi-metals...

Unfortunately, despite world-wide predictions of thousands of topological compounds, very few such compounds have been actually synthesized and proven to be topological.

The problem is due to the assumption the crystal structure must possess the highest symmetry possible, neglecting the possibility of energy-lowering symmetry breaking that destroys topological properties.

## II. INTERACTION POTENTIALS AND ELASTICITY

Robustness of crystal cohesion is checked with respect to departures from equilibrium distances due to disturbances such as vibrations induced by thermal fluctuations, presence of dislocations, impurities, defects, vacancies...
Considering crystal vibration energies, there are many approaches (11) for their description ranging from classical to a hybrid semi-classical mixture between classical and finally to fully quantum.

There exist varieties of stretching, bending, twisting, breathing, rocking... energies ranging from 1D to full 3D and the spatial interaction extent varies from nearest neighbors to remote ones at larger distances such as second, third order neighbors...

Vibrations depend generally on force constants that are measured experimentally from thermodynamic quantities such as specific heat and thermal conductivity. Mechanical measurements yield elastic constants, compressibility and flexural rigidity (bending stiffness or ability to resist bending) given by EI where E is Young modulus and I the moment of inertia. Structural measurements with X-Ray diffraction help in the determination of cohesion energy, equilibrium lattice parameters and angles (2).

Crystals like molecular or Van Der Waals behave like individual molecules and require radial deformations, whereas covalent semiconductors behave like large molecules implying the use of valence models with radial and angular deformations... Other crystals like metals require inclusion of interaction with electrons in the surrounding Fermi sea, magnetic crystals require a distinction between spin up and spin down electrons...

Deformation types are described below.

## A. Stretching deformation energy

Arbitrary deformation of a 3D spring is dealt through writing the elastic potential energy as:

$$V_e = \frac{1}{2} \kappa (\Delta\ell)^2 \qquad (26)$$

where $\Delta\ell$ is the spring elongation and $\kappa$ its elastic constant.

For a set of 3D springs, we generalize this formula as a pairwise ij sum:

$$V_e = \frac{1}{2} \sum_{i,j} \kappa_{ij} \left[ (\mathbf{u}_j - \mathbf{u}_i) \cdot \frac{(\mathbf{r}_i - \mathbf{r}_j)}{|\mathbf{r}_i - \mathbf{r}_j|} \right]^2 \qquad (27)$$

where $\mathbf{u}_i, \mathbf{u}_j$ are the displacements from equilibrium with respect to $\mathbf{r}_i, \mathbf{r}_j$.

In the simple case of a set of 1D aligned springs with different elastic constants, the potential energy is written as:

$$V_e = \frac{1}{2} \sum_{i,j} \kappa_{ij} [(\mathbf{u}_j - \mathbf{u}_i)]^2 \qquad (28)$$

where $\kappa_{ij}$ is the spring constant between sites $i, j$.

## B. Bending deformation energy

Evaluating bending energy entails handling angular distortion with an expansion that ought to be performed in the small distortion case as we did previously in the stretching case.

The simplest stretching and bending energy are given by the Born (11) model:

$$V_e = \frac{1}{2} \sum_{i,j} \alpha_{ij} \left[ (\mathbf{u}_j - \mathbf{u}_i) \cdot \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \right]^2 + \frac{1}{2} \sum_{i,j} \beta_{ij} (\mathbf{u}_j - \mathbf{u}_i)^2 \qquad (29)$$
where $\beta_{ij}$ are bending constants akin to the stretching constants $\alpha_{ij}$ and $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$.

The Keating (11) model contains a more elaborate stretching and bending energies:

$$V_e = \frac{3}{4} \alpha_K \sum_{i,j} \left[ (\mathbf{u}_j - \mathbf{u}_i) \cdot \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \right]^2 + \frac{3}{16} \beta_K \sum_{i,(j,k)} \left[ (\mathbf{u}_j - \mathbf{u}_i) \cdot \frac{\mathbf{r}_{ik}}{|\mathbf{r}_{ik}|} + (\mathbf{u}_i - \mathbf{u}_k) \cdot \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \right]^2 \quad (30)$$

where the indexed sum $i, (j, k)$ means it is performed over triplets of neighboring atoms. $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ and $\mathbf{r}_{ik} = \mathbf{r}_i - \mathbf{r}_k$.

In the case of graphene, another form of bending energy is possible between a triplet (1,2,3) of neighboring atoms given by:

$$\frac{\beta}{2} \left\{ \left[ (\mathbf{u}_2 - \mathbf{u}_1) \times \frac{\mathbf{r}_{12}}{|\mathbf{r}_{12}|^2} \right]_z - \left[ (\mathbf{u}_3 - \mathbf{u}_1) \times \frac{\mathbf{r}_{13}}{|\mathbf{r}_{13}|^2} \right]_z \right\}^2 \quad (31)$$

In fact, it originates from the equilibrium angular constraint of 120°. Notation $|z|$ means out-of-plane component $z$ of the mathematical quantity and $\beta$ is the main bending parameter.

Oshima et al. (12) made further progress by dividing the vibration energy into five terms: nearest neighbor stretching, next nearest neighbor stretching, in-plane bending, out of plane bending and finally twisting (12).

Tewary et al. (13) employed an even more sophisticated approach using Tersoff-Brenner method to rederive the phonon dispersion curves including up to fourth neighbor interactions obtaining results that compared well with experimental measurements (14; 15).

The Tersoff-Brenner (7; 8) potential (TB) has 14 parameters that are obtained from experimentally measured cohesive energy, lattice constant, elastic constants $C_{11}$ and $C_{66}$, flexural rigidity (or bending stiffness) EI (E is Young modulus and I is moment of inertia) and some phonon frequencies.

The TB (7; 8) potential is considered as one of the most accurate and faithful representation of Carbon-Carbon interactions in single, double and triple bonds. The TB potential showed its performance equally well in sp³ hybridized (like Diamond) and sp² hybridized systems (16) (like Graphite and Graphene).

Regarding Graphene elastic constants $C_{11}$ and $C_{66}$, the values obtained with the TB potential are respectively 846 GPa and 248 GPa whereas the values obtained by Tewary et al. (13) are $C_{11}$= 1060 GPa and $C_{66}$= 440 GPa.

### C. Pair potentials and elastic constants

**Note:** In this section, Einstein summation convention is used.

The extension of Hooke law $F = -kx$ of linear elasticity yields proportionality between $\sigma$ and $\epsilon$, $\sigma_{ij} = C_{ij,kl}\epsilon_{kl}$ obtained from the correspondence: $F \leftrightarrow \sigma, -k \leftrightarrow C_{ij,kl}, x \leftrightarrow \epsilon_{kl}$. For a general solid, the 1D spring elastic constant $k$ transforms into a rank-4 tensor $C_{ij,kl}$ containing elastic constants linking $\sigma$ and $\epsilon$ that are both rank-2 tensors.

Elastic energy writes: $U_E = \frac{1}{2}C_{ij,kl}\epsilon_{ij}\epsilon_{kl}$ extending the 1D spring energy definition $U_E = \frac{1}{2}kx^2$. Thus: $\sigma_{ij} = \frac{\partial U_E}{\partial \epsilon_{ij}}$ and since $\epsilon$ is symmetric ($\epsilon_{ij} = \epsilon_{ji}$) we infer $\sigma$ is as well ($\sigma_{ij} = \sigma_{ji}$). Hooke's law is extended to: $C_{ij,kl} = \frac{\partial^2 U_E}{\partial \epsilon_{ij} \partial \epsilon_{kl}}$.

Elastic energy $U_E = \frac{1}{2}C_{ij,kl}\epsilon_{ij}\epsilon_{kl}$ is invariant under $i \leftrightarrow j$ and $k \leftrightarrow l$ interchange. Moreover it is invariant under $\{ij\} \leftrightarrow \{kl\}$ interchange. Thus we infer

$$C_{ij,kl} = C_{ji,kl} = C_{ij,lk} = C_{kl,ij} \quad (32)$$

Thus we may replace a couple of indices $ij$ with a single index $I$ and replace rank-4 tensor $C_{ij,kl}$ with its matrix representation $C_{IJ}$ with $I, J = 1...6$.

The index replacement is done according to the recipe:

$$11 \rightarrow 1, 22 \rightarrow 2, 33 \rightarrow 3, 23 \rightarrow 4, 13 \rightarrow 5, 12 \rightarrow 6 \quad (33)$$
This can be written in a more compact way as:

$$ii \to i, ij \to 9 - (i + j) \text{ when } i \neq j \qquad (34)$$

The 6 × 6 matrix (36 components) represents completely $C_{ij,kl}$ with its $3^4 = 81$ components taking account of symmetry. This matrix is in fact symmetric $C_{ij,kl} = C_{kl,ij}$ originating from the property of elastic energy $\frac{1}{2}C_{ij,kl}\epsilon_{ij}\epsilon_{kl}$ providing another justification for practicality of Voigt notation.

We end up with 21 components (since $[36-6]/2 + 6$, i.e. $15+6=21$) for the triclinic crystal.

Cauchy relation $C_{12}/C_{44} = 1$ is satisfied in polycrystalline isotropic cubic crystals and often satisfied for van der Waals solids and ionic crystals. It is never valid for metals (e.g. $C_{12}/C_{44}$ is 1.5 for Cu, 1.9 for Ag, 3.7 for Au). This means that for van der Waals and ionic solids the elastic constants may be reasonably well described by the pair potential approximation. But for metals pair interaction may be used to represent only part of the total energy.

We are now in position to evaluate Young modulus along a given direction $n$. Applying a stress $\sigma_n$ i.e. a pressure $p$ along direction $n$ induces a strain $\epsilon_n$ along same direction yielding Young modulus $E(n)$:

$$\frac{1}{E(n)} = \frac{\epsilon_n}{\sigma_n} \qquad (35)$$

Inverse transform of $\sigma_{ij} = C_{ij,kl}\epsilon_{kl}$, into $\epsilon_{ij} = S_{ijkl}\sigma_{kl}$ defines $S_{ijkl}$ the compliance tensor yielding inverse of relation $C_{ij,kl}\sigma_{kl}$.

$S_{ijkl}$ is the inverse tensor of $C_{ij,kl}$ and inherits its symmetry properties.

This allows us to relate $\epsilon_n$ to $\sigma_n$:

$$\sigma_n = \sigma_{ij} n_i n_j, \quad \epsilon_n = \epsilon_{ij} n_i n_j \qquad (36)$$

Introducing $\sigma_{ij} = P n_i n_j$ with $P$ the pressure into $\sigma_n = \sigma_{ij} n_i n_j$, we get the Young modulus reciprocal as:

$$\frac{1}{E(n)} = \frac{\epsilon_n}{\sigma_n} = \frac{P S_{ijkl} n_i n_j n_k n_l}{P} = S_{ijkl} n_i n_j n_k n_l \qquad (37)$$

This illustrates how $E_\parallel$ Young modulus parallel to cleavage plane and $E_\perp$ Young modulus perpendicular to cleavage plane are evaluated for the measurement of the cleavage energy.

### III. LATTICE SUMS

Performing lattice sums has progressed tremendously (17) analytically and numerically to tackle particularly difficult alternating sums that suffer from extremely slow convergence.

Evjen method handled alternating sums by building neutral blocks whereas Ewald summation relied on error and Jacobi theta functions (18). More recently, multiple precision, Fourier based techniques and advanced analytical transforms (17) were developed with tremendous success.

Starting with periodic 1D linear chain, the alternating sum based Madelung constant reads:

<figure><img src="image_4.png" /><figcaption>Fig.3: (Color on-line) Periodic linear chain displaying first, second and third neutral cells (red boxes).</figcaption></figure>
$$ \alpha_M = - \sum_{i \neq 0} \frac{(-)^i}{p_{0i}} = 2 \left( 1 - \frac{1}{2} + \frac{1}{3} - \dots \right) \qquad (38) $$

where 0 is the central negative ion.

The above sum can be related to the expansion for $|x| < 1$:

$$ \frac{1}{1+x} = 1 - x + x^2 - x^3 \dots \qquad (39) $$

Integrating the expansion we get:

$$ \int dx \frac{1}{1+x} = \ln(1+x) = x - \frac{1}{2}x^2 + \frac{1}{3}x^3 \dots \qquad (40) $$

The expansion actually converges for $x = 1$ thus we get:

$$ \ln(2) = 1 - \frac{1}{2} + \frac{1}{3} \dots \qquad (41) $$

and the Madelung constant becomes $2 \ln 2 \sim 1.4$.

Numerically it is possible to transform the above alternating series into a non-alternating one to speed convergence.

The series:

$$ \ln \left( \frac{1+x}{1-x} \right) = 2 \left( x + \frac{x^3}{3} + \frac{x^5}{5} + \frac{x^7}{7} \dots \right) \qquad (42) $$

yields an approximate value of $2 \ln 2$ when $x = 1/3$ with a few terms.

Alternating sums have very slow convergence rate. Evjen (19) suggested to arrange charges into neutral cells to accelerate convergence.

Let us illustrate this with the periodic linear chain.

First neutral cell yields: $\alpha_M \sim 2 \times 1/2$. Second neutral cell yields: $\alpha_M \sim 2 \times (1/2 + 1/2 - 1/4)$. Third neutral cell gives: $\alpha_M \sim 2 \times (1/2 + 1/2 - 1/4 - 1/4 + 1/6) \sim 1.333$ which is close to the exact result $\sim 1.4$.

We apply this method to 2D and higher dimensional crystals.

For a square lattice, the Madelung constant given by:

$$ \alpha_M = - \sum_{i \neq j} \frac{(-)^{i+j}}{\sqrt{i^2 + j^2}} \qquad (43) $$

The primed sum means that the origin $i = j = 0$ must be excluded.

The first neutral cell gives: $4 \times 1/2 - 4 \times 1/4\sqrt{2} = 1.293$.

The second neutral cell yields $4 \times 1/2 - 4 \times 3/4\sqrt{2} - 4/(2.2) + 8 \times 1/2\sqrt{5} - 4 \times 1/4\sqrt{8} = 0.314$.

The third neutral cell gives in the same way: $3.6 \times 10^{-3}...$

Thus we get $\alpha_M \sim 1.61$. A more accurate value is 1.6155.

In the case of a 3D cubic lattice, Madelung constant is defined by:

$$ \alpha_M = - \sum_{i \neq j \neq k} \frac{(-)^{i+j+k}}{\sqrt{i^2 + j^2 + k^2}} \qquad (44) $$

and its numerical value is 1.7475 (see Table 7).
<figure><img src="image_5.png" /><figcaption>Fig.4: (Color on-line) Square ionic lattice. The first and second neutral cells are shown in red around the central negative charge.</figcaption></figure>

A. Ewald summation

In 1921 Ewald proposed in his thesis an efficient way to recast the Coulomb summation involving ion-ion interaction between charged particles in two fast converging series.

Consider a basic cell containing a positive and a negative charge as a single crystallographic element and make an infinite number of replicas around it. This is a periodic extension of the basic cell paving the way to using Fourier series techniques.

The entire system is then neutral and contains an infinite number of charges situated at points $\mathbf{r}_j^+$ and $\mathbf{r}_j^-$, respectively.

The total potential at the position of some ion residing in the basic cell is given by the finite difference of two infinite, diverging series (20):

$$ \phi(\mathbf{r}_i) = q \sum_{j=1}^{\infty} \frac{1}{|\mathbf{r}_i - \mathbf{r}_j^+|} - q \sum_{j=1}^{\infty} \frac{1}{|\mathbf{r}_i - \mathbf{r}_j^-|} $$

In order to avoid divergences, we represent these charges as delta-like charge densities,

$$ \rho(\mathbf{r}) = q \sum_{j=1}^{\infty} \delta(\mathbf{r} - \mathbf{r}_j^+) - q \sum_{j=1}^{\infty} \delta(\mathbf{r} - \mathbf{r}_j^-) $$

and expand the result in a Fourier series whose components are $\phi(\mathbf{k})$ the electrostatic potential.

Since the Fourier representation of a delta-function requires infinitely many terms, the Fourier space calculation would again lead to convergence problems.

The solution is to split up the potential in two well-behaved parts, one being represented in $\mathbf{r}$-space and the other in $\mathbf{k}$-space by rapidly converging series.

Let us illustrate this with a one-dimensional ion lattice.

We add to the delta-like point charges, Gaussian distributed charges of opposite sign, to form an auxiliary lattice:

$$ \rho'(\mathbf{r}) = -q_j \left( \frac{\eta_E^2}{\pi} \right)^{3/2} e^{-\eta_E^2 (\mathbf{r} - \mathbf{r}_j)^2} \quad (45) $$

where $\eta_E$ is an adjustable parameter introduced by Ewald.
Another lattice is introduced to compensate (20) the additional Gaussian charges, such that the contribution of each lattice to the potential are evaluated respectively in **r** and **k** spaces:

* Lattice 1:
  It is calculated in **r**-space, where the series converges very rapidly. The convergence is quicker for narrow Gaussians, i.e. when $\eta_E$ is large.

* Lattice 2:
  It is calculated in **k**-space. When the Gaussians are broad, i.e. when $\eta_E$ is small, a smaller number of Fourier components are required.

When $\eta_E$ is properly chosen, optimal convergence of both series may be achieved.

Extending this approach to three-dimensional model systems, we consider a cubic base cell with side length $L$ containing charges with wavevectors:

$$ \mathbf{k} = \frac{2\pi}{L} (k_x, k_y, k_z) $$

and interparticle vectors including all periodic images of the base cell, thus (20):

$$ \mathbf{r}_{i,j,n} = \mathbf{r}_j + n\mathbf{L} - \mathbf{r}_i \quad (i, j = 1, \dots, N) $$

where $n\mathbf{L}$ is a translation vector in the periodic lattice.

Ewald sum is obtained with the following operations:

$$ \phi(\mathbf{r}_i) = \frac{4\pi}{L^3} \sum_{j=1}^{N} q_j \left[ \sum_{\mathbf{k}} e^{-i\mathbf{k}\cdot\mathbf{r}_{ij}} k^{-2} e^{-k^2/4\eta_E^2} + \sum_{\mathbf{n}} \text{erfc}(\eta_E |\mathbf{r}_{i,j,n}|) \right] \quad (46) $$

with

$$ \text{erfc}(z) = \frac{2}{\sqrt{\pi}} \int_z^\infty e^{-t^2} dt \quad (47) $$

the complementary error function.

Considering these two corrections, the total Coulomb energy is:

$$ U = \frac{1}{2} \sum_{i=1}^{N} q_i \phi(\mathbf{r}_i) - \frac{\eta_E}{\sqrt{\pi}} \sum_{i=1}^{N} q_i^2 + \frac{2\pi}{3L^3} \left| \sum_{i=1}^{N} q_i \mathbf{r}_i \right|^2 \quad (48) $$

Thus it becomes possible to use Ewald summation to evaluate Madelung constants.

Table 3 displays values of the Madelung constant $\alpha_a$ normalized by the lattice parameter $a$ for several crystal structures and the occurrence of large values of $\alpha_a$.

## B. Theta function summation

Ewald introduced in his thesis in 1921 error functions as well as Jacobi theta functions allowing to calculate lattice sums analytically without having to calculate coordination numbers and distances as one goes from first to second, third... nearest neighbors.

We illustrate this approach with the SC structure.
<table><caption>**Table 3:** Madelung constants for several 3D Ionic crystals. *AR* is axial ratio 1 : *c* where *c* is the axial hexagonal lattice parameter. Note the large values of the Madelung constant in the Spinels (2) that are cubic ferrites (magnetic oxides). Madelung constant is normalized by the lattice parameter. In the Spinel (2) structure, *pe* and *qe* are the charges of the P and Q cations at the Tetrahedral and Octahedral sites. The anion charge is assumed to be 2e and *p* + 2q = 8. Adapted from Borwein *et al.* (17)</caption><thead><tr><th>Crystal</th><th>Variety</th><th>α<sub>a</sub></th></tr></thead><tbody><tr><td>NaCl</td><td></td><td>3.495129189</td></tr><tr><td>CsCl</td><td></td><td>2.035361510</td></tr><tr><td>ZnS</td><td></td><td>3.782926104</td></tr><tr><td>CaF<sub>2</sub></td><td></td><td>11.63657523</td></tr><tr><td>Cu<sub>2</sub>O</td><td></td><td>10.25945703</td></tr><tr><td>KZnF<sub>3</sub></td><td></td><td>12.37746803</td></tr><tr><td>LaAlO<sub>3</sub></td><td></td><td>44.55497526</td></tr><tr><td rowspan="3">CaWO<sub>4</sub></td><td>(AR = 1)</td><td>2.254775948</td></tr><tr><td>(AR = √2)</td><td>2.282508448</td></tr><tr><td>(AR = 2)</td><td>2.284666496</td></tr><tr><td>ReO<sub>3</sub></td><td></td><td>71.63183493</td></tr><tr><td>BaTiO<sub>3</sub></td><td></td><td>49.50987213</td></tr><tr><td>NaTaO<sub>3</sub></td><td></td><td>58.53549202</td></tr><tr><td>YOF</td><td></td><td>27.05607656</td></tr><tr><td>BiF<sub>3</sub></td><td></td><td>22.12196280</td></tr><tr><td>BaLiF<sub>3</sub></td><td></td><td>10.91770035</td></tr><tr><td>NbO</td><td></td><td>3.008539964</td></tr><tr><td>Li<sub>2</sub>O</td><td></td><td>10.77318449</td></tr><tr><td rowspan="2">PdO</td><td>(AR = √2)</td><td>3.649338711</td></tr><tr><td>(AR = 2)</td><td>4.539442445</td></tr><tr><td>ZrI<sub>4</sub></td><td></td><td>27.50917245</td></tr><tr><td rowspan="3">PtS</td><td>[Pt(0, 0, 0) S(1/4,1/4,1/4)]</td><td>2.636813487</td></tr><tr><td>[Pt(0, 0, 0) S(1/2,0,0)]</td><td>2.741365176</td></tr><tr><td>[Pt(0, 0, 0) S(1/2,0,1/2)]</td><td>2.254775948</td></tr><tr><td>K<sub>2</sub> PtCl<sub>6</sub></td><td></td><td>63.49038889</td></tr><tr><td>(NH<sub>4</sub>)<sub>3</sub> AlF<sub>6</sub></td><td></td><td>53.00500132</td></tr><tr><td>CuVS<sub>4</sub></td><td></td><td>66.7566350</td></tr><tr><td rowspan="8">Spinels</td><td>Normal</td><td></td></tr><tr><td>p = 2, q = 3</td><td>128.5671125</td></tr><tr><td>p = 3, q = 5/2</td><td>130.7743620</td></tr><tr><td>p = 4, q = 2</td><td>138.1991295</td></tr><tr><td>Inverse</td><td></td></tr><tr><td>p = 4, q = 2</td><td>132.5694531</td></tr><tr><td>p = 2, q = 3</td><td>131.7749472</td></tr></tbody></table>

Let us consider the following sum:

$$ \zeta_C(s) = \frac{1}{2} \sum'_{i,j,k} \frac{1}{(i^2 + j^2 + k^2)^{s/2}} \quad (49) $$

Using the Γ identity (21):

$$ \frac{1}{r^s} = \frac{1}{\Gamma(s/2)} \int_0^\infty dt \, t^{s/2-1} e^{-r^2 t} = \frac{1}{\Gamma(s/2)} \int_0^\infty \frac{dt}{t} t^{s/2} e^{-r^2 t} \quad (50) $$

where $r = \sqrt{i^2 + j^2 + k^2}$, the $\zeta_C(s)$ sum can be rewritten as:
$$ \zeta_C(s) = \frac{1}{2\Gamma(s/2)} \int_0^\infty \frac{dt}{t} t^{s/2} \left[ \sum_{i,j,k=-\infty}^\infty e^{-(i^2+j^2+k^2)t} - 1 \right] \quad (51) $$

The above sum can be expressed by a Jacobi theta (21) function $\theta_3$ as:

$$ \zeta_C(s) = \frac{1}{2\Gamma(s/2)} \int_0^\infty \frac{dt}{t} t^{s/2} [\theta_3^3(e^{-t}) - 1] \quad (52) $$

Jacobi function (18) $\theta_3(q)$ is defined as $\theta_3(q) = \sum_{n=-\infty}^{\infty} q^{n^2}$. $q \in ]-1, 1[$ is the theta function modulus.

In order to perform the lattice sum over a BCC lattice, we infer that it is made of two SC lattices shifted one versus other by half the lattice parameter along all three directions (2).

Consequently, the sum becomes:

$$ \zeta_B(s) = \frac{1}{2^{3/2}} [\zeta_C(s) + \psi(1/2, 1/2, 1/2; s)] \quad (53) $$

where the auxiliary function $\psi(x_1, x_2, x_3; s)$ is given by:

$$ \psi(x_1, x_2, x_3; s) = \frac{1}{2} \sum_{i,j,k=-\infty}^{\infty} \frac{1}{[(i-x_1)^2 + (j-x_2)^2 + (k-x_3)^2]^{s/2}} \quad (54) $$

Transforming with the $\Gamma$ identity, as previously, we get:

$$ \zeta_B(s) = \frac{1}{2^{3/2}} \left[ \zeta_C(s) + \frac{\pi^{s/2}}{2\Gamma(s/2)} \int_0^\infty \frac{dt}{t} t^{s/2} \theta_2^3(e^{-\pi t}) \right] \quad (55) $$

Jacobi function (18) $\theta_2(q)$ is defined by $\theta_2(q) = \sum_{n=-\infty}^{\infty} q^{(n-1/2)^2}$.

Madelung constants are calculated with Jacobi function (18) $\theta_4(q)$ defined by $\theta_4(q) = \sum_{n=-\infty}^{\infty} (-)^n q^{n^2}$.

As an example, we consider the $D$-dimensional (17) simple cubic lattice generalized (since the actual Madelung constant is evaluated with $s = 1/2$) Madelung sum:

$$ M_D(s) = \frac{1}{2} \sum'_{i_1, i_2, i_3, \dots, i_D} \frac{(-1)^{i_1+i_2+i_3+\dots+i_D}}{(i_1^2 + i_2^2 + i_3^2 + \dots + i_D^2)^s} \quad (56) $$

Using the $\Gamma$ identity allows expressing the result in terms of the $\theta_4(q)$ yielding (17) the following results for $D = 2, 4, 6; 8$:

$$ M_2(s) = -4\beta(s)\eta(s), \quad M_4(s) = -8\eta(s)\eta(s-1), \\ M_6(s) = 4\beta(s-2)\eta(s) - 16\beta(s)\eta(s-2) \quad M_8(s) = -16\zeta(s)\eta(s-3) \quad (57) $$

where $\zeta = \sum_{n=0}^{\infty} (n+1)^{-s}$ is the Riemann zeta function, $\eta(s) = \sum_{n=0}^{\infty} (-)^n (n+1)^{-s}$ and $\beta(s) = \sum_{n=0}^{\infty} (-)^n (2n+1)^{-s}$ are the Dirichlet eta and beta functions.

Thus it is possible to extend, in principle analytically, the Madelung sum $M_D(s)$ to any dimension and in particular to $D = 100$ as depicted in fig. 5 for $s = 1/2$:

The graph in fig. 5 is obtained with a combination of analytical and numerical procedures since large powers $n$ of the $\theta_4$ functions are required. It is possible to express $\theta_4^n(q)$, when $n$ is small, in other theta functions combined with Dirichlet eta and beta functions. Unfortunately, when $n$ is large we have to rely on numerical methods (22) to evaluate $\theta_4^n(q)$.
<figure><img src="image_6.png" /><figcaption>**Fig.5:** (Color on-line) Madelung constant $M_D(s = 1/2)$ for the simple cubic lattice as a function of space dimension $D$ up to 100. Adapted from Burrows *et al.* (22)</figcaption></figure>

# IV. CLASSICAL THERMODYNAMICS

## A. Thermodynamic potential functions

Thermodynamics is based on potential functions: Gibbs, Enthalpy, Energy and Free energy respectively $G, H, U$ and $F$ and their derivatives with respect to a set of variables: $P, T, S, V$ Pressure, Temperature, Entropy and volume respectively with other variables intervening depending on the physical system such as electric field and polarization or magnetic field and magnetization...

The Jacobian method allows one to derive easily any thermodynamic identity as long as one respects the canonical variables associated with some potential.

This is provided with the Guggenheim (23) table below where each potential Gibbs, Enthalpy, Energy and Free energy respectively $G, H, U$ and $F$ is flanked on each side by the appropriate canonical variables.

<figure><img src="image_7.png" /><figcaption>**Fig.6:** (Color on-line) At left, the table is constructed with a simple sentence: Great Physicists Have Studied Under Very Fine Teachers starting with $G$ at left of middle row and proceeding clockwise. The corner blue-highlighted variables make Maxwell relation table.</figcaption></figure>

If we want to derive the enthalpy differential $dH$ from the table, we write $dH = adP + bdS$ and $a, b$ are found by inspecting the element in box across the diagonal from the box containing the canonical variable. We have a (+) sign if we are moving down from the canonical variable to reach out for the diagonally opposite corner box and a (-) sign in the opposite case. Thus $a = V$ and $b = T$ yielding $dH = VdP + TdS$. In the energy differential case $dU$ we have $dU = adS + bdV$ with $a = T, b = -P$ yielding $dU = TdS - PdV$. In the case of free energy we have $dF(T, V) = -SdT - PdV$ whereas for Gibbs potential, $dG(P, T) = VdP - SdT$.

Maxwell identities are constructed as well from the table corners. They are given by the mnemonic matrix:

$$ \left| \begin{matrix} PS \\ TV \end{matrix} \right| $$ yielding: $(\frac{\partial P}{\partial T})_V = (\frac{\partial S}{\partial V})_T$ and each time we exchange variables diagonally we pick up a sign change as in the case $P \leftrightarrow V$:
$$ \left| \begin{array}{c} VS \\ TP \end{array} \right| \text{ equivalent to: } \left( \frac{\partial V}{\partial T} \right)_P = - \left( \frac{\partial S}{\partial P} \right)_T. $$

Note that if we want to enlarge the number of thermodynamic variables to the Grand Canonical ensemble case that includes $N, \mu$, we write:

$$ \bullet dU(S, V, N) = TdS - PdV + \mu dN, $$

$$ \bullet dH(P, S, N) = VdP + TdS + \mu dN, $$

$$ \bullet dF(T, V, N) = -SdT - PdV + \mu dN $$

$$ \bullet dG(P, T, N) = VdP - SdT + \mu dN. $$

From $dU(S, V, N)$. we obtain $\mu = -T(\frac{\partial U}{\partial N})_{S,V}$ or $\mu = -T(\frac{\partial S}{\partial N})_{V,U}$ by taking into $dU(S, V, N) = TdS - PdV + \mu dN$ $V, U$ as constants implying $dV = 0, dU = 0$.

Same can be done with the other thermodynamic potentials $H, F, G$.

In order to describe the use of Jacobians in Thermodynamics, let us recall some of their properties (24; 25):

$$ \bullet \frac{\partial(u,v)}{\partial(x,y)} = \begin{vmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{vmatrix} $$

$$ \bullet \frac{\partial(u,v)}{\partial(x,y)} = -\frac{\partial(v,u)}{\partial(x,y)} $$

$$ \bullet \frac{\partial(u,y)}{\partial(x,y)} = \left( \frac{\partial u}{\partial x} \right)_y $$

$$ \bullet \frac{\partial(u,v)}{\partial(x,y)} = \frac{\partial(u,v)}{\partial(s,t)} \frac{\partial(s,t)}{\partial(x,y)} $$

$$ \bullet \frac{\partial(u,v)}{\partial(x,y)} = 1 / \frac{\partial(x,y)}{\partial(u,v)} $$

This means that any term such as $\partial(x, y)$ behaves as if it were an elementary algebraic coefficient simplifying and speeding up enormously mathematical operations.

Nonetheless, there is a limit to express any physical quantity: It should depend on three experimentally measurable quantities (24):

1. Dilation coefficient $\alpha_P = \frac{1}{V} \left(\frac{\partial V}{\partial T}\right)_P$

2. Compressibility $\kappa_T = -\frac{1}{V} \left(\frac{\partial V}{\partial P}\right)_T$

3. Heat capacity $C_V = T \left(\frac{\partial S}{\partial T}\right)_V$ or $C_P = T \left(\frac{\partial S}{\partial T}\right)_P$

For illustration of these notions, let us evaluate several thermodynamic quantities:

1. Molar enthalpy of sublimation $H_s$:

$H_s = \frac{1}{2} z N_A \epsilon$ where $z$ is the number of first-nearest neighbors (Coordination Number), $N_A$ Avogadro number and $\epsilon$ the cleavage energy.

2. Linear dilation coefficient $\alpha_P$:

Using the picture of a particle moving in a 1D potential $\phi(x)$, we infer that if $\phi(x)$ is harmonic, the particle will always return to its equilibrium position at $x = 0$ with no dilation whatsoever.

Thus we consider a weakly anharmonic potential $\phi(x) = Cx^2 - gx^3 + fx^4$, with $C \gg f, g$, the average position $\langle x \rangle$, in the Canonical Ensemble (2), is given at temperature $T$ by:

$$ \langle x \rangle = \frac{\int_{-\infty}^{+\infty} x \exp(-\beta\phi(x))dx}{\int_{-\infty}^{+\infty} \exp(-\beta\phi(x))dx}; \quad \beta = \frac{1}{k_B T} \qquad (58) $$
<figure><img src="image_8.png" /><figcaption>**Fig.7:** (Color on-line) Harmonic potential $\phi(x) = Cx^2$ in green and anharmonic potential $\phi(x) = Cx^2 - gx^3 + fx^4$ in red displaying a local minimum at positive $x$ value implying dilation. The condition for local minimum existence is $9g^2 > 32Cf$.</figcaption></figure>

Since $C \gg f, g$ we expand the exponential to first order:

$$ \exp(-\beta\phi(x)) = \exp(-\beta C x^2) \exp(\beta g x^3 - \beta f x^4) \approx \exp(-\beta C x^2)[1 + \beta g x^3 - \beta f x^4] \quad (59) $$

The integration domain being symmetric, the even term that survives in the $\langle x \rangle$ numerator is:

$$ \beta g \int_{-\infty}^{+\infty} x^4 \exp(-\beta C x^2) dx \quad (60) $$

whereas in the $\langle x \rangle$ denominator the surviving term is:

$$ \int_{-\infty}^{+\infty} \exp(-\beta C x^2) dx \quad (61) $$

Making a change of variable $y = \beta C x^2$ yields:

$$ \langle x \rangle = \frac{\beta g \int_{0}^{+\infty} y^{\frac{3}{2}} \exp(-y) dy}{(\beta C)^2 \int_{0}^{+\infty} y^{-\frac{1}{2}} \exp(-y) dy} \quad (62) $$

Using $\Gamma(x)$ function properties (18):

$$ \Gamma(x) = \int_{0}^{\infty} e^{-t} t^{x-1} dt \quad [x > 0], \quad \Gamma(x+1) = x\Gamma(x), x \in \mathbb{R}; \quad \Gamma(n) = (n-1)! n \in \mathbb{N}, \quad \Gamma\left(\frac{1}{2}\right) = \sqrt{\pi} \quad (63) $$

we obtain:

$$ \langle x \rangle = \frac{\beta g \Gamma\left(\frac{5}{2}\right)}{(\beta C)^2 \Gamma\left(\frac{1}{2}\right)} = \frac{3g}{4C^2} k_B T, \quad \text{since} \quad \Gamma\left(\frac{5}{2}\right) = \frac{3}{2}\Gamma\left(\frac{3}{2}\right) = \frac{3}{4}\Gamma\left(\frac{1}{2}\right) \quad (64) $$
The dilation coefficient is given by:

$$ \alpha_P = \frac{1}{V} \left[ \frac{\partial V}{\partial T} \right]_P = \frac{\partial \langle x \rangle}{\partial T} = \frac{3g}{4C^2} k_B \quad (65) $$

In the Fe and Ni cases $\alpha_P \sim 10 \times 10^{-6}/\text{K}$, whereas in the Invar (Ni₃₆Fe₆₄ alloy) case, $\alpha_P$ remains smaller that $10^{-6}/\text{K}$ from low temperature up to 90° Celsius, making it the best alloy for the mechanical clock industry.

### 3. Difference between heat capacities $C_P - C_V$:

Let us now move on to find the expression of $C_P - C_V$. We start by writing:

$$ C_V = T \left( \frac{\partial S}{\partial T} \right)_V = T \frac{\partial(S, V)}{\partial(T, V)} = \frac{\partial(S, V)}{\partial(T, P)} \frac{\partial(T, P)}{\partial(T, V)} \quad (66) $$

which gives:

$$ C_P - C_V = -T \left( \frac{\partial V}{\partial T} \right)_P^2 / \left( \frac{\partial V}{\partial P} \right)_T = TV \alpha_P^2 / \kappa_T \quad (67) $$

In the ideal gas case we have the state equation $PV \propto T$ resulting in $\alpha_P = 1/T$ and $\kappa_T = 1/P$. Thus $C_P - C_V = PV/T = Nk_B$ where N is the number of particles.

### 4. Adiabatic bulk modulus $B_S$:

To find the expression of $B_S$, we start by writing: $B_S = -V \left(\frac{\partial P}{\partial V}\right)_S$ using Jacobians and Maxwell relations.

This leads to $B_S = -V \left(\frac{\partial P}{\partial V}\right)_S$ yielding:

$$ B_S = -V \frac{\partial(P, S)}{\partial(V, S)} = -V \frac{\partial(P, S)}{\partial(V, T)} \frac{\partial(V, T)}{\partial(V, S)} = -V \left| \begin{pmatrix} \frac{\partial P}{\partial V} \\ \frac{\partial S}{\partial V} \end{pmatrix}_T \begin{pmatrix} \frac{\partial P}{\partial T} \\ \frac{\partial S}{\partial T} \end{pmatrix}_V \right| \left( \frac{\partial T}{\partial S} \right)_V \quad (68) $$

Using Jacobian properties gives:

$$ B_S = -V \left| \begin{matrix} \frac{1}{V \kappa_T} \left( \frac{\partial P}{\partial T} \right)_V & T \\ \left( \frac{\partial S}{\partial V} \right)_T & \frac{C_V}{T} \end{matrix} \right| \quad (69) $$

After developing the determinant:

$$ B_S = -\frac{VT}{C_V} \left[ -\frac{C_V}{VT\kappa_T} - \left(\frac{\partial P}{\partial T}\right)_V \left(\frac{\partial S}{\partial V}\right)_T \right] \quad (70) $$

Maxwell relations yield:

$$ \left( \frac{\partial S}{\partial V} \right)_T = \left( \frac{\partial P}{\partial T} \right)_V $$

and the Jacobian method gives:

$$ \left(\frac{\partial P}{\partial T}\right)_V = -\left(\frac{\partial V}{\partial T}\right)_P \left(\frac{\partial P}{\partial V}\right)_T = \alpha_P V \frac{1}{VT\kappa_T} = \frac{\alpha_P}{\kappa_T} \quad (71) $$
Thus

$$B_S = \frac{VT}{C_V} \left[ \frac{C_V}{VT\kappa_T} + \left( \frac{\partial P}{\partial T} \right)_V^2 \right] \qquad (72)$$

and finally:

$$B_S = B_T + \frac{VT}{C_V} \alpha_P^2 B_T^2.$$

Since for a solid the dilation coefficient is relatively small, consequently: $B_S \sim B_T$.

## B. Non-ideal gas thermodynamics

In a real gas, attractive and repulsive interactions contribute to pressure, modifying the ideal gas equation of state $PV = nRT$.

In order to find the pressure in presence of interactions, we use the grand partition function written as a series in the fugacity (26) $z = \exp(\mu/k_B T)$ where $\mu$ is the chemical potential,

$$Z = \mathrm{Tr} \exp(-(H - \mu \hat{N})/k_B T) = \sum_{N=0}^{\infty} \exp(\mu N/k_B T) \mathrm{Tr} \exp(-H_N/k_B T) = \sum_{N=0}^{\infty} z^N Z_N, \qquad (73)$$

where $H_N$ is the $N$ particle Hamiltonian, and $Z_N = \mathrm{Tr} \exp(-H_N/k_B T)$ the canonical partition function.

The thermodynamic potential (26),

$$\Omega = -k_B T \ln Z = -k_B T \ln \left( \sum_{N=0}^{\infty} z^N Z_N \right) = -k_B T \ln \left( 1 + \sum_{N=1}^{\infty} z^N Z_N \right) \qquad (74)$$

Thus expansion of $\Omega$ (16) in fugacity powers is possible, moreover it can equally be done in density powers as in the Virial expansion.

Let us consider the thermodynamics of a real gas whose Hamiltonian contains pair interaction $\phi(r)$ depending on the distance $r$ between the gas atoms,

$$H_N = \sum_{n=1}^{N} \frac{p_n^2}{2m} + \sum_{n<m}^{N} \phi(|\mathbf{r}_n - \mathbf{r}_m|), \qquad (75)$$

where $p_n$ and $r_n$ are the momentum and position operators and $m$ the atom mass, and the sum over pairs of atoms avoiding double counting contain $N(N-1)/2$ terms.

The grand partition function (26):

$$\begin{aligned} Z(T, \mu, V) &= \sum_{N=0}^{\infty} Z_N(T, V) z^N \\ &= 1 + \sum_{N=1}^{\infty} \frac{M_N(T, V)}{N!} z^N, \end{aligned} \qquad (76)$$

where $M_N$ are moments defined by,

$$\frac{M_N(T, V)}{N!} = \mathrm{Tr} \exp(-H_N/k_B T) = \sum_n \exp(-E_n^{(N)}/k_B T), \qquad (77)$$
where $E_n^{(N)}$ are the energy eigenvalues of the N-particle Hamiltonian $H_N$.
Consequently the grand potential $PV$ of the dilute gas (16) is given by:

$$PV = k_B T \ln \Omega \qquad (78)$$

resulting in the equation of state for the grand potential $PV$ as a series function of $T, V, N$:

$$PV = Nk_B T \left( 1 + B_2(T) \frac{N}{V} + B_3(T) \frac{N^2}{V^2} + \dots \right) \qquad (79)$$

with the Virial coefficients $B_N(T)$. In particular, the second order virial coefficient is:

$$B_2(T) = -2\pi \int_0^\infty r^2[\exp\left(-\frac{\phi(r)}{k_B T}\right) - 1]dr \qquad (80)$$

As an example, for a hard-core potential $\phi_H(r)$ given by:

$$\phi_H(r) = \infty, \quad r < \sigma, \quad \phi_H(r) = 0, \quad r > \sigma \qquad (81)$$

where $\sigma$ is a minimal approach distance, the second order virial coefficient $B_2(T)$ is evaluated by:

$$B_2(T) = -2\pi \int_0^\infty r^2[\exp\left(-\frac{\phi_H(r)}{k_B T}\right) - 1]dr = 2\pi \int_0^\sigma r^2 dr = \frac{2\pi\sigma^3}{3} \qquad (82)$$

A minimal approach length $\sigma$ implies an excluded volume given by $B_2(T) = \frac{2\pi\sigma^3}{3}$.

The equation of state (to second order in density) is:

$$PV = Nk_B T \left( 1 + B_2(T) \frac{N}{V} \right). \qquad (83)$$

For a Lennard-Jones system (fig. 8), the second virial coefficient is evaluated as:

<figure><img src="image_9.png" /><figcaption>Fig.8: (Color on-line) At left, Lennard-Jones interaction potential $\phi(r)$ with its basic parameters $\sigma$ and $\epsilon$. Energy minimum $-\epsilon$ is at bonding distance $r^* = \sqrt{2}\sigma$. At right, Sutherland potential with hard-core repulsion and $-1/r^6$ attraction.</figcaption></figure>

$$B_2(T) = -2\pi \int_0^\infty r^2[\exp\left(-\frac{4\epsilon}{k_B T}\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]\right) - 1]dr \qquad (84)$$
<figure><img src="image_10.png" /><figcaption>**Fig.9:** (Color on-line) Example of $B_2(T)$ variation with temperature for Nitrogen considered as a Lennard-Jones system with $\epsilon/k_B = 95.5$K and $\sigma = 3.76$ Angströms. Theoretical values in red and experimental values in small circles. Adapted from Cowan (27).</figcaption></figure>

Making a change of variable $x = r/\sigma$, we write:

$$B_2(T) = -2\pi\sigma^3 \int_0^\infty x^2[\exp\{-4(x^{-12} - x^{-6})\epsilon/k_B T\} - 1]dx \quad (85)$$

The integral can be obtained analytically (21) as:

$$B_2(T) = -2\pi\sigma^3 \sqrt{2\pi} \left(\frac{\epsilon}{k_B T}\right)^{1/4} H_{1/2} \left(-\sqrt{\frac{4\epsilon}{k_B T}}\right) \quad (86)$$

where $H_{1/2}$ is Hermite polynomial (21).

Fig. 9 displays $B_2(T)$ variation with temperature for a Lennard-Jones system as an illustration.

Let us consider a system governed by the Sutherland potential (see fig. 8), a simplified form of the Lennard-Jones given by:

$$\phi(r) = \infty, \quad r < \sigma, \quad \phi(r) = -\epsilon \left(\frac{\sigma}{r}\right)^6, \quad r > \sigma \quad (87)$$

In the Sutherland case, the second virial coefficient is obtained as (28):

$$B_2(T) = \frac{2\pi\sigma^3}{3} \left(1 - \frac{\epsilon}{k_B T}\right) = \frac{V_0}{2} \left(1 - \frac{\epsilon}{k_B T}\right) \quad (88)$$

where $V_0$ is an excluded volume originating from the presence of a hard-core potential accounting for a better estimation of gas volume. It is eight (2³) times the atomic volume ($\sigma$ the hard-core minimal approach distance being twice the atomic radius as inferred from eq. 82 and fig. 8).

The equation of state (to second order in density):

$$PV = Nk_B T \left(1 + \frac{V_0}{2} \left(1 - \frac{\epsilon}{k_B T}\right) n\right). \quad (89)$$

where $n = N/V$ can then be rearranged as:

$$\frac{1}{k_B T} \left( P + \frac{\epsilon V_0}{2} n^2 \right) = n \left( 1 + n \frac{V_0}{2} \dots \right) \approx \frac{n}{1 - n \frac{V_0}{2}} = \frac{N}{V - N \frac{V_0}{2}} \quad (90)$$
This can be rewritten as the Van der Waals equation of state:

$$ \left[ P + \frac{\epsilon V_0}{2} \left( \frac{N}{V} \right)^2 \right] \left[ V - N \frac{V_0}{2} \right] = N k_B T \quad (91) $$

## V. MOLECULAR SOLIDS

These solids go by at least eight different names: Rare gas, Noble Gas, Inert Gas, Lennard-Jones, Van Der Waals, 6-12, Molecular, Induced dipoles...

These solids are thus formed with Noble gases or molecules having closed shells. Bonding may result from a spatial charge fluctuation inducing an electric dipole moment.

The physical reason is the polarization of electron shells of the atoms and resulting dipole-dipole interaction which behaves at large distances as:

$$ \phi(r \to \infty) \propto -1/r^6. \quad (92) $$

The forces are of the same nature as those operating in real gases containing interacting particles in contrast with non-interacting particles making an ideal gas.

### A. Van der Waals interaction

In these crystals the interaction potential is described by the Lennard-Jones formula

$$ \phi(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \quad (93) $$

The equilibrium point where $(d\phi/dr)_{r^*} = 0$ yields $r^* = \sqrt[6]{2}\sigma = 1.12\sigma$. Thus bonding distance $r^* \sim \sigma$.

It is possible to perform this operation to an entire molecular crystal made of pairs of atoms/molecules interacting in a Lennard-Jones fashion. There are several algorithms to perform the corresponding lattice sums (2; 17) and evaluate the total energy $U$ from which any thermodynamic quantity may be evaluated.

In an FCC structure an atom/molecule is surrounded by 12 nearest-neighbors with a pair energy given by: $4\epsilon \left[ \left(\frac{\sigma}{R}\right)^{12} - \left(\frac{\sigma}{R}\right)^{6} \right]$ at a distance $R$. There are 6 next-nearest-neighbors with a pair energy given by:

$4\epsilon \left[ \left(\frac{\sigma}{R\sqrt{2}}\right)^{12} - \left(\frac{\sigma}{R\sqrt{2}}\right)^{6} \right]$ with a distance of $R\sqrt{2}$... The total energy is given numerically (2) by:

$$ U = \frac{1}{2} N 4\epsilon \left[ A \left(\frac{\sigma}{R}\right)^{12} - B \left(\frac{\sigma}{R}\right)^{6} \right] \quad (94) $$

where $N$ is the number of atoms/molecules and $A=12.13188$ and $B=14.45392$

Minimizing the total energy $(dU/dR)_{R^*} = 0$ gives the bonding distance $R^* = 1.09\sigma$ which is close to $\sqrt[8]{2}\sigma$.

<table><caption>**Table 4:** Values of $R^*/\sigma$ for several inert gas solids.</caption><thead><tr><th>Ne</th><th>Ar</th><th>Kr</th><th>Xe</th></tr></thead><tbody><tr><td>1.14</td><td>1.11</td><td>1.10</td><td>1.09</td></tr></tbody></table>

In the HCP case, the numerical values of $A, B$ are $A=12.13229$ and $B=14.45489$. (see Table. 8). The total energy is obtained at $R^*$ as: $-2.15(4N\epsilon)$.
<table><caption>**Table 5:** L-J potential parameters for the inert gas crystals. $\epsilon \sim$ eV/100 and $\sigma \sim$ Angströms. Helium is liquid at 0K, thus no melting point.</caption><thead><tr><th>Element</th><th>ε (10<sup>-21</sup> J)</th><th>σ (Angströms)</th><th>Melting point (K)</th><th>Ionization potential of free atom (eV)</th></tr></thead><tbody><tr><td>He</td><td>0.141</td><td>2.56</td><td>—</td><td>24.58</td></tr><tr><td>Ne</td><td>0.492</td><td>2.74</td><td>24</td><td>21.56</td></tr><tr><td>Ar</td><td>1.70</td><td>3.40</td><td>84</td><td>15.76</td></tr><tr><td>Kr</td><td>2.30</td><td>3.65</td><td>117</td><td>14.00</td></tr><tr><td>Xe</td><td>3.10</td><td>4.07</td><td>161</td><td>12.13</td></tr></tbody></table>

For illustration, let us evaluate the adiabatic bulk modulus (see Classical Thermodynamics section):

$$B_S = -V \left(\frac{\partial P}{\partial V}\right)_S = V \left(\frac{\partial^2 U}{\partial V^2}\right)_S \quad \text{since} \quad P = -\left(\frac{\partial U}{\partial V}\right)_S \qquad (95)$$

After performing the second derivation of the total energy $U$ we get:

$$B_S = \frac{\sqrt{2}[(14.45)N^3\epsilon\sigma^6]^{5/2}}{[\frac{1}{2}(12.13)N^5\epsilon\sigma^{12}]^{3/2}} \qquad (96)$$

where all factors are displayed explicitly for transparency.

# VI. IONIC SOLIDS

The physical reason is near complete transfer of the electron from the anion to the cation. It is important for the alkali crystals NaCl, KCl, CsCl, etc. One can consider the interaction as the Coulomb one for point charges at the lattice sites. Because the nearest neighboring ions have opposite sign with respect to the central one, Coulomb interaction is attractive.

## A. Madelung constant

Crystal bonding is based on short-range Slater repulsion and long-range Coulomb attraction:

$$\phi(r_{ij}) = \begin{cases} \lambda e^{-\alpha R} - \frac{e_*^2}{4\pi\epsilon_0 R} & \text{for nearest neighbors,} \\ \pm \frac{e_*^2}{4\pi\epsilon_0 r_{ij}} & \text{otherwise} \end{cases} \qquad (97)$$

with $r_{ij} = R_{p_{ij}}$ where $p_{ij}$ represent distances for the lattice sites, $R$ is nearest-neighbor distance and $e_*$ is the effective charge. So the total energy $U = N_p\phi(r_{ij})$ where $N_p$ the number of interacting pairs, becomes:

$$U = N_p \left( z\lambda e^{-\alpha R} - \alpha_M \frac{e_*^2}{4\pi\epsilon_0 R} \right)$$

where $z$ is the number of nearest neighbors (CN or coordination number) and $\alpha_M$ is Madelung constant.

In order to find the bonding distance, we take the derivative of $U$ with respect to $R$. Afterwards we proceed as in the Molecular crystal case to evaluate with lattice sums the thermodynamic energy $U$ allowing us to estimate any thermodynamic quantity.

Typical values of $\alpha_M$ for 3D lattices are given in Table 3.

Some authors define the Madelung constant normalized by the lattice parameter instead of the nearest-neighbor distance. Thus one must make sure about the convention used as emphasized by Tosi (19) in Table. 7.
<table><caption>**Table 6:** Madelung constants normalized by the nearest neighbor distance for several crystals. CN is coordination number. Note: CsCl Madelung constant is sometimes given as 1.76 and sometimes as 2.03. The difference comes from the inclusion of repulsive (non-electrostatic) terms! The Electrostatic Madelung constant of CsCl is 2.03.
*Exact value dependent on details of structure.</caption><thead><tr><th>Radius Ratio<br>(Cation/Anion)</th><th>Lattice Type</th><th>Cation CN</th><th>Anion CN</th><th>Madelung<br>Constant α<sub>M</sub></th></tr></thead><tbody><tr><td>0.225-0.414</td><td>Wurtzite and Zinc Blende ZnS</td><td>4</td><td>4</td><td>1.63805</td></tr><tr><td>0.414-0.732</td><td>Rock salt (NaCl)</td><td>6</td><td>6</td><td>1.74756</td></tr><tr><td>0.732-1.000</td><td>CsCl</td><td>8</td><td>8</td><td>1.76167</td></tr><tr><td>0.225-0.414</td><td>β-quartz (SiO<sub>2</sub>)</td><td>4</td><td>2</td><td>2.201</td></tr><tr><td>0.414-0.732</td><td>Rutile (Tetragonal TiO<sub>2</sub>)</td><td>6</td><td>3</td><td>2.408*</td></tr><tr><td>0.732-1.000</td><td>Fluorite (Cubic CaF<sub>2</sub>)</td><td>8</td><td>4</td><td>5.03878</td></tr><tr><td>0.414-0.732</td><td>Corundum (Rhombohedral Al<sub>2</sub>O<sub>3</sub>)</td><td>6</td><td>4</td><td>4.1719*</td></tr></tbody></table>

<table><caption>**Table 7:** The Madelung constants $α_M$ and $α_a$ refer respectively to the nearest-neighbor distance and to the lattice parameter $a$. Note the large discrepancy in the values.</caption><thead><tr><th>Madelung constant</th><th>CsCl</th><th>NaCl</th><th>Zinc Blende</th><th>CaF<sub>2</sub></th><th>Cu<sub>2</sub>O</th></tr></thead><tbody><tr><td>α<sub>M</sub></td><td>1.76267</td><td>1.74756</td><td>1.63805</td><td>5.03878</td><td>4.44248</td></tr><tr><td>α<sub>a</sub></td><td>2.03536</td><td>3.49513</td><td>3.78293</td><td>11.63657</td><td>10.25946</td></tr></tbody></table>

## VII. COVALENT SOLIDS

The nature of this directional bonding is purely quantum akin to bonding in the H₂ molecule where the H atoms share two electrons with anti-parallel spins obeying Pauli principle with a strong bonding energy of 4.5 eV.

In fact, for most semiconductors, bonding is mixed meaning it is partly ionic and partly covalent. The values of ionicity numbers or effective charge is given in Table 8.

Covalent bonding depends both on atomic orbital and on distance. It decreases exponentially with distance.

### A. Hybridization

Hybridization helps explain molecule geometry, since angles between bonds are approximately equal to the angles between hybrid orbitals.

We describe below, hybridization sp, sp² and sp³ types applied to Carbon (akin to Silicon, Ge, GaAs...):

1. sp

<table><caption>**Table 8:** Ionicity numbers for semiconductor crystals.</caption><thead><tr><th>Crystal</th><th>Ionicity</th></tr></thead><tbody><tr><td>Si</td><td>0.0</td></tr><tr><td>SiC</td><td>0.18</td></tr><tr><td>Ge</td><td>0.0</td></tr><tr><td>ZnSe</td><td>0.63</td></tr><tr><td>ZnS</td><td>0.62</td></tr><tr><td>CdSe</td><td>0.70</td></tr><tr><td>InP</td><td>0.42</td></tr><tr><td>InAs</td><td>0.46</td></tr><tr><td>InSb</td><td>0.32</td></tr><tr><td>GaAs</td><td>0.31</td></tr><tr><td>GaSb</td><td>0.36</td></tr></tbody></table>
In this type, the 2s orbital is mixed with only one of the three p orbitals, $p_x, p_y, p_z$ resulting in two sp and two p orbitals. The chemical bonding in acetylene ($C_2H_2$) consists of sp-sp overlap between two C atoms forming a $\sigma$ bond and two additional $\pi$ bonds formed by p-p overlap. Each C also bonds to hydrogen in a $\sigma$ s-sp overlap at 180° angles.

## 2. sp²

In sp² hybridization the 2s orbital is mixed with only two of three available 2p orbitals $2p_x$ and $2p_y$. The third 2p orbital $2p_z$ is unhybridized, forming a total of three sp² orbitals with one remaining p orbital.

This is the $C_2H_4$ molecule case where C is sp² hybridized, because one $\pi$ bond is required for the double bond between C while only three $\sigma$ bonds are formed per C atom.

Quantum mechanically, the lowest energy is obtained if all three bonds are equivalent. That can be obtained with linear combinations of the valence-shell s and p wave functions, which are three sp² hybrids.

This may be viewed geometrically as follows. It suffices to take combinations of an s orbital with two p-orbitals $|p_x\rangle$ and $|p_y\rangle$. Since $|p_z\rangle$ is not used, we get a planar $x-y$ structure:

$$ 
\begin{aligned}
|\psi_1\rangle &= \frac{1}{\sqrt{3}} |s\rangle + \sqrt{\frac{2}{3}} |p_x\rangle \\
|\psi_2\rangle &= \frac{1}{\sqrt{3}} |s\rangle - \frac{1}{\sqrt{6}} |p_x\rangle + \frac{1}{\sqrt{2}} |p_y\rangle \\
|\psi_3\rangle &= \frac{1}{\sqrt{3}} |s\rangle - \frac{1}{\sqrt{6}} |p_x\rangle - \frac{1}{\sqrt{2}} |p_y\rangle
\end{aligned}
\qquad (98) $$

Normality and completeness are checked with summing the square of the coefficients horizontally and vertically: for instance, the first row gives: $(\frac{1}{\sqrt{3}})^2 + (\sqrt{\frac{2}{3}})^2 = 1$ whereas the first column yields: $(\frac{1}{\sqrt{3}})^2 + (\frac{1}{\sqrt{3}})^2 + (\frac{1}{\sqrt{3}})^2 = 1$.

<figure><img src="image_11.png" /><figcaption>Fig.10: (Color on-line) 2D trigonal structure with the 120° angle.</figcaption></figure>

Noting that the combination of $|p_x\rangle, |p_y\rangle$ orbitals can be rewritten as:

$$ 
\begin{aligned}
-\frac{1}{\sqrt{6}} (|p_x\rangle - \sqrt{3} |p_y\rangle) \\
-\frac{1}{\sqrt{6}} (|p_x\rangle + \sqrt{3} |p_y\rangle)
\end{aligned}
\qquad (99) $$
we deduce that

$$ \begin{aligned} \frac{1}{\sqrt{6}}(-|p_x\rangle + \sqrt{3}|p_y\rangle) &\sim [\bar{1}10] \\ \frac{1}{\sqrt{6}}(-|p_x\rangle - \sqrt{3}|p_y\rangle) &\sim [\bar{1}\bar{1}0] \end{aligned} \qquad (100) $$

where $[\bar{1}10], [\bar{1}\bar{1}0]$ are the equivalent Miller directions (see fig. 10).

Geometrically this is equivalent in turn to:

$$ \begin{aligned} \boldsymbol{x} &\sim [100] \\ \frac{1}{2}(-\boldsymbol{x} + \sqrt{3}\boldsymbol{y}) &\sim [\bar{1}10] \\ \frac{1}{2}(-\boldsymbol{x} - \sqrt{3}\boldsymbol{y}) &\sim [\bar{1}\bar{1}0] \end{aligned} \qquad (101) $$

where $(\boldsymbol{x}, \boldsymbol{y})$ are unit vectors along $x, y$ directions (see fig. 10).

The scalar product between each pair of the above normalized vectors gives $-1/2$ whose inverse cosine is $120^\circ$ as visible in fig. 10.

### 3. sp³

For a tetrahedrally coordinated C (like methane CH₄), the C atom should have four orbitals directed towards the four surrounding hydrogen atoms.

Four sp³ hybrid orbitals are overlapped by hydrogen 1s orbitals, yielding four σ bonds (that is, four single covalent bonds) of equal length and strength.

Let us show below how one can view tetrahedral bonding based on (sp³) hybridization.

Quantum mechanically, the lowest energy is obtained if the four bonds are equivalent, which requires that they are formed from equivalent orbitals in Silicon. A set of four equivalent orbitals can be obtained that are linear combinations of the valence-shell (core orbitals are almost never involved in bonding) s and p wave functions, which are the four sp³ hybrids.

This may be seen geometrically as follows. We consider the combinations of an s orbital with 3 p-orbitals. Thus we get:

$$ \begin{aligned} |\psi_1\rangle &= \frac{1}{2}(|s\rangle + |p_x\rangle + |p_y\rangle + |p_z\rangle) \sim [111] \\ |\psi_2\rangle &= \frac{1}{2}(|s\rangle + |p_x\rangle + |p_y\rangle + |p_z\rangle) \sim [\bar{1}11] \\ |\psi_3\rangle &= \frac{1}{2}(|s\rangle + |p_x\rangle + |p_y\rangle + |p_z\rangle) \sim [\bar{1}\bar{1}1] \\ |\psi_4\rangle &= \frac{1}{2}(|s\rangle + |p_x\rangle + |p_y\rangle + |p_z\rangle) \sim [11\bar{1}] \end{aligned} \qquad (102) $$

This construction satisfies the normalization and completeness criteria as explained in the sp² case.
<figure><img src="image_12.png" /><figcaption>**Fig.11:** (Color on-line) At left: Diamond structure displaying four tetrahedra inscribed inside the cube with their centers in blue. At right: Tetrahedron structure showing the different Miller directions.</figcaption></figure>

<table><thead><tr><th>Coordination number</th><th>Geometry</th><th>Hybridization</th><th>Bond Angle</th><th>Example</th></tr></thead><tbody><tr><td>2</td><td>1D Linear</td><td>sp hybridization</td><td>180°</td><td>CO<sub>2</sub></td></tr><tr><td>3</td><td>2D Trigonal</td><td>sp<sup>2</sup> hybridization</td><td>120°</td><td>BCl<sub>3</sub></td></tr><tr><td>4</td><td>3D Tetrahedral</td><td>sp<sup>3</sup> hybridization</td><td>109.5°</td><td>CCl<sub>4</sub></td></tr></tbody></table>

**Table 9:** In sp hybridization, bond angle is 180° due to σ bonds, whereas in sp² hybridization, the s orbital is mixed with only two of the three available p orbitals, thus bond angle is 120°. In sp³ bond angle $\theta = \cos^{-1}(-\frac{1}{3}) = 109.5°$ since we have the scalar product of the tetrahedral directions yield -1/3 as in fig. 11.

where [111], [111]... are the equivalent Miller directions (see fig. 11).

Geometrically this is equivalent in turn to:

$$ 
\begin{aligned}
\frac{1}{\sqrt{3}}(\mathbf{x} + \mathbf{y} + \mathbf{z}) &\sim [111] \\
\frac{1}{\sqrt{3}}(-\mathbf{x} + \mathbf{y} + \mathbf{z}) &\sim [\bar{1}11] \\
\frac{1}{\sqrt{3}}(\mathbf{x} - \mathbf{y} + \mathbf{z}) &\sim [1\bar{1}1] \\
\frac{1}{\sqrt{3}}(\mathbf{x} + \mathbf{y} - \mathbf{z}) &\sim [11\bar{1}]
\end{aligned}
\qquad (103) $$

The scalar product between each pair of the above normalized vectors gives -1/3 whose inverse cosine is 109.5° as visible in fig. 11.

The three hybridization modes are summarized in Table 9.

## VIII. METALS

Metals usually form closed packed FCC, BCC, or HCP structures where electrons are shared by all the atoms (Fermi sea). The bonding energy is determined by a balance between the negative energy of Coulomb interaction between electrons and positive ions (energy proportional to $e^2/a$) and positive kinetic energy of electron Fermi gas ($\propto n^{2/3} \propto 1/a^2$) where a is the lattice parameter and n the electron density in the Fermi sea.
References

[1] O. Madelung, *Introduction to solid state theory*, Springer-Verlag (1996).

[2] C. Kittel, *Introduction to Solid State Physics*, 8th edition, Wiley, New-York (2005).

[3] J. J. Gilman, J. App. Phys. **31**, 2308 (1960).

[4] P. M. Morse, Phys. Rev. **34**, 57 (1930).

[5] F.H. Stillinger and T. A. Weber, Phys. Rev. B **31**, 5262 (1985).

[6] L. M. Woods and G. D. Mahan Phys. Rev. B **61**, 10651 (2000).

[7] J. Tersoff, New empirical model for the structural properties of silicon, Phys. Rev. Lett. **56**, 632 (1986).

[8] D. W. Brenner, Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Phys. Rev. **B** **42**, 9458 (1990).

[9] A. Zunger, Nat. Rev. Chem. **2**, 0121 (2018).

[10] C. L. Kane, E. J. Mele, Phys. Rev. Lett. **95**, 226801 (2005).

[11] H. Bilz, and W. Kress, *Phonon Dispersion Relations in Insulators*, Springer-Verlag, Berlin (1979).

[12] T. Aizawa, R. Souda, S. Otani, Y. Ishizawa and C. Oshima, Phys. Rev. B **42**, 11469 (1990).

[13] V. K. Tewary and B. Yang, Phys. Rev. B **79**, 075442 (2009).

[14] M. Mohr, J. Maultzsch, E. Dobardžić, I. Milošević, M. Damnjanović, A. Bosak, M. Krish, and C. Thomsen, Phys. Rev. B **76**, 035439 (2007).

[15] K. H. Michel and B. Verberck, Phys. Rev. B **78**, 085424 (2008).

[16] D. A. McQuarrie *Statistical Mechanics*, Harper and Row, New York (1973).

[17] J. M. Borwein, M. L. Glasser, R. C. McPhedran, J. G. Wan and I. J. Zucker, Encyclopedia of mathematics and its applications, *Lattice Sums: Then and Now*, Cambridge University Press (2013).

[18] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, *NIST Handbook of Mathematical Functions* Cambridge University Press, New York, (2010).

[19] M. Tosi, Solid State Physics Vol. 16, F. Seitz and D. Turnbull editors, Academic Press (1965).

[20] F.J. Vesely, *Computational Physics: An Introduction*, Kluwer-Plenum (2001).

[21] I. S. Gradstein and I. M. Ryzhik, *Table of Integrals, Series and Products*, Academic Press, New-York (1980).

[22] A. Burrows, S. Cooper and P. Schwerdtfeger, *The Madelung constant in N dimensions*. Proc. R. Soc. A 478: 20220334 (2022).

[23] F. Reif, *Statistical and Thermal Physics*, McGraw-Hill, New-York (1985).

[24] H.B. Callen, *Thermodynamics and an introduction to thermostatistics*, 2nd edition, Wiley, New-York (1985).

[25] L. D. Landau and E. M. Lifshitz, *Statistical Physics, Part1*, Third Edition, Vol. 5, Pergamon, Oxford (1976).

[26] A. L. Fetter and J. D. Walecka, *Quantum Theory of Many-Particle Systems*, McGraw-Hill, New York (1971).

[27] B. Cowan, *Topics in Statistical Mechanics*, 2nd edition, World Scientific, Singapore (2021).

[28] M. Kardar, *Statistical Physics of Particles*, Cambridge University Press (2007).