# Predicting the mechanical properties of spring networks

Doron Grossman<EMAIL_ADDRESS>LadHyX, CNRS, Ecole
polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau Cedex, France
Arezki Boudaoud  LadHyX, CNRS, Ecole polytechnique, Institut Polytechnique de
Paris, 91128 Palaiseau Cedex, France

###### Abstract

The elastic response of mechanical, chemical, and biological systems is often
modeled using a discrete arrangement of Hookean springs, either modeling
finite material elements or even the molecular bonds of a system. However, to
date, there is no direct derivation of the relation between discrete spring
network, and a general elastic continuum. Furthermore, understanding the
networks’ mechanical response requires simulations that may be expensive
computationally. Here we report a method to derive the exact elastic continuum
model of any discrete network of springs, requiring network geometry and
topology only. We identify and calculate the so-called ”non-affine”
displacements. Explicit comparison of our calculations to simulations of
different crystalline and disordered configurations, shows we successfully
capture the mechanics even of auxetic materials. Our method is valid for
residually stressed systems with non-trivial geometries, is easily
generalizable to other discrete models, and opens the possibility of a
rational design of elastic systems.

Ever since its introduction in the 19th century [1, 2], the theory of
elasticity has been phenomenological. That is - to date, it has never been
derived from first principles as continuum limit, and the elastic properties
of a material whose microscopic characteristics are known, could not be
computed, in general. A notable exception to this statement is the continuum
limit of crystals [3, 4]. However, by definition, these are essentially flat,
ordered, systems i.e - a very limited (though common) class of materials.
Despite this missing link between the microscopic scale and the macroscopic,
it is widely accepted that in essence, linear elasticity may be described
using spring-like interactions between constituents (e.g. first order
approximation of intermolecular forces around equilibrium). In fact,
elasticity is often described using a spring network either for computational
or analytical aspects [5, 6] in a plethora of different systems and cases -
from modeling the shape of self assembled membranes [7, 8, 9, 10], through
biological systems [11, 12], to modeling crack propagation [13] and various
bio-inspired, and meta materials [14, 15, 16]. Typically, calculation of the
network’s elastic response can only be done via direct simulation of a loading
scheme (i.e. simulating a mechanical load and the response to it).

In this paper we directly derive a generalized elastic continuum limit of any
triangulated spring network, with arbitrary reference lengths and spring
constants, in two and three dimensions, thus solving an age-long question. The
resulting continuum limit depends only on the network geometry and topology,
as expressed by reference lengths, spring constants, and bonds, not on
anything else. From this description, any macroscopic elastic quantity can be
extracted, such as Poisson’s ratio. We then demonstrate the strength of the
approach, by calculating Poisson’s ratios for different ordered, and
disordered, systems, recovering even negative Poisson’s ratio values in
isotropic, disordered, foam - like configurations. We identify the so-called
”non-affine” displacements that are local deformations which deviate from the
local average deformation, and are responsible of the wide range of responses
seen in disordered elastic media. The results are valid for residually
stressed elastic systems as well.

The continuum limit we derive is formulated within the theory of incompatible
elasticity,[17] which is a modern formulation of elasticity, that successfully
describes residually stressed elastic systems [18, 19, 20]. It is a
generalization of standard theory of elasticity. In this formulation, an
elastic material is described by a metric $\mathbf{g}$ with elements
$g_{\mu\nu}$, which describes actual distances between material elements, and
a reference metric $\mathbf{\bar{g}}$ with elements $\bar{g}_{\mu\nu}$
describing ideal distances. The elastic energy then depends on the squared
difference of $\mathbf{g}-\mathbf{\bar{g}}$,
$E_{el}\propto\left\|\bf{g-\bar{g}}\right\|^{2}$, for some proper choice (yet
to be defined) of the norm $\|\cdot\|^{2}$, through the elastic (four-indexed)
tensor $\bar{\bar{A}}$ (with elements $A^{\mu\nu\alpha\beta}$).

This description, via use of metrics, is very similar in essence and form to
the classical description of Hookean springs. It is independent from the
assumption about the existence of a rest configuration. It is this
independence that enables the treatment of residual stresses. In the
following, we will consider a discrete network of springs and show how such
formulation naturally arises. We will then coarse grain the network and will
identify the non-affine quantities, and show how they contribute to the
elastic continuum energy.

##  I Results

###  I.1 Framework

The theory of incompatible elasticity [17] is the framework to which the
results of this paper are anchored to. Within it, the elastic energy is given
by:

$\displaystyle E_{el}=\int
A^{\mu\nu\alpha\beta}\left(g_{\mu\nu}-\bar{g}_{\mu\nu}\right)\left(g_{\alpha\beta}-\bar{g}_{\alpha\beta}\right)\,\text{d}V_{\bar{g}}$
(1)

Where $g_{\mu\nu}$ is the actual metric, describing distances between
neighboring material points, $\bar{g}_{\mu\nu}$ is the reference metric,
describing ideal distances. $\,\text{d}V_{\bar{g}}=\sqrt{\bar{g}}d^{D}x$ is
the volume element in $D$ dimensions, $\bar{g}=\det{\bf\bar{g}}$.
$A^{\mu\nu\alpha\beta}$ is the elastic tensor. Einstein’ summation is assumed
for repeating upper and lower Greek indices. Greek indices refer to
coordinated within the surface/volume of the $D$ dimensional manifold.

In an isotropic material,
$A^{\mu\nu\alpha\beta}=\frac{Y}{16(1-\nu^{2})}\left[\frac{1}{2}(1-\nu)\left(\bar{g}^{\mu\alpha}\bar{g}^{\nu\beta}+\bar{g}^{\nu\alpha}\bar{g}^{\mu\beta}\right)+\nu\bar{g}^{\mu\nu}\bar{g}^{\alpha\beta}\right]$,
where $\bar{g}^{\mu\nu}$ is the inverse reference metric, $Y$ is Young’s
modulus, setting the rigidity scale of the system, and $\nu$ is Poisson’s
ratio, describing the amount a material contracts in one axis, when the other
is stretched (negative values indicate expansion). In a non isotropic
materials Poisson’s ratio is orientation dependent, and the expression of
$A^{\mu\nu\alpha\beta}$ will typically depend on additional terms.

The elastic stress is given by the variation $\sigma^{\mu\nu}=\frac{\delta
E_{el}}{\delta g_{\mu\nu}}$, and the material satisfies the usual force
balance equation:

$\displaystyle\bar{\nabla}_{\mu}\sigma^{\mu\nu}+\sigma^{\mu\alpha}\left(\Gamma_{\mu\alpha}^{\nu}-\bar{\Gamma}_{\mu\alpha}^{\nu}\right)=f_{ext}^{\nu}$
(2)

where $f_{ext}^{\nu}$ are the external forces acting on the systems, and
$\sigma^{\mu\nu}=\frac{\delta E_{el}}{\delta g_{\mu\nu}}$ is the elastic
stress, $\bar{\nabla}_{\mu}$ is the covariant derivative with respect to
$\bar{g}_{\mu\nu}$ and
$\Gamma^{\alpha}_{\beta\gamma},\bar{\Gamma}^{\alpha}_{\beta\gamma}$ are the
christoffel symbols of $g_{\mu\nu}$ and $\bar{g}_{\mu\nu}$ respectively.

###  I.2 Analytical Derivation

We begin by considering a triangulated mesh of springs, each with reference
length $\ell_{e}$, spring constant $k_{e}$ and an actual length $l_{e}$, where
the index $e$ enumerates the springs. The elastic energy of the systems is
exactly given by

$\displaystyle
E_{el}=\sum_{e}\frac{1}{2}k_{e}\left(l_{e}-\ell_{e}\right)^{2}.$ (3)

Since the network is triangulated we may easily divide the sum unto specific
simplexes (or cells). In three dimensions these are tetrahedrons, and in two
dimensions these are simple triangles

$\displaystyle E_{el}=\frac{1}{2}\sum_{s}\sum_{e\in
s}\frac{1}{2}k_{e}\left(l_{e}-\ell_{e}\right)^{2}.$ (4)

Here, the index $s$ enumerates simplexes, the marking $e\in s$ means summation
over all the springs associated with the simplex $s$. Regardless of parameter
values, this network, when left to relax, assumes some configuration (not
necessarily unique) $\\{\vec{f}_{v}\\}$ in $\mathbb{R}^{n}$ for every vertex.
Setting coordinates $x^{\mu}_{v}$ to each vertex $v$ ($\mu$ is the coordinate
component), and given actual lengths $\\{l_{e}\\}$, we may uniquely define a
”local metric” $g^{(s)}_{\mu\nu}$ associated with a cell $s$, so that

$\displaystyle l_{e}^{2}=g^{(s)}_{\mu\nu}\Delta x^{\mu}_{e}\Delta
x^{\nu}_{e}~{}~{}~{}\forall e\in s.$ (5)

Where $\Delta x^{\mu}_{(1,2)}=x_{2}^{\mu}-x_{1}^{\mu}$ is the ”coordinate
difference” of the edge $e=(1,2)$ connecting vertexes $1$ and $2$. Note that
(5) is not an approximation, rather it is an exact definition of the local
quantity $g^{(s)}_{\mu\nu}$, over the whole simplex.

A given simplex, therefore, defines uniquely a local metric. Thus, every cell,
$s$, has an associated actual metric $g^{(s)}_{\mu\nu}$. For a given set of
coordinates, a physical systems is constrained such that any two local metrics
$\mathbf{g^{(i)}}$ and $\mathbf{g^{(j)}}$ with a shared edge $\Delta x_{e}$,
agree on its length- $l_{e}[g^{(i)}]=l_{e}[g^{(j)}]$, where $l_{e}[g]$ is the
edge’s length, as measured using the metric $\mathbf{g}$. We can now rewrite
the energy -

$\displaystyle E_{el}=\frac{1}{4}\sum_{s}\sum_{e\in
s}k_{e}\left(\sqrt{{g^{(s)}}_{\mu\nu}\Delta x_{e}^{\mu}\Delta
x_{e}^{\nu}}-\ell_{e}\right)^{2}.$ (6)

In order to advance, the following assumptions are introduced: First we assume
that the reference lengths $\\{\ell_{e}\\}$ are compatible in the sense that
in the absence of other parts of the network, a simplex can assume the shape
described by the lengths $\ell_{e}$. This means that the reference lengths
locally define a reference metric ${\bar{g}^{(s)}}_{\mu\nu}$. Under these
assumptions we can write $\ell_{e}=\sqrt{\bar{g}^{(s)}_{\mu\nu}\Delta
x_{e}^{\mu}\Delta x_{e}^{\nu}}$.

Second, we assume that in an equilibrium configuration (again, not necessarily
unique), deviations of actual lengths from the reference lengths are small.

$\displaystyle
l_{e}-\ell_{e}=\frac{l_{e}^{2}-\ell_{e}^{2}}{l_{e}+\ell_{e}}\simeq\frac{1}{2\ell_{e}}\left(l_{e}^{2}-\ell_{e}^{2}\right)+\cdots$
(7)

The right hand side can be further expanded in higher terms of
$l_{e}^{2}-\ell_{e}^{2}$ if required. Thus

$\displaystyle E_{el}=\sum_{s}\sum_{e\in
s}\frac{k_{e}}{16\ell_{e}^{2}}\left(g^{(s)}_{\mu\nu}\Delta x_{e}^{\mu}\Delta
x_{e}^{\nu}-\bar{g}^{(s)}_{\mu\nu}\Delta x_{e}^{\mu}\Delta
x_{e}^{\nu}\right)^{2}+\cdots$ (8)

Keeping to leading order, we expand the energy:

$\displaystyle E_{el}=$
$\displaystyle\sum_{s}\left(g^{(s)}_{\mu\nu}-\bar{g}^{(s)}_{\mu\nu}\right)\left(g^{(s)}_{\alpha\beta}-\bar{g}^{(s)}_{\alpha\beta}\right)\sum_{e\in
s}\frac{k_{e}\Delta x_{e}^{\mu}\Delta x_{e}^{\nu}\Delta x_{e}^{\alpha}\Delta
x_{e}^{\beta}}{16\ell_{e}^{2}}.$ (9)

Marking the local elastic tensor

$A_{(s)}^{\mu\nu\alpha\beta}=\sum_{e\in s}\frac{k_{e}\Delta x_{e}^{\mu}\Delta
x_{e}^{\nu}\Delta x_{e}^{\alpha}\Delta x_{e}^{\beta}}{16\ell_{e}^{2}},$ (10)

we write -

$\displaystyle E_{el}=$
$\displaystyle\sum_{s}A_{(s)}^{\mu\nu\alpha\beta}\left(g^{(s)}_{\mu\nu}-\bar{g}^{(s)}_{\mu\nu}\right)\left(g^{(s)}_{\alpha\beta}-\bar{g}^{(s)}_{\alpha\beta}\right)$
(11)

This equation is very similar to (1), and may be considered as a discrete
version of that equation.

Finally, a third assumption in introduced - that $\bar{g}_{s}$ varies slowly
on some large enough region. Without this assumption the continuum limit
cannot hold (though an effective continuum may be derived, in principle).
Defining an average metric, ${g}_{\mu\nu}$, on some (possibly small)
neighborhood, we may expand $g^{(s)}_{\mu\nu}={g}_{\mu\nu}+\delta
g^{(s)}_{\mu\nu}$. Formally, $\delta g^{(s)}_{\mu\nu}$ describe the ”non-
affine” deformations. The energy then reads -

$\displaystyle
E_{el}=\sum_{\Omega}\left(\sum_{s\in\Omega}A_{(s)}^{\mu\nu\alpha\beta}\Delta
g_{\mu\nu}\Delta
g_{\alpha\beta}+\sum_{s\in\Omega}A_{(s)}^{\mu\nu\alpha\beta}\Delta
g_{\mu\nu}\delta
g^{(s)}_{\alpha\beta}+\sum_{s\in\Omega}A_{(s)}^{\mu\nu\alpha\beta}\delta
g_{\mu\nu}\Delta
g^{(s)}_{\alpha\beta}+\sum_{s\in\Omega}A_{(s)}^{\mu\nu\alpha\beta}\delta
g^{(s)}_{\mu\nu}\delta g^{(s)}_{\alpha\beta}\right),$ (12)

where $\Delta g_{\mu\nu}={g}_{\mu\nu}-\bar{{g}}_{\mu\nu}$, and $\sum_{\Omega}$
is the sum over all the neighborhoods in which $\mathbf{g}$ and
$\mathbf{\bar{g}}$ may be regarded as constant. Since, under our assumptions,
if ${g}_{\mu\nu}={\bar{g}}_{\mu\nu}$, $\delta g^{(s)}_{\mu\nu}=0\,\forall s$,
then for small deviations from $\bar{{g}}_{\mu\nu}$, $\delta
g^{(s)}_{\mu\nu}=W_{(s)\,\mu\nu}^{\alpha\beta}\Delta g_{\alpha\beta}$. While
mathematically a different entity, we identify the proportion tensors
$W_{(s)\,\mu\nu}^{\alpha\beta}$, with ”non affine” deformations of each
simplex, which are yet unknown.

The elastic energy $E_{el}^{\Omega}$ within a single neighborhood
($E_{el}=\sum_{\Omega}E_{el}^{\Omega}$) then reads

$\displaystyle E_{el}^{\Omega}=$
$\displaystyle\sum_{s}\left(A_{(s)}^{\mu\nu\alpha\beta}+A_{(s)}^{\mu\nu\rho\sigma}W_{(s)\,\rho\sigma}^{\alpha\beta}+A_{(s)}^{\alpha\beta\rho\sigma}W_{(s)\,\rho\sigma}^{\mu\nu}+A_{(s)}^{\tau\lambda\rho\sigma}W_{(s)\,\rho\sigma}^{\alpha\beta}W_{(s)\,\tau\lambda}^{\mu\nu}\right)\Delta
g_{\mu\nu}\Delta g_{\alpha\beta}$ (13)
$\displaystyle+\chi^{\alpha\beta}\sum_{s}W_{(s)\alpha\beta}^{\mu\nu}\Delta
g_{\mu\nu}$

The second line is a Lagrange term forcing the requirement that
$\sum_{s}\delta g_{(s)\mu\nu}=0$. We note, that as $g_{\mu\nu}$ is an average
metric, $\sum_{s}\delta g^{(s)}_{\mu\nu}=0$. Therefore our requirement
translates to $\sum_{s}W_{(s)\,\mu\nu}^{\alpha\beta}=0$. By defining $n$ as
the number of simplexes in the neighborhood $\Omega$, and the averages

$\displaystyle A^{\mu\nu\alpha\beta}$
$\displaystyle=\frac{1}{n}\sum_{s}A^{\mu\nu\alpha\beta}_{(s)}$ (14)
$\displaystyle\tilde{A}^{\mu\nu\alpha\beta}$
$\displaystyle=\frac{1}{n}\sum_{s}\left(A_{(s)}^{\mu\nu\alpha\beta}+A_{(s)}^{\mu\nu\rho\sigma}W_{(s)\,\rho\sigma}^{\alpha\beta}+A_{(s)}^{\alpha\beta\rho\sigma}W_{(s)\,\rho\sigma}^{\mu\nu}+A_{(s)}^{\tau\lambda\rho\sigma}W_{(s)\,\rho\sigma}^{\alpha\beta}W_{(s)\,\tau\lambda}^{\mu\nu}\right)$
(15) $\displaystyle\delta A^{\mu\nu\alpha\beta}_{(s)}$
$\displaystyle={A}^{\mu\nu\alpha\beta}_{(s)}-A^{\mu\nu\alpha\beta},$ (16)

we may find the $W_{(s)\,\mu\nu}^{\alpha\beta}$ for a finite (small) strain by
solving the coupled set of equations

$\displaystyle\bar{\nabla}_{\mu}\sigma^{\mu\nu}+\sigma^{\mu\alpha}\left(\Gamma_{\mu\alpha}^{\nu}-\bar{\Gamma}_{\mu\alpha}^{\nu}\right)=f_{ext}^{\nu}$
(17)
$\displaystyle\left(A_{(s)}^{\mu\nu\alpha\beta}+A_{(s)}^{\lambda\tau\alpha\beta}W_{(s)\,\lambda\tau}^{\mu\nu}\right)\Delta
g_{\mu\nu}+\frac{1}{2}\chi^{\alpha\beta}=0$

The first line is actually the elastic equation(2). However, at this point is
enough to note that $\sigma^{\mu\nu}=\frac{\delta E_{el}}{\delta g_{\mu\nu}}$,
where $E_{el}$ is given by eq. (13), and that $\Gamma^{\alpha}_{\beta\gamma}$
and $\bar{\Gamma}^{\alpha}_{\beta\gamma}$ are the Christoffel symbols
associates with the (now well defined) metrics $\bf{g}$ and $\bf{\bar{g}}$,
and $\bar{\nabla}_{\mu}$ is the connection (covariant derivative) associated
with the metric $\bf{\bar{g}}$. By marking $\langle
A_{(s)}W_{(s)}\rangle^{\mu\nu\alpha\beta}=\frac{1}{n}\sum_{s}A_{(s)}^{\lambda\tau\alpha\beta}W_{(s)\,\lambda\tau}^{\mu\nu}$,
and using the fact that $\langle
A_{(s)}W_{(s)}\rangle^{\mu\nu\alpha\beta}=\langle\delta
A_{(s)}W_{(s)}\rangle^{\mu\nu\alpha\beta}$, we may rewrite the second line of
(17), after a little algebra

$\displaystyle\delta
A_{(s)}^{\mu\nu\alpha\beta}+A_{(s)}^{\lambda\tau\alpha\beta}W_{(s)\,\lambda\tau}^{\mu\nu}-\langle\delta
A_{(s)}W_{(s)}\rangle^{\mu\nu\alpha\beta}=0.$ (18)

This linear equation for the non-affine deformation terms, $W$, depends solely
on the network’s local topology, as given by the different
$A_{(s)}^{\mu\nu\alpha\beta}$ values. It is solved by mapping the tensor
components and indexes unto a multi - index notation and using the symmetries
of the tensors (i.e in two dimensions $W_{(s)\,\mu\nu}^{\alpha\beta}$ has only
9 independent entries, while in three dimension it has 36)

$\displaystyle\delta
A_{S}+\sum_{S^{\prime}}\left(A_{SS^{\prime}}-B_{SS^{\prime}}\right)W_{S^{\prime}}=0$
(19)

Where $\delta A_{S}$, $A_{SS^{\prime}}$ and $B_{SS^{\prime}}$ are merely
reorganizations of the elements of $\\{A_{s}^{\mu\nu\alpha\beta}\\}$ and
$\\{\delta A_{s}^{\mu\nu\alpha\beta}\\}$ into matrices compatible with the new
multi index (see appendix A). The solution is now simple -

$\displaystyle
W_{S}=-\sum_{S^{\prime}}{\left[A-B\right]^{-1}}_{SS^{\prime}}\delta
A_{S^{\prime}}.$ (20)

We may now write the elastic energy -

$\displaystyle E_{el}=$
$\displaystyle\sum_{\Omega}n\tilde{A}^{\mu\nu\alpha\beta}\left(g_{\mu\nu}-\bar{g}_{\mu\nu}\right)\left(g_{\alpha\beta}-\bar{g}_{\alpha\beta}\right).$
(21)

This equation is essentially already coarse grained, it is easier to see it by
adding a unit term inside the sum -

$\displaystyle E_{el}=$
$\displaystyle\sum_{\Omega}n\tilde{A}^{\mu\nu\alpha\beta}\left(g_{\mu\nu}-\bar{g}_{\mu\nu}\right)\left(g_{\alpha\beta}-\bar{g}_{\alpha\beta}\right)\frac{V^{\Omega}_{\mathbf{\bar{g}}}}{V^{\Omega}_{\mathbf{\bar{g}}}}=\sum_{\Omega}\rho_{\Omega}\tilde{A}^{\mu\nu\alpha\beta}\left(g_{\mu\nu}-\bar{g}_{\mu\nu}\right)\left(g_{\alpha\beta}-\bar{g}_{\alpha\beta}\right)\int_{\Omega}\,\text{d}V_{\mathbf{\bar{g}}}$
(22) $\displaystyle=$
$\displaystyle\int\tilde{A}^{\mu\nu\alpha\beta}\left(g_{\mu\nu}-\bar{g}_{\mu\nu}\right)\left(g_{\alpha\beta}-\bar{g}_{\alpha\beta}\right)\,\text{d}V_{\bar{g}}$

where $V^{\Omega}_{\mathbf{\bar{g}}}$ is the volume of the neighborhood
$\Omega$, $\rho_{\Omega}=n/V^{\Omega}_{\mathbf{\bar{g}}}$ is the local density
(which we absorb into the definition of
$\tilde{A}^{\mu\nu\alpha\beta}$),$\int_{\Omega}\,\text{d}V_{\mathbf{\bar{g}}}$
is an integral over the region $\Omega$ and we use the fact the
$\sum_{\Omega}\int_{\Omega}\,\text{d}V_{\mathbf{\bar{g}}}=\int\,\text{d}V_{\mathbf{\bar{g}}}$
over all the network.

Eqs. (22), and (20) form the central result of this work. Together with the
definitions of
$A^{\mu\nu\alpha\beta},\tilde{A}_{s}^{\mu\nu\alpha\beta},W_{(s)\,\lambda\tau}^{\mu\nu}$,
they fully describe the response of the network and offer a novel way of
computing it directly from the network geometry, without the need to consider
any specific load. Under this view the metrics $\mathbf{g}$ is the the actual,
coarse grained, metric of the system, $\mathbf{\bar{g}}$ describes the
reference geometry of the system, and $\tilde{A}^{\mu\nu\alpha\beta}$ is the
coarse grained elastic tensor, governing the mechanical response (as opposed
to the local or ”bare” term $A^{\mu\nu\alpha\beta}_{(s)}$).
$W^{\mu\nu}_{(s)\,\alpha\beta}$ describe the non-affine displacements.

###  I.3 Comparison to simulation

Results were tested numerically by comparing the expected Poisson’s ratio
using the above scheme, to that of simulated two dimensional triangulated
spring networks. In general, we find a very good agreement between theory and
simulation.

We used our assumptions that $\bf\bar{g}$ is well defined on a large enough
region to limit our numerical and analytical solution for the case
$\bar{g}_{\mu\nu}=\delta_{\mu\nu}$, as we can always work in a locally flat
frame. This condition is sufficient as we want to isolate the effects of the
non-trivial structure of the network itself, not the whole (non uniform)
mechanical response of a complex, possibly residually stressed structure. We
compared the results of 3 test cases, ordered (non isotropic),foam like
(following procedure in [15]), and a honeycomb. The latter can be calculated
analytically, rather than simulated. Exact simulation protocols and details,
can be found in the methods section.

In parallel to simulation, for every network architecture we calculated
$\tilde{A}^{\mu\nu\alpha\beta}$ and using it, forced $g_{yy}=1+\epsilon$
component of the actual metric. $\epsilon$ being the forced strain. Using the
elastic equation (17), and working in a geometric mean-field approximation, we
solved the other terms $g_{xx}$ and $g_{xy}$ and calculate Poisson’s ratio
$\nu=-\frac{g_{xx}-1}{\epsilon}$ (see appendix B for details).

####  I.3.1 Ordered netwroks

In the ordered case we simulated a triangular lattice of varying unit cells’
shapes, and computed the angle dependence Poisson’s ratio. In this case all of
the non-affine tensors $W_{(s)\,\mu\nu}^{\alpha\beta}=0$ identically vanish,
leading to a simple calculation using eq.(10) (detailed analysis in appendix
C). In Figure 1, we see comparison of the analytical solution, and the
numerical estimation. Insets show the lattice structure.

Figure 1: Ordered networks. Simulation (yellow) vs analytical estimation
(blue) of Poisson’s ratio as a function of the angle depending on different
shape parameters $(\phi,\psi)$. Begining in the middle row, left image, and
advancing clockwise: $(\phi=1,\psi=1)$, $(2,2)$, $(3,0.9)$, $(1.5,1)$,
$(2,0.5)$, $(1,0.5)$. Insets- the spring network shape (up to rotations)

All unit cells generated from a triad of vertices
$\\{(0,0),(\phi/2,\psi\sqrt{3}/2),(0,1)\\}$. $\phi$ is called the ”shear”
factor, $\psi$ is the ”elongation” factor. The case $\phi=\psi=1$ corresponds
to an equilateral triangle with edge length of 1.

####  I.3.2 Foam-like

Following [15] we simulate a random, foam-like, network, exhibiting an auxetic
behavior at certain parameter range. To do this we produce a regular
triangular lattice with lattice constant 1. We then move each vertex in a
random direction by a fixed amount $0<\eta<0.5$. We then directly compute
$\nu(\eta)$ using the formalism described here, and, independently, using a
simulation. Both when computing, and when simulating, for each $\eta$ we
calculate several times to average the results. Our results are consistent
with the numerical simulations in [15] \- $\nu$ decreases as a function of
$\eta$ reaching $\nu=0$ at $\eta\simeq 0.46$ and reaching $\nu=-0.1$ when
$\eta\rightarrow 0.5$. In fig. 2 we show the comparison between the network
simulation (discrete triangles) and the semi-analytical computation described
in this paper.

Figure 2: Random foam-like networks. Comparison between theoretical
calculation (solid line) and numerical simulation (triangles). Solid line
represents the average of 7 different realizations of about 150 vertexes each
(shaded region is typical deviation). Numerical simulations were done over 10
realizations of 676 vertexes each, error bars mark the deviations. Note that
the last triangle is pointing down, indicating the average is beyond plot
boundaries. Inset- an example of a $\eta=0.45$ realization.

####  I.3.3 Hexagonal (honeycomb) network

We consider a honeycomb network, in which the basic hexagonal unit can vary
between regular and a re-entrant hexagon, continuously with a diameter
(distance between two opposing vertexes) of $2\ell$ $0<\ell<1$. (see inset in
figure). In such a case, Poisson’s ratio is analytically given [21]

$\displaystyle\nu(\ell)$
$\displaystyle=-\frac{\left(2-4\ell\right)\left(\ell+1/2\right)}{3+4\ell-4\ell^{2}}$
(23)

Figure 3: Hexagonal networks. Theoretical (solid line) and computational
(points) results for a honeycomb made out of uniform hexagons with diameter
$2\ell$. To avoid singular expressions, each hexagon was divided to triangles
(as indicated by dashed lines). Such that the added edges had negligible, but
finite, rigidity ($k_{dashed}=10^{-3}k_{solid}$), insets (from left t right -
re-entrant hexagon, general hexagon, regular hexagon)

In order to calculate the elastic response, we use a triangularized hexagon,
with a vertex at the center. We then set the spring constant of the radial
springs connecting the center with each corner of the hexagon to a very small
value ($1/1000$’th of the peripheral springs). This was done as the original
formulation becomes singular when the spring constant vanishes completely.
Results are described in fig. 3 showing a very high degree of agreement
between the analytical result and our formulation, despite the use of a large
difference between springs constants. This points at the strength an relevance
of our approach.

##  II Conclusions

In this work we give, to our knowledge, the first analytical derivation of the
effective, coarse grained, elastic description of a spring network. Comparison
of computational results stemming from this derivation to known/numerical
results shows a high degree of agreement. Additionally, we identified the
”non-affine” deformations, and have shown how they affect the resulting
continuum elastic model. In systems such as granular media, these quantities
play a role in local stress release by means of plastic deformations [22].

While derived over a spring-network, the results derivation shown is relevant
to many fields and systems in two ways. First the interaction between elements
is almost always approximated as that of a simple spring, especially at small
deformation. This is true for mechanical systems such as meta materials [14,
15, 16] and coarse grained mechanical models[13], chemical systems such as
self assemblies [7, 8, 9, 10] , and many biological systems as well [11, 12].
As such, the resulting theory, as is,is relevant to engineers, physicists,
chemists, and biologist, and opens the possibility of rational design of
materials.

Second, the coarse graining process described here can be generalized to
other, more complicated interaction, and is not limited to point masses
connected by linear springs. Nonlinearities can be addressed by using higher
orders terms (shape - related nonlinearities are actually addressed by the
usage of a metric description). Formulated correctly it could apply to complex
molecules, cell-cell interaction, and to polymer-networks. Activity may be
involved in it as well.

Finally, the introduction of the new, $W$, quantities invites further
investigation as to the nature of the solutions eq.(20), both analytically,
and numerically. It is known, for example, that the non-affine deformations
have a characteristic scale [22], this formulation may allow further insight
into their scale dependence. Other usages route would involve intelligent
design - relating the required mechanical behavior to the non-affine
deformations, and from that to the network structure.

##  III Acknowledgements

D.G would want to thank Amos Grossman for his help, patience and useful
discussions.

## References

  * Germain [1821] S. Germain,  _Recherches sur la théorie des surfaces élastiques_ (V. Courcier, 1821).
  * Cauchy [1826] A. L. B. Cauchy,  _Exercises de mathématiques_ , Vol. 1 (Chez De Bure Frères, 1826).
  * Martin  _et al._ [1972] P. C. Martin, O. Parodi, and P. S. Pershan, Unified hydrodynamic theory for crystals, liquid crystals, and normal fluids, Physical Review A 6, 2401 (1972).
  * Miserez  _et al._ [2022] F. Miserez, S. Ganguly, R. Haussmann, and M. Fuchs, Continuum mechanics of nonideal crystals: Microscopic approach based on projection-operator formalism, Physical Review E 106, 054125 (2022).
  * Gompper and Kroll [1996] G. Gompper and D. M. Kroll, Random surface discretizations and the renormalization of the bending rigidity, Journal de Physique I 6, 1305 (1996).
  * Kot  _et al._ [2014] M. Kot, H. Nagahashi, and P. Szymczak, Elastic moduli of simple mass spring models, The Visual Computer 31, 1339 (2014).
  * Kantor and Nelson [1987] Y. Kantor and D. R. Nelson, Phase transitions in flexible polymeric surfaces, Physical Review A 36, 4020 (1987).
  * Seung and Nelson [1988] H. S. Seung and D. R. Nelson, Defects in flexible membranes with crystalline order, Physical Review A 38, 1005 (1988).
  * Bailey and Wilson [2021] M. H. Bailey and M. Wilson, Self assembly of model polymers into biological random networks, Computational and Structural Biotechnology Journal 19, 1253 (2021).
  * Underhill and Doyle [2004] P. T. Underhill and P. S. Doyle, On the coarse-graining of polymers into bead-spring chains, Journal of Non-Newtonian Fluid Mechanics 122, 3 (2004).
  * Moshe  _et al._ [2018] M. Moshe, M. J. Bowick, and M. C. Marchetti, Geometric frustration and solid-solid transitions in model 2d tissue, Physical Review Letters 120, 268105 (2018).
  * Chen and Boyle [2014] M. Chen and F. J. Boyle, Investigation of membrane mechanics using spring networks: Application to red-blood-cell modelling, Materials Science and Engineering: C 43, 506 (2014).
  * Bolander and Saito [1998] J. Bolander and S. Saito, Fracture analyses using spring networks with random geometry, Engineering Fracture Mechanics 61, 569 (1998).
  * Rocks  _et al._ [2017] J. W. Rocks, N. Pashine, I. Bischofberger, C. P. Goodrich, A. J. Liu, and S. R. Nagel, Designing allostery-inspired response in mechanical networks, Proceedings of the National Academy of Sciences 114, 2520 (2017).
  * Liu  _et al._ [2019] J. Liu, Y. Nie, H. Tong, and N. Xu, Realizing negative poisson's ratio in spring networks with close-packed lattice geometries, Physical Review Materials 3, 055607 (2019).
  * Luo  _et al._ [2023] Y. Luo, C.-L. Ho, B. R. Helliker, and E. Katifori, Flow-network-controlled shape transformation of a thin membrane through differential fluid storage and surface expansion, Physical Review E 107, 024419 (2023).
  * Efrati  _et al._ [2009] E. Efrati, E. Sharon, and R. Kupferman, Elastic theory of unconstrained non-euclidean plates, Journal of the Mechanics and Physics of Solids 57, 762 (2009).
  * Zhang  _et al._ [2019] M. Zhang, D. Grossman, D. Danino, and E. Sharon, Shape and fluctuations of frustrated self-assembled nano ribbons, Nature Communications 10, 10.1038/s41467-019-11473-6 (2019).
  * Levin  _et al._ [2021] I. Levin, E. Siéfert, E. Sharon, and C. Maor, Hierarchy of geometrical frustration in elastic ribbons: Shape-transitions and energy scaling obtained from a general asymptotic theory, Journal of the Mechanics and Physics of Solids 156, 104579 (2021).
  * Efrati  _et al._ [2013] E. Efrati, E. Sharon, and R. Kupferman, The metric description of elasticity in residually stressed soft materials, Soft Matter 9, 8187 (2013).
  * Gibson  _et al._ [1982] L. J. Gibson, M. F. Ashby, G. S. Schajer, and C. I. Robertson, The mechanics of two-dimensional cellular materials, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 382, 25 (1982).
  * Tsamados  _et al._ [2009] M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Barrat, Local elasticity map and plasticity in a model lennard-jones glass, Physical Review E 80, 026112 (2009).

##  Appendix A Supplemental Material

##  Appendix B Methods

###  B.1 Simulation

The simulation was created for the purpose of this research. In each run we
simulated a strip with length to width ratio of 4:1. A total of about
$13\times 13\times 4=676$ vertexes, corresponding to about $1000$ edges,
depending on the exact details of each simulation.

When creating a lattice, vertexes were positioned using the base vectors -

$\displaystyle v_{1}=$ $\displaystyle\left(1,0\right)$ (24) $\displaystyle
v_{2}=$ $\displaystyle\left(\phi\frac{1}{2},\psi\frac{\sqrt{3}}{2}\right)$

where $0<\phi,\psi$ are the shear and elongation parameters, respectively, and
are used to control the shape of the triangles. $\psi=\phi=1$ corresponds to
an equilateral triangle, and any $\phi=1$ is an isosceles triangle. The strip
was created by keeping all vertexes whose coordinates satisfy $0\leq x\leq 13$
and $0\leq y\leq 52$. Different orientations lattices was created by rotating
the base vectors.

The set of vertexes was then used to create the list of edges via
triangulation, and extracting the list of neighbors. The energy of each edge
was directly calculated from the positions of its vertexes. In the simulation,
the coordinates of the top and bottom vertexes is held constant and all other
vertexes are allowed to move in order to minimize the total energy.

The two vertexes initially closest to $x=0,y=26$ and $x=13,y=26$ were
identified to measure the strain between them $\delta=\frac{\Delta
x_{final}}{\Delta x_{intial}}$, where $\Delta x_{final}$ is the final
$x$-coordinate difference between the two vertexes, and $\Delta x_{initial}$
is the initial difference. After setting the top vertexes at
$y=52(1+\epsilon)$ ($\epsilon=0.01$). And letting the system to relax
elastically, Poisson’s ratio was calculated via $\nu=-\delta/\epsilon$. And
averaged over several simulations, if required (in the more stochastic
simulations).

Simulating the foam-like structure, is stochastic in essence. We used the same
initialization process, with the following differences. After generating a
triangular lattice strip with $\phi=\psi=1$, and triangulation, we changed the
position of each vertex by an amount $0<\eta<0.5$ at a random direction, and
used the resulting distances as the reference lengths of each vertex. We then
followed the regular procedure by stretching the strip, and letting the system
relax (with the reference lengths calculated just a moment before).

###  B.2 Calculation through eq.(20)

A square patch was generated, independently of the simulation. Generation of
the network itself was done similar to the way described in the simulation.
However, once that calculated, instead of stretching the network we calculate
$A^{\mu\nu\alpha\beta}$ using the $\\{W_{(s)\,\lambda\tau}^{\mu\nu}\\}$ which
are calculated using eq.(20). Poisson’s ratio is calculate in the mean field
approximation as described in appendix B.

###  B.3 Appendix A: Multi - index notation

Here we write explicitly the expressions of the multi-index notation of the
matrices $A^{SS^{\prime}}$ and $B^{SS^{\prime}}$ that appear in the text. One
way to directly build them is two first write the matrices $W$ in terms of
their independent terms (labeled $W_{n}$) and similarly for the matrix $A$.
Then calculate the multiplication
$M^{\mu\nu\alpha\beta}=A^{\mu\nu\rho\sigma}W^{\alpha\beta}_{\rho\sigma}$. At
this point we vectorize $W$, and in using the same scheme to vectorize $M$
(which has the same number of independent components as $W$), and $\delta A$.
We now define $A^{LL^{\prime}}=\frac{\partial M_{L}}{\partial
W_{L^{\prime}}}$.

E.g. in 2D - marking

$\displaystyle A^{IJ}_{s}=$
$\displaystyle\left(\begin{array}[]{cccc}A_{s}^{1111}&A_{s}^{1112}&A_{s}^{1211}&A^{1212}\\\
A_{s}^{1121}&A_{s}^{1122}&A_{s}^{1221}&A^{1222}\\\
A_{s}^{2111}&A_{s}^{2112}&A_{s}^{2211}&A^{2212}\\\
A_{s}^{2121}&A_{s}^{2122}&A_{s}^{2221}&A^{2222}\\\
\end{array}\right)=\left(\begin{array}[]{cccc}A_{s}^{1}&A_{s}^{2}&A_{s}^{2}&A_{s}^{3}\\\
A_{s}^{2}&A_{s}^{3}&A_{s}^{3}&A_{s}^{4}\\\
A_{s}^{2}&A_{s}^{3}&A_{s}^{3}&A_{s}^{4}\\\
A_{s}^{3}&A_{s}^{4}&A_{s}^{4}&A_{s}^{5}\\\ \end{array}\right)$ (33)
$\displaystyle W^{IJ}_{s}=$
$\displaystyle\left(\begin{array}[]{cccc}W_{(s)\,11}^{11}&W_{(s)\,12}^{11}&W_{(s)\,11}^{12}&W_{(s)\,12}^{12}\\\
W_{(s)\,21}^{11}&W_{(s)\,22}^{11}&W_{(s)\,21}^{12}&W_{(s)\,22}^{12}\\\
W_{(s)\,11}^{21}&W_{(s)\,12}^{21}&W_{(s)\,11}^{22}&W_{(s)\,12}^{22}\\\
W_{(s)\,21}^{21}&W_{(s)\,22}^{21}&W_{(s)\,21}^{22}&W_{(s)\,22}^{22}\\\
\end{array}\right)=\left(\begin{array}[]{cccc}W_{s}^{1}&W_{s}^{2}&W_{s}^{4}&W_{s}^{5}\\\
W_{s}^{2}&W_{s}^{3}&W_{s}^{5}&W_{s}^{6}\\\
W_{s}^{4}&W_{s}^{5}&W_{s}^{7}&W_{s}^{8}\\\
W_{s}^{5}&W_{s}^{6}&W_{s}^{8}&W_{s}^{9}\\\ \end{array}\right)$ (42)

Where repeating markings indicate that these elements are equal. Additionally
the numbering of the $W_{s}$ matrix elements suggest a possible vectorization
scheme (which is used in this text).

It now follows by direct calculation that after vectorization

$\displaystyle A^{LL^{\prime}}_{s}=$ $\displaystyle\left(\begin{array}[]{ccc
ccc ccc}A_{s}^{1}&0&0&2A_{s}^{2}&0&0&A_{s}^{3}&0&0\\\
0&A_{s}^{1}&0&0&2A_{s}^{2}&0&0&A_{s}^{3}&0\\\
0&0&A_{s}^{1}&0&0&2A_{s}^{2}&0&0&A_{s}^{3}\\\
A_{s}^{2}&0&0&2A_{s}^{3}&0&0&A_{s}^{4}&0&0\\\
0&A_{s}^{2}&0&0&2A_{s}^{3}&0&0&A_{s}^{4}&0\\\
0&0&A_{s}^{2}&0&0&2A_{s}^{3}&0&0&A_{s}^{4}\\\
A_{s}^{3}&0&0&2A_{s}^{4}&0&0&A_{s}^{5}&0&0\\\
0&A_{s}^{3}&0&0&2A_{s}^{4}&0&0&A_{s}^{5}&0\\\
0&0&A_{s}^{3}&0&0&2A_{s}^{4}&0&0&A_{s}^{5}\\\ \end{array}\right)$ (52)

We can now finally define $A^{SS^{\prime}}$ and $B^{SS^{\prime}}$
(independently of dimension)-

$\displaystyle A^{SS^{\prime}}=$
$\displaystyle\left[diag_{n}(A^{LL^{\prime}}_{s})\right]^{S,S}=\left(\begin{array}[]{cccc}A^{LL^{\prime}}_{1}&0&0&\\\
0&A^{LL^{\prime}}_{2}&0&\cdots\\\ 0&0&A^{LL^{\prime}}_{3}&\\\
&\vdots&&\ddots\end{array}\right)$ (57) $\displaystyle B^{SS^{\prime}}=$
$\displaystyle\frac{1}{n}\left[\delta A^{LL^{\prime}}_{s}\otimes
ones_{n}\right]^{S,S}=\frac{1}{n}\left(\begin{array}[]{cccc}A^{LL^{\prime}}_{1}&A^{LL^{\prime}}_{2}&A^{LL^{\prime}}_{3}&\\\
A^{LL^{\prime}}_{1}&A^{LL^{\prime}}_{2}&A^{LL^{\prime}}_{3}&\cdots\\\
A^{LL^{\prime}}_{1}&A^{LL^{\prime}}_{2}&A^{LL^{\prime}}_{3}&\\\
&\vdots&&\ddots\end{array}\right)$ (62)

In three dimensions the same logic follows. The matrices look different, as
they are much larger.

###  B.4 Appendix B: Mean field approximation

Within a geometric mean field, $\bf{g}$ is assumed constant in space. We begin
by considering one edge lies at coordinate $y=0$, the other one at coordinate
$y=y_{m}$, and similarly there are edges at $x=0$ and $x=x_{m}$. At time zero
we deform the network in real space $\vec{R}=\left\\{X(x,y),Y(x,y)\right\\}$
so that $\left|Y(y=y_{m})-Y(y=0)\right|=L$. We take this constraint into
account by introducing an effective energy using a Lagrange multiplier
$\lambda$ which is the stress on the boundary causing the stretching of the
tissue.

$\displaystyle E_{ext}$
$\displaystyle=\int\lambda\left(\left|\vec{R}(y=y_{m})-\vec{R}(y=0)\right|-L\right)\,\text{d}x=\int\,\text{d}x\lambda\left(\left|\int\limits_{0}^{y_{m}}\partial_{y}\vec{R}\,\text{d}y\right|-L\right)\,\text{d}x$
(63)
$\displaystyle=\int\lambda\left(\sqrt{\int\int\vec{R}^{\prime}(y)\cdot\vec{R}^{\prime}(y^{\prime})\,\text{d}y^{\prime}\,\text{d}y}-L\right)\,\text{d}x.$

Within the mean field approximation,
$\vec{R}^{\prime}(y)=\partial_{y}\vec{R}(y)=const.$ So that
$\vec{R}^{\prime}(y)\cdot\vec{R}^{\prime}(y^{\prime})=g_{22}$ we can write

$\displaystyle E_{eff}$
$\displaystyle=\int\lambda\left(\sqrt{y_{m}^{2}g_{22}}-L\right)\,\text{d}x=\int\lambda\left(y_{m}\sqrt{g_{22}}-y_{m}\sqrt{G}\right)\,\text{d}x=\int{\lambda}y_{m}\left(\sqrt{g_{22}}-\sqrt{G}\right)\,\text{d}x$
(64)
$\displaystyle=\int{\lambda}\left(\sqrt{g_{22}}-\sqrt{G}\right)\,\text{d}y\,\text{d}x$

where the length scale of the stretched network is $L=\sqrt{G}x_{y}$.

Taking the variation of the total energy with respect to $g$ an $\lambda$ we
derive the following equations:

$\displaystyle\sqrt{g_{22}}-\sqrt{G}=0$ (65)
$\displaystyle\sigma^{\mu\nu}=\sigma_{el}^{\mu\nu}+\frac{\lambda}{2\sqrt{g_{22}}}\left(\begin{array}[]{cc}0&0\\\
0&1\end{array}\right)$ (68)

where $\sigma_{el}^{\mu\nu}$ is the elastic stres In principle, in order to
find the metric $g_{\mu\nu}$ minimizing the elastic energy we must find the
stress so that on the boundaries $y=0$ and $y=y_{m}$, the stress balances the
force. In the mean field approximation the conditions at the boundary impose
$\sigma^{\mu\nu}=0$. Using $G=1+\epsilon$ we get the required expression.

###  B.5 Appendix C: Ordered networks

In an ordered network, calculation is fairly simple. As
$\tilde{A}^{\mu\nu\alpha\beta}={A}^{\mu\nu\alpha\beta}={A}^{\mu\nu\alpha\beta}_{s}$
is uniform and the non-affine coefficients are 0. Where, with a little abuse
of notation

$A^{\mu\nu\alpha\beta}(x)=\langle
A^{\mu\nu\alpha\beta}\rangle=\frac{1}{N_{s}(\Omega)}\sum_{s\in\Omega(x)}A^{\mu\nu\alpha\beta}_{(s)}=\frac{1}{N_{s}(\Omega)}\sum_{s\in\Omega(x)}\sum_{e\in
s}\frac{k_{e}\Delta x_{e}^{\mu}\Delta x_{e}^{\nu}\Delta x_{e}^{\alpha}\Delta
x_{e}^{\beta}}{16\ell_{e}^{2}},$

and $\Omega(x)$ is some region around the point, x.

For an ordered systems, it is enough to consider just a small, finite sum of a
few simplex edges (in two dimension these are just triangles). We used
$k_{e}=1$, and thus
$A^{\mu\nu\alpha\beta}=\ell_{1}^{2}\left(\cos\theta_{1}\right)^{a}\left(\sin\theta_{1}\right)^{b}+\ell_{2}^{2}\left(\cos\theta_{2}\right)^{a}\left(\sin\theta_{2}\right)^{b}+\ell_{3}^{2}\left(\cos\theta_{3}\right)^{a}\left(\sin\theta_{3}\right)^{b}$,
where $a=8-\mu-\nu-\alpha-\beta$, $b=4-a$, $\ell_{1},\ell_{2},\ell_{3}$ are
the reference lengths of the edges, $\theta_{1},\theta_{2},\theta_{3}$ are the
angles of the edges with relative to the ”1” direction, and
$\mu,\nu,\alpha\beta\in(1,2)$. We thus solve for $g_{11}-1=-\nu(g_{22}-1)$.

We modeled each triangle as composed of the following vertexes (up tp a global
rotation angle $\Theta$ of the triangles relative to the stretching direction)

$\\{(0,0),(\phi/2,\psi\sqrt{3})/2,(0,1)\\}$

$\phi$ corresponds to ”shearing” of the lattice unit cell, while $\psi$
corresponds to ”elongation”. This results with the edge list

$\\{(0,1),(\phi/2,\psi\sqrt{3}/2),(1-\phi/2,-\psi\sqrt{3}/2)\\}$

and the reference lengths
$\\{1,\sqrt{\phi^{2}+3\psi^{2}}/2,\sqrt{(1-\phi)^{2}+3\psi^{2}}/2\\}$

The theoretical calculation results with a complex expression that can be
easily evaluated

$\displaystyle\nu=$
$\displaystyle-\frac{P(\phi,\psi,\Theta)}{4Q(\phi,\psi,\Theta)}$ (69)
$\displaystyle P(\phi,\psi,\Theta)=$
$\displaystyle\cos(4\theta)\left\\{27\psi^{6}-(\phi-2)^{2}\phi^{2}[\phi^{2}-2\phi+4]+9\psi^{4}(\phi^{2}-2\phi)-3\psi^{2}[(\phi^{2}-2\phi)(\phi^{2}-2\phi-20)-16]\right\\}$
$\displaystyle+\sin(4\theta)\left\\{2\sqrt{3}\psi(\phi-1)\left[9\psi^{4}+6\psi^{2}(\phi^{2}-2\phi-2)+(\phi^{2}-2\phi)(\phi^{2}-2\phi+8)\right]\right\\}-48\psi^{2}$
$\displaystyle+\left(\phi^{2}-2\phi-3\psi^{2}\right)\left[9\psi^{4}+6\psi^{2}(\phi^{2}-2\phi)+(\phi^{2}-2\phi)(\phi^{2}-2\phi+4)\right]$
$\displaystyle Q(\phi,\psi,\Theta)=$ $\displaystyle
36\psi^{4}\sin^{4}(\theta)+2\cos^{4}(\theta)\left[6\psi^{2}(\phi^{2}-2\phi)^{2}+(\phi^{2}-2\phi)^{2}(\phi^{2}-2\phi+4)+9\psi^{4}(\phi^{2}-2\phi+2)\right]$
$\displaystyle+6\psi^{2}\sin^{2}(\theta)\cos^{2}(\theta)\left(9\psi^{4}+\phi^{4}\right)+3\psi^{2}\sin^{2}(2\theta)\left[3\psi^{2}(\phi^{2}-2\phi+2)-2(\phi-2)(\phi^{2}-4\phi+2)\right]$
$\displaystyle-4\sqrt{3}\psi(\phi-1)\sin(\theta)\cos^{3}(\theta)\left(9\psi^{4}+6\psi^{2}(\phi^{2}-2\phi)+(\phi^{2}-2\phi)(\phi^{2}-2\phi+8)\right)-48\sqrt{3}\psi^{3}(\phi-1)\sin^{3}(\theta)\cos(\theta).$

For each $\phi$ $\psi$ values the results are then plotted as a function of
$\Theta$ and are compared to simulation.


#  Ab initio lattice dynamical studies of silicon clathrate frameworks and
their negative thermal expansion

Ville J. Härkönen  Antti J. Karttunen<EMAIL_ADDRESS>University of
Jyväskylä, PO Box 35, FI-40014, Finland

###### Abstract

The thermal and lattice dynamical properties of seven silicon clathrate
framework structures are investigated with ab initio density functional
methods (frameworks I, II, IV, V, VII, VIII, and H). The negative thermal
expansion (NTE) phenomenon is investigated by means of quasiharmonic
approximation and applying it to equal time displacement correlation
functions. The thermal properties of the studied clathrate frameworks,
excluding the VII framework, resemble those of the crystalline silicon diamond
structure. The clathrate framework VII was found to have anomalous NTE
temperature range up to 300 K and it is suitable for further studies of the
mechanisms of NTE. Investigation of the displacement correlation functions
revealed that in NTE, the volume derivatives of the mean square displacements
and mean square relative displacements of atoms behave similarly to the
vibrational entropy volume derivatives and consequently to the coefficients of
thermal expansion as a function of temperature. All studied clathrate
frameworks, excluding the VII framework, possess a phonon band gap or even two
in the case of the framework V.

negative thermal expansion, phonon dispersion, clathrate, ab initio lattice
dynamics, DFTP, displacement correlation function, momentum correlation
function, mean square displacement, vibrational entropy, silicon.

###### pacs:

63.20.dk, 63.70.+h, 65.40.De, 65.60.+a, 65.80.-g

##  I Introduction

In addition to the crystalline silicon diamond structure (d-Si), various other
crystal structures contain tetrahedrally bonded silicon atoms. Examples of
such structures are the experimentally known semiconducting clathrates,Kasper
_et al._ (1965) which are typically obtained as host-guest compounds such as
$\textrm{Na}_{24}\textrm{Si}_{136}$, where a porous Si clathrate II framework
(sometimes denoted as Si34 or Si136) is partially filled by Na guest
atoms.Shevelkov and Kovnir (2011) However, silicon clathrates with a
(practically) guest-free Si framework structure have been prepared, as
well.Gryko  _et al._ (2000); Ammar  _et al._ (2004) The semiconducting
clathrates have been investigated intensively due to their interesting
properties such as anomalous electronic structure (differing from that of
d-Si),Adams  _et al._ (1994); Saito and Oshiyama (1995)
superconductivityKawaji  _et al._ (1995) and relatively high thermoelectric
efficiency.Nolas  _et al._ (1998); Sootsman  _et al._ (2009); Christensen  _et
al._ (2010) In particular, the most intense research efforts on the
thermoelectric properties of clathrates have so far been directed at
germanium-based clathrates such as
$\textrm{Ba}_{8}\textrm{Ga}_{16}\textrm{Ge}_{30}$.Kuznetsov  _et al._ (2000);
Toberer  _et al._ (2008)

Tang et al. have investigated the thermal properties of the clathrate II
framework by means of ab initio calculations and experimental measurements and
have shown that the clathrate II framework exhibits negative thermal expansion
(NTE) in the temperature range of 10–140 K.Tang  _et al._ (2006) NTE is a
phenomenon in which a material contracts instead of expanding when heated to
higher temperatures. From the point of view of practical applications, NTE
materials have been used to design composite materials with zero thermal
expansion.Evans  _et al._ (1997) The NTE materials have been studied actively
after the discovery of large temperature range NTE materials such as
$\textrm{ZrW}_{2}\textrm{O}_{8}$.Evans  _et al._ (1996); Mary  _et al._
(1996), which shows NTE in a wide temperature range from nearly 0 K up to 1080
K. The underlying mechanisms causing the NTE behavior are not fully
understood. The NTE phenomenon and the different NTE mechanisms have been
reviewed for example by Barrera et al.Barrera  _et al._ (2005) d-Si also shows
NTE within temperature range from 25 K up to 120 K.Reeber and Wang (1996)
However, to our knowledge, there are no elemental semiconductors with
anomalously large NTE temperature ranges comparable to the oxides such as
$\textrm{ZrW}_{2}\textrm{O}_{8}$.

The NTE behavior of materials can be experimentally studied by means of
Extended X-Ray Absorption Fine Structure (EXAFS) spectroscopy, which enables
the measurement of local dynamics of the crystal lattice.Sayers  _et al._
(1971); Lee and Pendry (1975); Hayes  _et al._ (1976); Tröger  _et al._ (1994)
Quantities which can be obtained from EXAFS spectra are, for example, mean
square displacements (MSD) and mean square relative displacements (MSRD).
Consequently, EXAFS has also been used to study the NTE phenomenon.Sanson  _et
al._ (2006); el All  _et al._ (2012) However, computational studies of the NTE
phenomenon from the point of view of MSRD are not common. To our knowledge,
there are no previous studies where the NTE behavior of materials is
investigated by means of full ab initio calculation of MSRDs.

In this paper we investigate the thermal properties and negative thermal
expansion of seven different silicon clathrate frameworks by means of ab
initio lattice dynamics. The thermal properties of the materials are
determined within quasiharmonic approximation (QHA).Huang and Born (1954) The
results are compared with available experimental data and previous
computational predictions for d-Si to establish the accuracy of the
computational methods. We use equal time displacement correlation functions to
investigate the mechanisms of the thermal expansion and to elucidate the
negative thermal expansion behavior of the silicon clathrate frameworks.

##  II Theory, computational methods and studied structures

###  II.1 Lattice dynamics and thermal expansion

To obtain the various thermal properties, the phonon eigenvalues and phonon
eigenvectors have been extracted with ab initio methods, by applying density
functional perturbation theory (DFPT).Baroni  _et al._ (2001) The DFTP
formalism enables the calculation of dynamical matrix (introduced by the
relations below), after which the phonon eigenvalues and eigenvectors can be
solved from an eigenvalue equation written within the harmonic approximation
asMaradudin  _et al._ (1971)

$\omega^{2}_{j}\left(\mathbf{q}\right)e_{\alpha}\left(\kappa|\mathbf{q}j\right)=\sum_{\kappa^{\prime},\beta}D_{\alpha\beta}\left(\kappa\kappa^{\prime}|\mathbf{q}\right)e_{\beta}\left(\kappa^{\prime}|\mathbf{q}j\right),$
(1)

where $\alpha$ and $\beta$ are cartesian indices, $\kappa$ and
$\kappa^{\prime}$ are atom indices within the primitive unit cell,
$\omega_{j}\left(\mathbf{q}\right)$ is phonon eigenvalue (referred as phonon
frequency from now on) for phonon branch $j$ and wave vector $\mathbf{q}$,
$e_{\alpha}\left(\kappa|\mathbf{q}j\right)$ are the components of the phonon
eigenvector and $D_{\alpha\beta}\left(\kappa\kappa^{\prime}|\mathbf{q}\right)$
are the components of the dynamical matrix. The dynamical matrix is the
Fourier transform of the real space force constant matrix
$\Phi_{\alpha\beta}\left(l\kappa;l^{,}\kappa^{,}\right)$

$\displaystyle
D_{\alpha\beta}\left(\kappa\kappa^{\prime}|\mathbf{q}\right)=\left(M_{\kappa}M_{\kappa^{\prime}}\right)^{-1/2}\sum_{l^{\prime}}\Phi_{\alpha\beta}\left(l\kappa;l^{\prime}\kappa^{\prime}\right)$
$\displaystyle\times{}e^{-i\mathbf{q}\cdot\left[\mathbf{x}\left(l\right)-\mathbf{x}\left(l^{\prime}\right)\right]},$
(2)

where $M_{\kappa}$ and $M_{\kappa^{\prime}}$ are the atomic masses of atoms
$\kappa$ and $\kappa^{\prime}$ respectively, and $\mathbf{x}\left(l\right)$ is
the position vector of the $l_{\textrm{th}}$ unit cell. The components of the
eigenvector $\mathbf{e}\left(\kappa|\mathbf{q}j\right)$ are chosen to satisfy
the orthonormality and closure conditionsHuang and Born (1954)

$\sum_{\kappa,\alpha}e_{\alpha}\left(\kappa|\mathbf{q}j^{\prime}\right)e^{*}_{\alpha}\left(\kappa|\mathbf{q}j\right)=\delta_{jj^{\prime}},$
(3)
$\sum_{j}e_{\alpha}\left(\kappa|\mathbf{q}j\right)e^{*}_{\beta}\left(\kappa^{\prime}|\mathbf{q}j\right)=\delta_{\alpha\beta}\delta_{\kappa\kappa^{\prime}},$
(4)

where $\delta_{\alpha\beta}$ is the Kronecker delta. The volumetric
coefficient of thermal expansion (CTE) $\alpha_{V}$ is defined as

$\alpha_{V}=\frac{1}{V_{0}}\frac{\partial{V}}{\partial{T}},$ (5)

where $V$ is the volume, $V_{0}$ is the equilibrium volume and $T$ the
temperature. The linear CTE is a measure of the change in particular dimension
of the crystal with respect to temperature. In the cubic crystal system,
linear CTE $\alpha_{L}$ in the direction of the principal axis can be written
as

$\alpha_{L}=\frac{1}{L_{0}}\frac{\partial{L}}{\partial{T}}=\frac{1}{3V_{0}}\frac{\partial{V}}{\partial{T}}=\frac{1}{3}\alpha_{V}.$
(6)

In non-cubic crystals, the anisotropy of the crystal structure must be taken
in account. That is, the different components of the thermal expansion along
different principal axes must be calculated separately to obtain the linear
CTEs. In this work only volumetric CTEs have been calculated for non-cubic
structures. Using QHA, the expression for the volumetric CTE can be written in
terms of isothermal bulk modulus $B$, phonon mode Grüneisen parameters
$\gamma_{j}\left(\mathbf{q}\right)$ and phonon mode heat capacities at
constant volume $c_{v,j}\left(\mathbf{q}\right)$Ashcroft and Mermin (1976)

$\alpha_{V}=\frac{1}{B}\sum_{\mathbf{q},j}c_{v,j}\left(\mathbf{q}\right)\gamma_{j}\left(\mathbf{q}\right),$
(7)

The heat capacity can be expressed in terms of Bose-Einstein distribution
function $\overline{n}_{\mathbf{q},j}$ and phonon frequencies
$\omega_{j}\left(\mathbf{q}\right)$

$c_{v,j}\left(\mathbf{q}\right)=\frac{\omega_{j}\left(\mathbf{q}\right)}{V_{0}}\frac{\partial}{\partial{T}}\overline{n}_{\mathbf{q},j},$
(8)

and the Grüneisen parameters can be expressed as

$\gamma_{j}\left(\mathbf{q}\right)=-\frac{V_{0}}{\omega_{0,j}\left(\mathbf{q}\right)}\frac{\partial{\omega_{j}}\left(\mathbf{q}\right)}{\partial{V}},$
(9)

where $\omega_{0,j}\left(\mathbf{q}\right)$ is the phonon frequency for the
phonon mode $j$ and wave vector $\mathbf{q}$ at the equilibrium volume
$V_{0}$. The Grüneisen parameters can also be expressed asTaylor (1998)

$\displaystyle\gamma_{j}\left(\mathbf{q}\right)=-\frac{1}{2}\frac{V_{0}}{\omega^{2}_{0,j}\left(\mathbf{q}\right)}\sum_{\alpha,\beta}\sum_{\kappa,\kappa^{\prime}}e^{*}_{\alpha}\left(\kappa|\mathbf{q}j\right)$
$\displaystyle\times{}\frac{\partial
D_{\alpha\beta}\left(\kappa\kappa^{\prime}|\mathbf{q}\right)}{\partial{V}}e_{\beta}\left(\kappa^{\prime}|\mathbf{q}j\right).$
(10)

Eq. 7 shows that when NTE occurs, the Grüneisen parameters must have negative
values since the bulk modulus and the heat capacities have positive values.
Furthermore, from Eq. 9 it can be seen that to have negative Grüneisen
parameter values, the phonon frequencies must decrease with decreasing volume.
Analogously, Eq. 10 implies that in NTE the components of the dynamical matrix
have smaller values with decreasing volume. Finally, an alternative way to
express $\alpha_{V}$ is to make use of thermodynamical relations and to write
its expression in terms of vibrational entropyBarrera  _et al._ (2005)

$\alpha_{V}=\frac{1}{B}\sum_{\mathbf{q},j}\frac{\partial{S_{j}\left(\mathbf{q}\right)}}{\partial{V}},$
(11)

From Eq. 11 it can be seen that the vibrational entropy $S$ must increase with
decreasing volume for NTE to occur.

###  II.2 Correlation functions

Mean square displacements (MSD) are a special case of the displacement time
correlation functions derived by for example Maradudin.Maradudin  _et al._
(1971) Equal time displacement correlation function can be written as

$\displaystyle\left\langle
u_{\alpha}\left(l\kappa\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle=\frac{\hbar}{2N\left(M_{\kappa}M_{\kappa^{\prime}}\right)^{1/2}}$
$\displaystyle\times{}\sum_{\mathbf{q},j}\frac{e_{\alpha}\left(\kappa|\mathbf{q}j\right)e^{*}_{\beta}\left(\kappa^{\prime}|\mathbf{q}j\right)}{\omega_{j}\left(\mathbf{q}\right)}$
$\displaystyle\times{}e^{i\mathbf{q}\cdot\left[\mathbf{x}\left(l\right)-\mathbf{x}\left(l^{\prime}\right)\right]}\coth\left[\frac{\hbar\omega_{j}\left(\mathbf{q}\right)}{2k_{B}T}\right],$
(12)

where $N$ is the number of unit cells (and $\mathbf{q}$-points) and $k_{B}$ is
the Boltzmann constant. In the special case of $l=l^{\prime}$ and
$\kappa=\kappa^{\prime}$, Eq. 12 is called the equal time displacement
autocorrelation function, vibrational amplitude, or MSD. From now on, all
discussed displacement correlation functions are equal time displacement
correlation functions and corresponding autocorrelation functions, if not
otherwise mentioned. To analyze the NTE phenomenon, it is useful to know how
the atoms move relative to each other. Mean square relative displacement
(MSRD) in the interatomic direction between atoms $\left(l\kappa\right)$ and
$\left(l^{\prime}\kappa^{\prime}\right)$, also known as the parallel MSRD, is
defined asBeni and Platzman (1976)

$\left\langle
u^{2}_{\parallel}\right\rangle=\left\langle\left[\left(\mathbf{u}\left(l\kappa\right)-\mathbf{u}\left(l^{\prime}\kappa^{\prime}\right)\right)\cdot\mathbf{\hat{r}}\right]^{2}\right\rangle,$
(13)

where $\mathbf{\hat{r}}$ is unit vector in the direction of vector between
atoms $\left(l\kappa\right)$ and $\left(l^{\prime}\kappa^{\prime}\right)$.
Only the difference of $l$ and $l^{\prime}$ matters, not their absolute
values, thus $l=0$ is chosen. Writing out the component form of Eq. 13, the
parallel MSRD is

$\displaystyle\left\langle
u^{2}_{\parallel}\right\rangle=\sum_{\alpha,\beta}r_{\alpha}r_{\beta}\left.[\left\langle
u_{\alpha}\left(0\kappa\right)u_{\beta}\left(0\kappa\right)\right\rangle\right.$
$\displaystyle\left.+\left\langle
u_{\alpha}\left(l^{\prime}\kappa^{\prime}\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle-2\left\langle
u_{\alpha}\left(0\kappa\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle\right.].$
(14)

Perpendicular MSRD is defined asFornasini (2001)

$\left\langle u^{2}_{\perp}\right\rangle=\left\langle
u^{2}_{R}\right\rangle-\left\langle u^{2}_{\parallel}\right\rangle,$ (15)

where

$\left\langle
u^{2}_{R}\right\rangle=\left\langle\left\|\mathbf{u}\left(0\kappa\right)-\mathbf{u}\left(l^{\prime}\kappa^{\prime}\right)\right\|^{2}\right\rangle,$
(16)

which is the total MSRD of atoms $\left(l\kappa\right)$ and
$\left(l^{\prime}\kappa^{\prime}\right)$. Writing Eq. 16 in the component form
and substituting to Eq. 15, perpendicular MSRD can be written as

$\displaystyle\left\langle
u^{2}_{\perp}\right\rangle=\sum_{\alpha,\beta}\left(\delta_{\alpha\beta}-r_{\alpha}r_{\beta}\right)\left.[\left\langle
u_{\alpha}\left(0\kappa\right)u_{\beta}\left(0\kappa\right)\right\rangle\right.$
$\displaystyle\left.+\left\langle
u_{\alpha}\left(l^{\prime}\kappa^{\prime}\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle-2\left\langle
u_{\alpha}\left(0\kappa\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle\right.].$
(17)

Similar calculations involving the displacement correlation functions have
been previously carried out by Nielsen and Weber, who used empirical
interatomic potentials to investigate the correlation functions for
d-Si.Nielsen and Weber (1980); Weber (1987) Finally, the total MSD of atom
$\left(l\kappa\right)$ is the sum of the autocorrelation functions along the
different cartesian components

$\left\langle
u^{2}\right\rangle=\left\langle\left\|\mathbf{u}\left(l\kappa\right)\right\|^{2}\right\rangle=\sum_{\alpha}\left\langle
u_{\alpha}\left(l\kappa\right)u_{\alpha}\left(l\kappa\right)\right\rangle.$
(18)

In a cubic crystal, the total MSD (within numerical accuracy) is due to
symmetry reasons equal to $3\left\langle
u_{\alpha}\left(l\kappa\right)u_{\alpha}\left(l\kappa\right)\right\rangle$,
for any $\alpha=1,2,3$.

Because the preceding displacement correlation functions are valid only within
harmonic approximation, they cannot be used to describe the thermal expansion
phenomenon as such. If QHA is assumed to be valid, the relation (Eq. 9) for
the volume dependence of the phonon frequencies can be used. Inspection of Eq.
12 shows that the two terms including the hyperbolic cotangent and the inverse
of the frequency are always positive. As in Sec. II.1, from Eq. 7 and Eq. 9 it
can be seen that in NTE, the change in the frequencies must be in the same
direction as the change in the volume. In general, the change of Eq. 12 with
respect to volume also depends on the behaviour of the phonon eigenvector and
its components $e_{\alpha}\left(\kappa|\mathbf{q}j\right)$ as a function of
volume. We show first with heuristic deduction that the volume dependence of
Eq. 18 in the case of d-Si is only due to the volume dependence of the
frequencies. Eq. 18 reads

$\left\langle
u^{2}\right\rangle=C\sum_{\alpha}\sum_{\mathbf{q},j}\left|e_{\alpha}\left(\kappa|\mathbf{q}j\right)\right|^{2}\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right],$
(19)

where

$C=\frac{\hbar}{2NM_{\kappa}},$ (20)

and

$\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right]=\omega_{j}\left(\mathbf{q}\right)^{-1}\coth\left[\frac{\hbar\omega_{j}\left(\mathbf{q}\right)}{2k_{B}T}\right].$
(21)

The frequencies are not dependent on the cartesian index $\alpha$, and
$\xi\left(\omega_{j}\left(\mathbf{q}\right)\right)$ can be moved out of the
sum over $\alpha$

$\left\langle
u^{2}\right\rangle=C\sum_{\mathbf{q},j}\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right]\sum_{\alpha}\left|e_{\alpha}\left(\kappa|\mathbf{q}j\right)\right|^{2}.$
(22)

From Eq. 3 it follows that the sum over $\alpha$ in Eq. 22 with any chosen $j$
and $\mathbf{q}$ must have value $1/2$, in the case of d-Si with two
equivalent atoms in the primitive unit cell. Eq. 18 for d-Si is then

$\left\langle
u^{2}\right\rangle_{\textit{d}-Si}=\frac{C}{2}\sum_{\mathbf{q},j}\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right].$
(23)

Preceding is also true for structures where the lattice constants have been
displaced in a way that preserves the original symmetry and it is also
confirmed for d-Si by means of a direct calculation. This implies that in NTE,
the MSD of a particular atom $\left(l\kappa\right)$ and a particular mode $j$
increases with decreasing volume. The reason for this behaviour becomes
clearer when the hyperbolic cotangent is written in terms of the Bose-Einstein
distribution function $\overline{n}_{\mathbf{q}j}$. Using the definition of
the hyperbolic cotangent, the last term in Eq. 23 is

$\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right]=\omega_{j}\left(\mathbf{q}\right)^{-1}\left(2\overline{n}_{\mathbf{q},j}+1\right),$
(24)

In NTE, Eq. 11 and Eq. 23 actually have their values because of the same
underlying reason. In the following, the similarity between these quantities,
MSD for d-Si and $S$, in NTE is considered when the temperature increases. If
NTE occurs, the phonon frequency of the state $\left(\mathbf{q}j\right)$
contributing to NTE must have lower values when the volume decreases. The same
effect occurs, when the considered phonon mode corresponds to a positive CTE
and the volume increases. For the phonon modes corresponding to NTE, the
preceding implies that there are more lower frequency states for phonons to
occupy ($\overline{n}_{\mathbf{q}j}$ have higher values for lower
$\omega_{j}\left(\mathbf{q}\right)$ in a particular temperature $T$).
Therefore, the entropy increases and according to Eq. 11, NTE occurs when the
opposite effect of the phonon modes corresponding to positive CTE is weaker.
Similar conclusion applies to Eq. 23 in NTE. If the state
$\left(\mathbf{q}j\right)$ has lower values of frequency when the volume
decreases, there are more low frequency states that the phonons can occupy and
$\xi\left[\omega_{j}\left(\mathbf{q}\right),T\right]$ has larger values,
resulting in larger MSD values for the phonon modes contributing to NTE. The
opposite is true for the phonon modes contributing to positive CTE. For both
negative and positive CTE, the total MSD increases with increasing
temperature. As a summary, in NTE, the MSD and entropy of the phonon modes
corresponding to positive CTE increase less than the MSD and entropy of the
phonon modes corresponding to NTE.

To study the effect of thermal expansion phenomenon to the MSD as a function
of temperature, we define the following dimensionless quantity

$\left\langle
u^{2}\right\rangle_{r}\equiv\frac{1}{N_{a}}\left(\frac{V_{0}}{\left\langle
u^{2}\right\rangle_{V_{0}}}\frac{\partial\left\langle
u^{2}\right\rangle}{\partial V}-\left\langle
u^{2}\right\rangle_{r,T_{0}}\right),$ (25)

which includes the volume derivatives of the total MSD, and where

$\left\langle u^{2}\right\rangle_{r,T_{0}}=\frac{V_{0}}{\left\langle
u^{2}\right\rangle_{V_{0},T_{0}}}\frac{\partial\left\langle
u^{2}\right\rangle_{T_{0}}}{\partial V}.$ (26)

In Eq. 25, $N_{a}$ is the number of atoms within the primitive unit cell,
$V_{0}$ is the volume of the equilibrium structure, $\left\langle
u^{2}\right\rangle_{V_{0}}$ is the total MSD for the equilibrium structure,
and $T_{0}$ indicates the value of the function at $T=0$K. Eq. 25 measures the
product of the total MSD volume derivative and the ratio of the equilibrium
volume and the equilibrium values of the total MSD. The values of Eq. 25 are
therefore negative for modes for which the values of the frequency decrease as
the volume decreases, at a particular temperature $T$. This is based on the
discussions above and the fact that the MSD (Eq. 18) have only positive values
for all values of $T$. The values of the first term within the brackets in Eq.
25 differ from zero when $T\rightarrow 0$ K. For clarity, a normalization
factor is defined (the second term within the brackets in Eq. 25), which has a
constant value for any chosen state $\left(\mathbf{q}j\right)$. The
normalization factor ensures that Eq. 25 is equal to zero when $T\rightarrow
0$ K. The preceding normalization works properly since the first term within
the brackets in Eq. 25 increases (positive CTE) or decreases (negative CTE)
with all relevant $T>0$. In addition, a second normalizing factor $N_{a}$ is
introduced to facilitate comparisons between structures with different-sized
primitive unit cells. $N_{a}$ is related to the number of phonon modes $N_{j}$
by the relation $N_{j}=3N_{a}$. The same scheme is implemented for the
parallel and perpendicular MSRD to study changes in the relative displacements
with respect to volume and their role in the NTE phenomenon. They are defined
as

$\left\langle
u^{2}_{\parallel}\right\rangle_{r}\equiv\frac{1}{N_{a}}\left(\frac{V_{0}}{\left\langle
u^{2}_{\parallel}\right\rangle_{V_{0}}}\frac{\partial\left\langle
u^{2}_{\parallel}\right\rangle}{\partial V}-\left\langle
u^{2}_{\parallel}\right\rangle_{r,T_{0}}\right),$ (27)

and

$\left\langle
u^{2}_{\perp}\right\rangle_{r}\equiv\frac{1}{N_{a}}\left(\frac{V_{0}}{\left\langle
u^{2}_{\perp}\right\rangle_{V_{0}}}\frac{\partial\left\langle
u^{2}_{\perp}\right\rangle}{\partial V}-\left\langle
u^{2}_{\perp}\right\rangle_{r,T_{0}}\right).$ (28)

The preceding quantities in Eqs. 25 and 27–28 are named as relative MSD and
MSRD volume derivatives, and these quantities are denoted as MSD-VD and MSRD-
VD, respectively.

Properties analogous to the MSD-VD and MSRD-VD introduced above can be derived
also for equal time momentum correlation functions, which can be written
within harmonic approximation asMaradudin  _et al._ (1971)

$\displaystyle\left\langle
p_{\alpha}\left(l\kappa\right)p_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle=\frac{\hbar}{2N}\left(M_{\kappa}M_{\kappa^{\prime}}\right)^{1/2}$
$\displaystyle\times{}\sum_{\mathbf{q},j}e_{\alpha}\left(\kappa|\mathbf{q}j\right)e^{*}_{\beta}\left(\kappa^{\prime}|\mathbf{q}j\right)\omega_{j}\left(\mathbf{q}\right)$
$\displaystyle\times{}e^{i\mathbf{q}\cdot\left[\mathbf{x}\left(l\right)-\mathbf{x}\left(l^{\prime}\right)\right]}\coth\left[\frac{\hbar\omega_{j}\left(\mathbf{q}\right)}{2k_{B}T}\right]=$
$\displaystyle
M_{\kappa}M_{\kappa^{\prime}}\omega^{2}_{j}\left(\mathbf{q}\right)\left\langle
u_{\alpha}\left(l\kappa\right)u_{\beta}\left(l^{\prime}\kappa^{\prime}\right)\right\rangle,$
(29)

where $p_{\alpha}\left(l\kappa\right)$ is the momentum for the atom
$\left(l\kappa\right)$ in the direction $\alpha$. The difference is that the
change of the mean square momentum (MSM) with respect to the volume, derived
from Eq. 29, is opposite in direction to the MSD-VD and MSRD-VD. This means
that if the volume decreases and the temperature increases for example by 1 K
(NTE occurs), the phonon modes contributing to NTE (positive CTE) resist
(assist) the increase of the MSM, but since the changes in the frequencies are
small (
$\left|\Delta\omega_{j}\left(\mathbf{q}\right)\right|/\omega_{j}\left(\mathbf{q}\right)<<1$),
MSM increases in both cases.

Finally, few remarks about the method of calculating CTE from Eq. 7 or 11 and
the related properties calculated from Eq. 25 are necessary. In general, the
Grüneisen parameters are not independent of the temperature, which also means
that $\omega_{j}\left[\mathbf{q},V\left(T\right),T\right]$. When the
temperature independence is assumed it follows from Eq. 7 and 11 that the
phonon modes corresponding to positive CTE (NTE) have positive (negative)
contribution to CTE at all finite $T>0$. The same features result for the MSD-
VD and MSRD-VD and the related quantities derived for momentum.

###  II.3 Studied structures and computational details

The structural characteristics of the studied seven silicon clathrate
frameworks are illustrated in Fig. 1 and described in Table 1.

Figure 1: Schematic figures of the silicon clathrate structures included in this work. The vertices of the polyhedral cages present silicon atoms. Crystallographic unit cell edges are drawn in black. For a more comprehensive description of the structural characteristics, see.Karttunen  _et al._ (2010) (Color online) Table 1: Structural data and computational details for the studied structures. Structure | Space Group | Atoms/cell111The number of atoms in the primitive cell. | Elect. (k1,k2,k3)222The mesh used for the electronic $\mathbf{k}$-sampling. | Phon. (q1,q2,q3)333The meshes used for phonon calculations, PDOS calculations, and correlation function calculations, respectively. The meshes for the PDOS and correlation function calculations were Fourier interpolated from the meshes used in the phonon calculations. | PDOS (q1,q2,q3)c | Corr. (q1,q2,q3)c
---|---|---|---|---|---|---
d-Si | $Fd\bar{3}m\left(227\right)$ | 2 | 16,16,16 | 8,8,8 | 80,80,80 | 20,20,20
I | $Pm\bar{3}n\left(223\right)$ | 46 | 4,4,4 | 4,4,4 | 30,30,30 | -
II | $Fd\bar{3}m\left(227\right)$ | 34 | 4,4,4 | 4,4,4 | 30,30,30 | 8,8,8
IV | $P6/mmm\left(191\right)$ | 40 | 4,4,4 | 3,3,3 | 30,30,30 | -
V | $P6_{3}/mmc\left(194\right)$ | 68 | 4,4,3 | 4,4,3 | 30,30,30 | -
VII | $Im\bar{3}m\left(229\right)$ | 6 | 8,8,8 | 8,8,8 | 60,60,60 | -
VIII | $I\bar{4}3m\left(217\right)$ | 23 | 6,6,6 | 4,4,4 | 40,40,40 | -
H | $P6/mmm\left(191\right)$ | 34 | 4,4,4 | 3,3,3 | 30,30,30 | -

The semiconducting clathrate structures are classified according to the
polyhedral cages they are composed of. The structural properties of the
clathrates have been reviewed in e.g. Rogl (2006); Kovnir and Shevelkov
(2004); Shevelkov and Kovnir (2011)

The ab initio density functional calculations to optimize the structures and
to calculate the phonon dispersion relations were carried out with the Quantum
Espresso program package (QE, version 5.0.3).Giannozzi  _et al._ (2009) The Si
atoms were described by Martins-Troullier-type norm-conserving
pseudopotentials Troullier and Martins (1991); Dal Corso (2012) and the Local
Density Approximation (LDA) was used as the exchange-correlation energy
functional.Perdew and Zunger (1981) In all calculations, the following kinetic
energy cutoffs have been used: 40 Ry for wavefunctions and 160 Ry for charge
densities and potentials. The applied $\mathbf{k}$\- and $\mathbf{q}$-sampling
for each studied structure is listed in Table 1. The $\mathbf{q}$-meshes for
the phonon density of states (PDOS) and correlation function calculations were
Fourier interpolated from the mesh used in the corresponding phonon
calculation (QE module matdyn.x).Baroni  _et al._ (2001) We carried out
convergence tests for both SCF and phonon frequency calculations with
different $\mathbf{k}$-meshes. The lattice constant for the optimized
structure was practically identical for (8,8,8) and (16,16,16)
$\mathbf{k}$-meshes and the final energy was about 5*10-3 eV lower with the
(16,16,16) mesh. The maximum differences between the phonon frequencies
obtained with the (8,8,8) and (16,16,16) $\mathbf{k}$-meshes are about 1.5%.
Since a reasonable convergence was found for the (8,8,8) $\mathbf{k}$-mesh in
the case of d-Si, the $\mathbf{k}$-meshes listed in in Table 1 for the
clathrate frameworks were chosen as a compromise between accuracy and
computational cost. For d-Si, the use of the denser (16,16,16)
$\mathbf{k}$-mesh was still computationally feasible. We have also confirmed
the convergence of the calculated phonon frequencies with respect to the
applied $\mathbf{q}$-meshes. Both the lattice constants and the atomic
positions of the studied structures were fully optimized (applying the space
group symmetries listed in Table 1). The Grüneisen parameters were calculated
by displacing the structures from the equilibrium lattice constant by
$\pm$0.5%, optimizing the atomic positions, and calculating the numerical
derivatives using central differences. Bulk moduli were calculated by
displacing the structures from the equilibrium lattice constant up to $\pm$2%
with a step size of $\pm$0.5% and calculating the bulk modulus using the
Murnaghan’s equation of state (QE module ev.x.). Phonon density of states
values were calculated using the tetrahedron method (QE module
matdyn.x).Blöchl  _et al._ (1994) In all cases where a summation over
$\mathbf{q}$ is involved, the appropriate weights to normalize the sums have
been used. The algorithms for calculating the thermodynamic quantities and the
correlation functions were implemented as Matlab routines.

Calculation of the correlation functions involves some numerical issues, which
need to be discussed. To our knowledge, there is no way to label the phonon
modes uniquely in a point of degeneracy when the diagonalization of the
dynamical matrix in Eq. 1 is done numerically. This is true also for the QE
program package and probably results in minor numerical inaccuracies in the
calculation of the correlation functions. Some phonon labeling could be
carried out by continuity arguments for the phonon eigenvectors in particular
direction when using relatively dense $\mathbf{q}$-meshes, but in the absence
of a rigorous general approach, the labeling of the phonon modes is left to
the default algorithm in QE. The numerical accuracy is also sensitive to the
applied Fast Fourier Transform (FFT) mesh. In this work, the default settings
of QE were used to choose the FFT grid.

##  III Results and discussion

###  III.1 Optimized structures and phonon dispersion relations

Table 2 lists the lattice constants, the relative energies $\Delta E$ with
respect to d-Si, and the bulk moduli for the optimized silicon clathrate
frameworks.

Table 2: Lattice constants $a,c$, relative energies ($\Delta E$), and bulk moduli ($B$) for the optimized structures. Structure |  $a$(Å) |  $c$(Å) |  $\Delta E$ (eV/atom) |  $B$ (GPa)
---|---|---|---|---
d-Si | 5.38 | $-$ | 0.00 | 96.5
I | 10.08 | $-$ | 0.09 | 84.1
II | 14.52 | $-$ | 0.08 | 82.4
IV | 10.07 | 10.23 | 0.11 | 78.5
V | 10.26 | 16.80 | 0.08 | 82.1
VII | 6.63 | $-$ | 0.27 | 72.6
VIII | 9.95 | $-$ | 0.10 | 86.6
H | 10.36 | 8.43 | 0.11 | 79.6

As could be expected, the structure II, which has been synthetized as a nearly
guest-free silicon framework, shows the lowest relative energy among the
studied clathrate frameworks. It is followed closely by structure V, which is
a hypothetical hexagonal modification of the structure II.Karttunen  _et al._
(2010) The predicted relative energies are in good agreement with previous
results obtained with the PBE0 hybrid density functional method and localized
Gaussian type basis sets (GTO).Karttunen  _et al._ (2010) The structure VII
shows the largest difference in $\Delta E$ in comparison to the PBE0 results
(0.03 eV/atom), while for the other structures the difference is only 0.01
eV/atom. The structure VII is a rather strained one with a high relative
energy and it is not experimentally as relevant as the other studied
structures. The predicted bulk moduli are systematically somewhat smaller in
comparison to previous results obtained with the PBE0 functional,Karttunen
_et al._ (2011) but the result for d-Si is similar to previous LDA-
predictions.Nielsen and Martin (1985) The previous LDA predictions for the
clathrate frameworks I and II (87 and 87.5 GPa, respectively) are also
comparable to the present results.San-Miguel  _et al._ (2002) LDA slightly
underestimates the bulk modulus ($\approx$1.3% in the case of d-Si) in
comparison to experiment, the experimental value for d-Si being 97.8 GPa.Hall
(1967) The difference between calculated and experimental lattice constant
value is approximapdftely 1%.Basile  _et al._ (1994) The experimental value of
the bulk modulus for the silicon clathrate framework II has been found to be
$90\pm 5$ GPa.San-Miguel  _et al._ (1999)

Figures 2 and 3 illustrate the phonon dispersion relations in the studied
structures and the corresponding phonon density of states (PDOS).

Figure 2: Phonon dispersion relations along high symmetry paths in the first
Brillouin zone for a$)$ d-Si (experimental values indicated by diamondsNilsson
and Nelin (1972)) and b$)$ for structure VII. Figure 3: Phonon dispersion
relations along high symmetry paths in the first Brillouin zone for structures
I, II, IV, V, VIII, and H.

Fig. 2 also shows the experimental phonon dispersion data for d-Si.Nilsson and
Nelin (1972) The LDA results for d-Si are in relatively good agreement with
the experiment, the differences being the largest for large wave vectors (the
max. difference is approximately 7%). The estimated experimental error is
$\pm$1 cm-1.Nilsson and Nelin (1972) The LDA phonon dispersion data for the
clathrate framework II are also in agreement with previous computational
results.Tang  _et al._ (2006) The PDOS of the silicon clathrate frameworks I
and II have also been determined experimentally using inelastic neutron
scattering for structures containing some potassium and sodium impurities
(e.g. $\textrm{Na}:\textrm{Si}$ ratio 1:34). Despite the impurities, the
experimental results show similar main features when compared to the results
shown in Fig. 3, such as the high PDOS values at the frequencies of
approximately 175 cm-1.Mélinon  _et al._ (1999) In the structure VII (Fig.
2b), a significant decrease of the lowest optical modes towards the acoustic
modes can be observed, in particular for small wave vectors. This indicates
strong anomalies in the interatomic forces. The other clathrate frameworks
show very similar phonon dispersion relations compared to each other. The
structures other than d-Si and structure VII have a phonon band gap within the
frequency range of 400–475 cm-1, the structure V even showing two phonon band
gaps at frequencies of about 440 and 390 cm-1 (Fig. 3).

###  III.2 Thermal properties from the quasiharmonic approximation

The calculated heat capacities at constant volume are shown in Fig. 4.

Figure 4: Heat capacities at constant volume for studied structures. Only the
structure VII is slightly different from the other structures (indicated with
a dash-dotted line).

The results are very similar for all structures, the strained structure VII
showing the largest difference to the other structures. The results for the
structure II are in agreement with previous experimental and computational
results.Biswas  _et al._ (2008) The reason for the deviation observed for the
heat capacity of structure VII can be seen from Fig. 5, where Eq. 8 has been
plotted.

Figure 5: Heat capacity at constant volume as a function of temperature and
frequency in QHA. Every separate line represents the value of Eq. 8 for a
particular phonon frequency as a function of the temperature.

At low $T$, only the low frequency modes affect the heat capacity
significantly. When $T$ increases, the contribution of the low frequency modes
does not increase as rapidly as the contribution of the higher frequency
modes. Since the structure VII has more low frequency modes and less higher
frequency modes in comparison to the other structures, it has a slightly
higher heat capacity at low temperature and lower heat capacity at higher
temperatures. This also affects the CTE of the structure VII. The very similar
heat capacities suggest that in the studied structures, and within QHA, the
Grüneisen parameters and the differences in the values of the bulk moduli have
the greatest effect on the magnitude of CTE.

Fig. 6 shows the Grüneisen parameters for the studied structures as a function
of the phonon frequency.

Figure 6: Grüneisen parameters for the studied structures, each dot represents
one mode at a particular $\mathbf{q}$-point.

For every structure, approximately 75000 frequency values have been used to
generate the Grüneisen parameter data. The Grüneisen parameters predicted for
d-Si and the structure II are in agreement with previous computational
results.Tang  _et al._ (2006) The structure VII has distinctly different
Grüneisen parameters from the other structures, showing negative values that
are approximately 10 times larger than for the other structures. An
interesting feature is that the Grüneisen parameters only have negative values
when the phonon frequency is under a certain threshold that is 220 cm-1 for
d-Si and VII and 190 cm-1 for the other structures. However, in every
structure, the acoustic modes with frequencies significantly below 100 cm-1 do
have positive Grüneisen parameter values at some $\mathbf{q}$. For the
longitudal phonon mode of d-Si, the Grüneisen parameter values are positive
for all wave vectors. On the other hand, the two transverse acoustic phonon
modes have positive values of Grüneisen parameter only at few wave vectors.
The Grüneisen parameters of the clathrate frameworks show features similar to
amorphous silicon,Fabian and Allen (1997) where the Grüneisen parameter values
are spread out at low frequency values. This is in contrast to d-Si, where the
Grüneisen parameters show a rather symmetric shape as a function of frequency
in comparison to the studied clathrate frameworks and amorphous silicon.

The results in Figures 4 and 6 suggest that only the structure VII should show
significantly different CTE in comparison to the other studied structures.
This is confirmed by the linear and volumetric CTE values shown in Fig. 7 for
all studied structures.

Figure 7: Top: Linear CTE for the studied cubic structures. The experimental
values for d-Si are indicated with circles.Reeber and Wang (1996) Bottom:
Volumetric CTE for the studied hexagonal structures, d-Si is included as a
reference.

For d-Si, the experimental values for the linear CTE are also shown.Reeber and
Wang (1996) The CTE values obtained from Eq. 7 and Eq. 11 are identical. In
comparison to experiment, LDA overestimates the magnitude of the NTE for d-Si
and underestimates the CTE at positive values. In all studied structures,
excluding the structure VII, the CTE values are very similar. The structure
VII with distinctly different Grüneisen parameters shows NTE behavior up to
approximately 300 K and the NTE is nearly five times larger than for d-Si. As
discussed above, the structure VII has a very high relative energy and it is
not experimentally as relevant as the other structures. However, it is still
an interesting test case for detailed studies on the NTE phenomenon due to the
small size of its primitive unit cell and its anomalous lattice dynamical
properties.

Fig. 8 illustrates the change in entropy $\Delta S$ with respect to the change
in the volume $\Delta V$ for each phonon mode as a function of temperature for
d-Si, and the clathrate frameworks II and VIII.

Figure 8: $\Delta S/\Delta V$ for each phonon mode in the d-Si, II, and VIII
structures as a function of the temperature. The zero level is at the point
where the lines representing the phonon modes coincide at $T=0$ K. The
acoustic modes are indicated with dash-dotted lines.

In d-Si, the longitudinal acoustic mode has a positive $\Delta S/\Delta V$ at
all temperatures, having a positive contribution to CTE (Eq. 11). The
transverse acoustic modes have a negative contribution to CTE, giving rise to
the experimentally confirmed NTE behavior of d-Si. In the structures II and
VIII, according to QHA, all acoustic modes have a negative contribution to CTE
at all considered temperatures. The reason for this was discussed at the end
of Sec. II.2. An interesting feature of the structures II and VIII is that
several optical modes (all optical modes in the structure II) have larger
contributions to NTE than the acoustic modes when $T>$100 K. Similar crossing
occurs between optical phonon modes having a positive contribution to CTE, as
shown in Fig. 8. The difference is that the crossing occurs at $T\approx
200$K. For the studied structures, the maximum value of $\Delta S/\Delta V$,
when the temperature changes by $\pm 1$K, is six orders of magnitude smaller
than the absolute value of $S$ at the same temperature.

###  III.3 Thermal expansion and displacement correlation functions

Fig. 9 shows the displacement autocorrelation functions (Eq. 12) for all atoms
within the primitive unit cells of d-Si and clathrate II.

Figure 9: Displacement autocorrelation functions for d-Si and the clathrate
framework II (Eq. 12). Each line represents the displacement autocorrelation
function for one atom in the primitive unit cell.

The results for d-Si are in good agreement with previous results of Nielsen
and Weber obtained using empirical potentials,Nielsen and Weber (1980) and
with the results of Strauch et al. obtained using ab initio methods.Strauch
_et al._ (1996) The differences between d-Si and the clathrate framework II
are small. In the case of the structure II, the values are separated in three
different groups corresponding to the three different Wyckoff positions in the
crystal structure. The primitive unit cell of the structure II is shown in
Fig. 10 to facilitate the discussion on the results related to the
displacement correlation functions.

Figure 10: The primitive unit cell of the clathrate framework II shown in a)
[100] b) [111] and c) general direction. Unit cell edges are drawn in black.

The highest MSD values correspond to atoms 1 and 2 in Fig. 10 (Wyckoff
position 8a), the intermediate values to atoms 3-10 (32e), and the smallest
values to atoms 11-34 (96g). The MSD values calculated for the symmetry
equivalent atoms in the structure II show a minor deviation from each other.
The highest deviations correspond to atoms 11-34, when the max. difference
within this group is approximately 1.8%, which occurs at $T=500$ K. One
possible reason for this is the issue of phonon mode labeling for degenerate
modes during the diagonalization of the dynamical matrix, as discussed in the
end of Sec. II.3. The total MSD (Eq. 18) calculated for the structure VII is
approximately 25% higher than the total MSDs calculated for d-Si and the
structure II.

The parallel and perpendicular displacement correlation functions between the
nearest neighbours in d-Si ($l\kappa=l1;l\kappa^{\prime}=l2$) are presented in
Fig. 11.

Figure 11: Parallel and perpendicular correlation functions for d-Si within
same unit cell, i.e. for the nearest neigbours. Figure 12: Parallel and
perpendicular MSRD within the same unit cell for the structure II. a)
Normalized parallel MSRD. b) Parallel MSRD and c) perpendicular MSRD.
Different lines represent different atoms $\kappa$, the second atom
$\kappa^{\prime}=2$ in all cases (see Fig. 10).

The perpendicular MSRDs have larger values than the parallel MSRDs. Similar
behavior has been confirmed with experimental EXAFS studies and computational
results for d-Ge,Dalba  _et al._ (1999),Nielsen and Weber (1980) and also for
other materials such as CdTe with zincblende structure using EXAFS.el All  _et
al._ (2012) In cubic crystal structures, the perpendicular MSRD values must be
larger than the parallel MSRD values due to symmetry reasons. If a three-
dimensional Cartesian coordinate system is considered to have one of its
orthonormal basis vectors along the vector between the nearest neighbour atoms
in d-Si, the difference in the magnitudes of perpendicular and parallel MSRD
becomes more evident, as the perpendicular vibrations can be represented as
linear combinations of the two other basis vectors. Because the MSDs of the
atoms have the same values for each cartesian component, it is expected that
the relative displacements are distributed among the perpendicular and
parallel MSRD as shown in Fig. 11.

Fig. 12 shows the normalized parallel and perpendicular MSRD values for the
clathrate framework II. All results presented in Fig. 12 are for
$\kappa^{\prime}=2$ and $\kappa$ runs over all atoms within the same unit
cell. The normalized parallel correlation functions have again similar values
as in the previous works for d-Si obtained using empirical potentials.Nielsen
and Weber (1980) Both the parallel and perpendicular MSRD have similar
characteristics, that is, the nearest atoms have the smallest MSRD values and
vice versa. The atoms $\kappa=7-10$ have the smallest values and the atoms
$\kappa=1,3$ the largest values in all cases in Fig. 12. The atoms
$\kappa=1,3$ are on the opposite side of the cavity with respect to atom
$\kappa^{\prime}=2$ (see Fig. 10), which may partly explain the higher values
of the MSRD in comparison to the other atoms. The perpendicular MSRD values do
not form as distinct groups for the three symmetry equivalent atoms as in the
case of the parallel MSRD.

Fig. 13 shows the MSD-VD and MSRD-VD values for d-Si. The values in Fig. 13
for the different modes and for the total volume derivatives have similar
characteristics as CTE calculated from Eq. 11.

Figure 13: MSD-VD- and MSRD-VD values for d-Si calculated from Eq. 25, Eq. 27,
and Eq. 28. The MSRD-VD values are for the nearest neighbours. The acoustic
modes are indicated with a dash-dotted line.

The reason for this similarity was discussed at the end of Sec. II.2. The MSD-
VD values for the phonon modes corresponding to NTE are negative within the
whole temperature range similar to the entropies shown in Fig. 8. Like in the
case of the entropy, the MSD-VD values are negative for two acoustic phonon
modes, while one acoustic mode and the optical modes have a positive
contribution to CTE. The similarity of the MSRD-VD and CTE results could be
due to the fact that in d-Si only the relative motion of the two atoms within
the primitive unit cell atoms determines the CTE. According to our results,
the parallel and perpendicular MSRD-VD have nearly identical values.

The MSD-VD values for the clathrate framework II shown in Fig. 14 correspond
relatively well to the values for d-Si.

Figure 14: MSD-VD for three different atoms in the structure II calculated
from Eq. 25. Acoustic modes are indicated with a dash-dotted line.

Only three atoms representing the different symmetry independent positions in
the crystal structure are shown. Comparison of Fig. 8 and Fig. 14 implies that
the MSD-VD and $\Delta S/\Delta V$ behave similarly as a function of
temperature also for the structure II. The MSD-VD results also agree with the
results obtained for CTE in the sense that the structure II has slightly lower
MSD-VD values in both positive CTE and NTE temperature ranges when comparing
it to d-Si. The proper treatment of MSRD-VD values in larger structures
requires symmetry considerations.

##  IV Conclusions

We have investigated the thermal properties of silicon clathrate frameworks by
combining ab initio DFT lattice dynamics with QHA. The computational results
for d-Si were in relatively good agreement with the experimental values for
the phonon dispersion relations, heat capacities, and CTE. Studied clathrate
frameworks have similar phonon spectra, only the structure VII differing from
the others significantly. All studied clathrate frameworks were found to
possess a phonon band gap and similar threshold frequencies confining the
phonons that can have negative Grüneisen parameter values. A similar threshold
frequency was also found for d-Si. Furthermore, the clathrate framework V has
two phonon band gaps. The NTE temperature range of the studied clathrate
frameworks is slightly wider than for d-Si. The MSD and MSRD values calculated
for d-Si were in agreement with the previous computational results. The
results obtained from the displacement correlation functions for the clathrate
framework II were similar to d-Si. Equal time displacement correlation
functions were used to study the NTE phenomenon. The MSD-VD and MSRD-VD
quantities defined here behave almost identically to the vibrational entropy
volume derivatives in the studied structures. The parallel and perpendicular
MSRD-VD results for d-Si indicate that both quantities behave in similar
fashion as a function of temperature. Further application of the volume
derivatives of the displacement correlation functions may produce useful
information on the mechanisms of the NTE phenomenon.

###### Acknowledgements.

We thank Prof. Andrea Dal Corso (SISSA, Trieste) for providing modified QE
routines to enable efficient distributed calculation of the phonon dispersion
relations. We also thank Prof. Robert van Leeuwen and Dr. Gerrit Groenhof
(University of Jyväskylä) for useful discussions on various aspects of the
present work. The computing resources for this work were provided by the
Finnish Grid Infrastructure (FGI). We gratefully acknowledge funding from the
Foundation for Research of Natural Resources in Finland (grant 17591/13) and
the Academy of Finland (grant 138560/2010).

## References

  * Kasper  _et al._ (1965) J. S. Kasper, P. Hagenmuller, M. Pouchard, and C. Cros, Science 150, 1713 (1965).
  * Shevelkov and Kovnir (2011) A. Shevelkov and K. Kovnir,  _Zintl Clathrates_ (Springer Berlin Heidelberg, 2011) pp. 97–142.
  * Gryko  _et al._ (2000) J. Gryko, P. F. McMillan, R. F. Marzke, G. K. Ramachandran, D. Patton, S. K. Deb, and O. F. Sankey, Phys. Rev. B 62, R7707 (2000).
  * Ammar  _et al._ (2004) A. Ammar, C. Cros, M. Pouchard, N. Jaussaud, J.-M. Bassat, G. Villeneuve, M. Duttine, M. Ménétrier, and E. Reny, Solid State Sci. 6, 393 (2004).
  * Adams  _et al._ (1994) G. B. Adams, M. O’Keeffe, A. A. Demkov, O. F. Sankey, and Y.-M. Huang, Phys. Rev. B 49, 8048 (1994).
  * Saito and Oshiyama (1995) S. Saito and A. Oshiyama, Phys. Rev. B 51, 2628 (1995).
  * Kawaji  _et al._ (1995) H. Kawaji, H.-o. Horie, S. Yamanaka, and M. Ishikawa, Phys. Rev. Lett. 74, 1427 (1995).
  * Nolas  _et al._ (1998) G. Nolas, J. Cohn, G. Slack, and S. Schujman, Appl. Phys. Lett. 73, 178 (1998).
  * Sootsman  _et al._ (2009) J. Sootsman, D. Chung, and M. Kanatzidis, Angew. Chem. Int. Ed. 48, 8616 (2009).
  * Christensen  _et al._ (2010) M. Christensen, S. Johnsen, and B. B. Iversen, Dalton Trans. 39, 978 (2010).
  * Kuznetsov  _et al._ (2000) V. L. Kuznetsov, L. A. Kuznetsova, A. E. Kaliazin, and D. M. Rowe, J. Appl. Phys. 87, 7871 (2000).
  * Toberer  _et al._ (2008) E. S. Toberer, M. Christensen, B. B. Iversen, and G. J. Snyder, Phys. Rev. B 77, 075203 (2008).
  * Tang  _et al._ (2006) X. Tang, J. Dong, P. Hutchins, O. Shebanova, J. Gryko, P. Barnes, J. K. Cockcroft, M. Vickers, and P. F. McMillan, Phys. Rev. B 74, 014109 (2006).
  * Evans  _et al._ (1997) J. Evans, T. Mary, and A. Sleight, Physica B 241-–243, 311 (1997).
  * Evans  _et al._ (1996) J. S. O. Evans, T. A. Mary, T. Vogt, M. A. Subramanian, and A. W. Sleight, Chem. Mater. 8, 2809 (1996).
  * Mary  _et al._ (1996) T. A. Mary, J. S. O. Evans, T. Vogt, and A. W. Sleight, Science 272, 90 (1996).
  * Barrera  _et al._ (2005) G. D. Barrera, J. A. O. Bruno, T. H. K. Barron, and N. L. Allan, J. Phys. 17, R217 (2005).
  * Reeber and Wang (1996) R. R. Reeber and K. Wang, Mater. Chem. Phys. 46, 259 (1996).
  * Sayers  _et al._ (1971) D. E. Sayers, E. A. Stern, and F. W. Lytle, Phys. Rev. Lett. 27, 1204 (1971).
  * Lee and Pendry (1975) P. A. Lee and J. B. Pendry, Phys. Rev. B 11, 2795 (1975).
  * Hayes  _et al._ (1976) T. M. Hayes, P. N. Sen, and S. H. Hunter, J. Phys. C 9, 4357 (1976).
  * Tröger  _et al._ (1994) L. Tröger, T. Yokoyama, D. Arvanitis, T. Lederer, M. Tischer, and K. Baberschke, Phys. Rev. B 49, 888 (1994).
  * Sanson  _et al._ (2006) A. Sanson, F. Rocca, G. Dalba, P. Fornasini, R. Grisenti, M. Dapiaggi, and G. Artioli, Phys. Rev. B 73, 214305 (2006).
  * el All  _et al._ (2012) N. A. el All, G. Dalba, D. Diop, P. Fornasini, R. Grisenti, O. Mathon, F. Rocca, B. T. Sendja, and M. Vaccari, J. Phys. 24, 115403 (2012).
  * Huang and Born (1954) K. Huang and M. Born,  _Dynamical Theory of Crystal Lattices_ (Clarendon Press Oxford, 1954) p. 298.
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Maradudin  _et al._ (1971) A. Maradudin, E. W. Montroll, G. H. Weiss, and I. P. Ipatova,  _Theory of The Lattice Dynamics in The Harmonic Approximation_ , Vol. Supplement 3 (Academic Press, 1971) p. 62.
  * Ashcroft and Mermin (1976) N. Ashcroft and D. Mermin,  _Solid State Physics_ (Harcourt, 1976).
  * Taylor (1998) R. E. Taylor,  _Thermal Expansion of Solids_ , Vol. 4 (ASM international, 1998).
  * Beni and Platzman (1976) G. Beni and P. M. Platzman, Phys. Rev. B 14, 1514 (1976).
  * Fornasini (2001) P. Fornasini, J. Phys. 13, 7859 (2001).
  * Nielsen and Weber (1980) O. H. Nielsen and W. Weber, J. Phys. C 13, 2449 (1980).
  * Weber (1987) W. Weber, J. Phys. F: Metal Phys. 17, 27 (1987).
  * Karttunen  _et al._ (2010) A. J. Karttunen, T. F. Fässler, M. Linnolahti, and T. A. Pakkanen, Inorg. Chem. 50, 1733 (2010).
  * Rogl (2006) P. Rogl,  _Thermoelectrics Handbook: Macro to Nano_ (CRC Press, 2006) pp. 32–1.
  * Kovnir and Shevelkov (2004) K. A. Kovnir and A. V. Shevelkov, Russ. Chem. Rev. 73, 923 (2004).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. 21, 395502 (2009).
  * Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  * Dal Corso (2012) A. Dal Corso, http://qe-forge.org/gf/project/pslibrary/ (2012).
  * Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Blöchl  _et al._ (1994) P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 16223 (1994).
  * Karttunen  _et al._ (2011) A. J. Karttunen, V. J. Härkönen, M. Linnolahti, and T. A. Pakkanen, J. Phys. Chem. C 115, 19925 (2011).
  * Nielsen and Martin (1985) O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3792 (1985).
  * San-Miguel  _et al._ (2002) A. San-Miguel, P. Mélinon, D. Connétable, X. Blase, F. Tournus, E. Reny, S. Yamanaka, and J. P. Itié, Phys. Rev. B 65, 054109 (2002).
  * Hall (1967) J. J. Hall, Phys. Rev. 161, 756 (1967).
  * Basile  _et al._ (1994) G. Basile, A. Bergamin, G. Cavagnero, G. Mana, E. Vittone, and G. Zosi, Phys. Rev. Lett. 72, 3133 (1994).
  * San-Miguel  _et al._ (1999) A. San-Miguel, P. Kéghélian, X. Blase, P. Mélinon, A. Perez, J. P. Itié, A. Polian, E. Reny, C. Cros, and M. Pouchard, Phys. Rev. Lett. 83, 5290 (1999).
  * Nilsson and Nelin (1972) G. Nilsson and G. Nelin, Phys. Rev. B 6, 3777 (1972).
  * Mélinon  _et al._ (1999) P. Mélinon, P. Kéghélian, A. Perez, B. Champagnon, Y. Guyot, L. Saviot, E. Reny, C. Cros, M. Pouchard, and A. J. Dianoux, Phys. Rev. B 59, 10099 (1999).
  * Biswas  _et al._ (2008) K. Biswas, C. W. Myles, M. Sanati, and G. S. Nolas, J. Appl. Phys. 104, 033535 (2008).
  * Fabian and Allen (1997) J. Fabian and P. B. Allen, Phys. Rev. Lett. 79, 1885 (1997).
  * Strauch  _et al._ (1996) D. Strauch, P. Pavone, N. Nerb, K. Karch, W. Windl, G. Dalba, and P. Fornasini, Physica B 219–220, 436 (1996).
  * Dalba  _et al._ (1999) G. Dalba, P. Fornasini, R. Grisenti, and J. Purans, Phys. Rev. Lett. 82, 4240 (1999).