A phenomenological expression to describe the temperature dependence of pressure-induced softening in negative thermal expansion materials

Perhaps the simplest model to have NTE is the one shown in figure 2, where each rigid-square unit rotates against each other on heating while the unit itself is totally rigid for thermal excitation, resulting in the reduction of distance between the near square centers hence the thermal contraction of the whole system.

Zoom In Zoom Out Reset image size

The energy of such a system in an independent-object model where coupled terms are ignored (as suggested by [1, 4]) can be written as,

$$\begin{equation}\label {eq:energy} E = E_{{pL}} + E_{\mathrm {rot}} + E_{\mathrm {st}} \end{equation} \tag{ 2 }$$

with

$$\begin{equation}\label {eq:energypl} E_{{pL}} = - 2p\left ({1 - \cos \theta } \right ). \end{equation} \tag{ 3 }$$

Because of the cell reduction by a factor of (1 −cosθ) due to rotation of the square with a tilting angle of θ, this term lowers the energy for compression or rotations. E_rot describes the energy cost of the rotation with even powers of θ and K₁, K₂ > 0,

$$\begin{equation}\label {eq:energyrot} E_{\mathrm {rot}} = K_1 \theta ^2 + \frac {{K_2 }}{{12}}\theta ^4 +\cdots . \end{equation} \tag{ 4 }$$

The fourth-order terms of θ in E_pL and E_rot are related to the anharmonicity of the thermally excited rotation modes. E_st is the strain energy inside the squares under pressure p, where the bonds from the centre to the apex of the square can for example correspond to Zr/W–O bonds in ZrW₂O₈ or Zn–C/N bonds in Zn(CN)₂.

With a partition function calculated using the energy in equation (2), the expectation value of θ² at temperature τ = k_BT is,

$$\begin{eqnarray} \langle {\theta ^2 } \rangle _{\tau } &=& \frac {{\int _0^{{\pi /2 } } {\theta ^2 \exp \left ({ - E/\tau } \right ) \mathrm {d}\theta }}}{{\int _0^{{\pi /2 } } {\exp \left ({ - E/\tau } \right )\mathrm {d}\theta } }} \nonumber \\ &=& \frac {{\int _0^{{\pi /2} } {\theta ^2 \exp \left \{{ - \left [{\left ({K_1 - p} \right ) \left ({\theta ^2 /2} \right ) + \left ({K_2 + p} \right ) \left ({\theta ^4 /12} \right )} \right ]/\tau } \right \}\mathrm {d}\theta } }}{{\int _0^{{\pi /2 }} {\exp \left \{{ - \left [{\left ({K_1 - p} \right )\left ({\theta ^2 /2} \right ) + \left ({K_2 + p} \right ) \left ({\theta ^4 /12} \right )} \right ]/\tau } \right \}\mathrm {d}\theta } }},\nonumber \\ &&\label {eq:expecttheta} \end{eqnarray} \tag{ 5 }$$

where the strain energy E_st is independent from θ. cosθ is expanded as cosθ = 1 − θ²/2! + θ⁴/4! . Using a mean-field theory [5], one can approximate the fourth-order term of θ as

$$\begin{equation}\label {eq:anharmtheta} \theta ^4 = 3\langle {\theta ^2} \rangle \theta ^2 = 12C\tau \theta ^2, \end{equation} \tag{ 6 }$$

where the factor of 3 comes from the numerical test that 〈θ⁴〉 = 3〈θ²〉², and we put 〈θ²〉≈ 4Cτ in the classical approximation with positive coefficient C (the factor of 4 is for derivation convenience). By inserting this in equation (5) and putting

$$\begin{eqnarray}&&\label {eq:coefficient} \begin{array}{c} \left ({K_1 - p} \right ) = \chi _1 \\ [4pt] C\left ({K_2 + p} \right ) = \chi _2, \end{array} \end{eqnarray} \tag{ 7 }$$

one can rewrite the equation as

$$\begin{eqnarray}\label {eq:expecttheta2} \langle {\theta ^2 } \rangle _{\tau } &=& \frac {{\int _0^{{\pi /2 }} {\theta ^2 \exp \left \{{ - \left ({\chi _1 + \chi _2 \tau } \right ) \theta ^2 /\tau } \right \}\mathrm {d}\theta } }}{{\int _0^{{\pi /2 }} {\exp \left \{{ - \left ({\chi _1 + \chi _2 \tau } \right )\theta ^2 /\tau } \right \} \mathrm {d}\theta } }} \nonumber \\ &=& \frac {{\left [{\left ({\sqrt \pi /4} \right )c^{ - 3/2} } \right ]}} {{\left [{\left ({\sqrt \pi /2} \right )c^{ - 1/2} } \right ]}} \nonumber \\ &=& \frac {1}{2}\frac {\tau }{{\chi _1 + \chi _2 \tau }}, \end{eqnarray} \tag{ 8 }$$

where c = χ₁ + χ₂τ /τ. Note that the positive coefficient C (see equation (6)) and hence χ₂ defined by equation (7) provide a measure for the magnitude of the mode anharmonicity. We put the equilibrium cell length as unit 1 for convenience so that the pressure p and the force constants K₁ as well as K₂ (in equation (7)) all have energy dimensions in the model. The coefficient C has dimension of 1/energy (in equation (6)). χ₁ has energy dimension and χ₂ is dimensionless.

In such a model, for a small pressure p, the reduction of cell length (ΔL: length, area and volume reductions in one, two and three dimensions, respectively) at temperature τ is proportional to the expectation value 〈θ²〉_τ plus the compressed strain e ≈−κp < 0 of the bonds inside the square,

$$\begin{equation}\label {eq:cellchange} \Delta L \approx - \langle {\theta ^2 } \rangle _{\tau } + e \approx - \frac {1}{2}\frac {\tau }{{\chi _1 + \chi _2 \tau }} - \kappa p \end{equation} \tag{ 9 }$$

with

$$\begin{equation}\label {eq:kappa} \kappa \approx a - bp \quad\qquad {(a > 0,b > 0)} \end{equation} \tag{ 10 }$$

as the compressibility of the bond that decreases on compression. a and b (in equation (10)) have dimensions of 1/energy and 1/energy², respectively.

With the model being established, we can now discuss the temperature dependence of the bulk modulus (B₀) and its pressure derivative ($B^{\prime }_0$). The bulk modulus can be calculated from equation (9) as,

$$\begin{eqnarray}\label {eq:bulkmodulus} B &=& - L \frac {{\partial p}}{{\partial L}} \nonumber \\ &=& \frac {{\left ({\chi _1 + \chi _2 \tau } \right )^2 }} {{\tau - C\tau ^2 + \kappa \left ({\chi _1 + \chi _2 \tau } \right )^2 }}, \end{eqnarray} \tag{ 11 }$$

where we have used ∂χ₁/∂p = −1 and ∂χ₂/∂p = C according to equation (7).

For zero anharmonicity, when χ₂ = C = 0, equation (11) becomes

$$\begin{equation}\label {eq:bulkmharm} B = \frac {{\chi _1 ^2 }}{{\tau + \kappa \chi _1^2 }} \end{equation} \tag{ 12 }$$

with an asymptotic form of

$$\begin{equation}\label {eq:bulkmharmasymp} B\left ({\tau ^{\prime } \to 0} \right ) = \frac {1}{\kappa } \left ({1 - \frac {{\tau _{\mathrm {Z}} + \tau ^{\prime }}}{{\kappa \chi _1 ^2 }}} \right ), \end{equation} \tag{ 13 }$$

where we have put τ = τ_Z + τ^' with τ_Z a small constant representing the contribution from the quantum zero-point energy. Thus, B₀ (defined as B at p = 0) with χ₁ = K₁ and κ = a starts as the modulus of the bond (1/κ) at zero temperature and decreases as 1/τ on heating. The schematic plot of B₀(T) is shown as the dash curve in figure 3(a).

Zoom In Zoom Out Reset image size

The first derivative of B with respect to pressure can be calculated from equation (12) as

$$\begin{equation}\label {eq:Bdharm} B^{\prime } = \frac {{\partial B}}{{\partial p}} = \frac {{b\chi _1 ^4 - 2\chi _1 \tau }}{{\left ({\tau + \kappa \chi _1 ^2 } \right )^2 }}. \end{equation} \tag{ 14 }$$

At zero temperature, $B^{\prime }_0$ (defined as $B^{\prime }_0$ at p = 0) with χ₁ = K₁ and κ = a is

$$\begin{equation}\label {eq:Bdzerot} B^{\prime }_0 \left ({\tau ^{\prime } \to 0} \right ) \approx \frac {b}{{a^2 }} \end{equation} \tag{ 15 }$$

which shows that the value of $B^{\prime }_0$ at zero temperature is positive and equal to the pressure derivative of the bond modulus (see equation (10)). However, $B^{\prime }_0$ will become negative on heating as shown in equation (14). It also has an asymptotic form of

$$\begin{equation}\label {eq:Bdharmasymp} B^{\prime }_0 \left ({\tau \to \infty } \right ) \approx - \frac {{2K_1 }}{\tau }, \end{equation} \tag{ 16 }$$

meaning that $B^{\prime }_0$ tends to zero with τ→∞. The schematic plot of $B^{\prime }_0(T)$ is shown as the dash curve in figure 3(b).

In the case of small anharmonicity, such that C and hence χ₂ are both small quantities with negligible high-order terms of O(C²) and $O(\chi _2^2)$, the bulk modulus from equation (11) becomes

$$\begin{eqnarray}B &=& \frac {{\chi _1 ^2 + 2\chi _1 \chi _2 \tau }}{{\tau + \kappa \chi _1 ^2 + 2\chi _1 \chi _2 \tau - C\tau ^2 }} \nonumber \\ &\approx & \frac {{\chi _1 ^2 }}{{\tau + \kappa \chi _1 ^2 }}\left ({1 - \frac {{2\kappa \chi _1 \chi _2 \tau - C\tau ^2 }}{{\tau + \kappa \chi _1 ^2 }}} \right ) + \frac {{2\chi _1 \chi _2 \tau }}{{\tau + \kappa \chi _1 ^2 }}.\label {eq:Banharm} \end{eqnarray} \tag{ 17 }$$

Note that B₀ ≈ 1/κ at zero temperature. The value of B₀ will increase linearly with τ→∞. These are shown by the solid curve in figure 3(a).

From equation (17), one can further obtain,

$$\begin{eqnarray} &&B^{\prime } = \frac {{\partial B}}{{\partial p}} = \frac {{b\chi _1 ^4 - 2\chi _1 \tau }} {{\left ({\tau + \kappa \chi _1 ^2 } \right )^2 }} \left ({1 - \frac {{2\kappa \chi _1 \chi _2 \tau - C\tau ^2 }} {{\tau + \kappa \chi _1 ^2 }}} \right ) \nonumber \\ &&-\, \frac {{\chi _1 ^2 }}{{\tau + \kappa \chi _1 ^2 }} \left [{ \frac {{\left ({2\kappa \chi _1 C\tau - 2b\chi _1 \chi _2 \tau -2\kappa \chi _2\tau } \right )\left ({\tau + \kappa \chi _1 ^2 } \right )}}{{\left ({\tau + \kappa \chi _1 ^2 } \right )^2 }}} \right ] \nonumber \\ &&-\, \frac {{\chi _1 ^2 }}{{\tau + \kappa \chi _1 ^2 }} \left [{ \frac {{\left ({b\chi _1 ^2 + 2\kappa \chi _1 } \right )\left ({2\kappa \chi _1 \chi _2 \tau - C\tau ^2 } \right )}} {{\left ({\tau + \kappa \chi _1 ^2 } \right )^2 }}} \right ] \nonumber \\ &&+\, \frac {{\left ({ - 2 \chi _2 \tau + 2\chi _1 C\tau } \right )\left ({\tau + \kappa \chi _1 ^2 } \right ) + 2\chi _1 \chi _2 \tau \left ({b\chi _1 ^2 + 2\kappa \chi _1 } \right )}}{{\left ({\tau + \kappa \chi _1 ^2 } \right )^2 }}.\label {eq:Bdanharm} \end{eqnarray} \tag{ 18 }$$

At zero temperature, $B^{\prime }_0$ is the same as in equation (15). However, with the anharmonic terms, the value of $B^{\prime }_0$ will tend to zero more rapidly with temperature than the one without anharmonicity. For the case of large mode anharmonicity, the value may even become positive again. These are shown by the solid and the dash-dot curves in figure 3(b) compared to the dash curve without anharmonicity.

To sum up, the value of $B^{\prime }_0$ at zero temperature is positive and equal to the pressure derivative of the modulus of the bond, b/a² (see equation (10)). With elevated temperature, $B^{\prime }_0$ becomes negative and tends to zero with τ→∞ if the mode anharmonicity is not large. Thus, in both cases with and without anharmonicity, a phenomenological expression as illustrated in figure 1 is a good yet simple description for the temperature dependence of $B^{\prime }_0$. At zero temperature this form gives a positive value $\beta$, and at high temperature this expression shows an asymptotic behaviour similar to −exp(−T) that tends to zero.

Besides its simplicity, the advantage of the phenomenological expression of equation (1) is that its three variables can be assigned with physical meanings. $\beta$ is the value of $B^{\prime }_0$ at zero temperature. It is the pressure derivative of the bond modulus. T₀ is the temperature at which the value of $B^{\prime }_0$ drops to zero, corresponding to τ₀ that makes equation (18) equal to zero,

$$\begin{equation}\label {eq:T0} T_0 = \tau _0 \approx \frac {{bK_1 ^3 }}{2} + C\frac {(bK_1^3)^2}{2}. \end{equation} \tag{ 19 }$$

Note that with zero anharmonicity, i.e. C = K₂ = 0, T₀ is equal to the solution of equation (14) as expected. In the model, since the bond modulus is 1/κ = 1/(a − bp), 1/a is its value at zero pressure and equal to the force constant k of the bond potential. Thus, one may rewrite equation (19) as

$$\begin{equation}\label {eq:T02} T_0 \approx \left ({\frac {\beta }{2}}\right ) \left ({\frac {{K_1 }}{k}} \right )^2 K_1 + C \left ({\frac {\beta ^2}{2}} \right ) \left ({\frac {{K_1 }}{k}} \right )^4 K_1 ^2, \end{equation} \tag{ 20 }$$

where the ratio between the force constants of the rotational term of the rigid unit modes (RUMs) and the bond term, K₁/k, gives a measure of the relative stiffness of the two. The force constant K₁ can correspond to a mean frequency of the RUMs that have significant negative mode Grüneisen parameters leading to NTE of the system. Therefore, the value of T₀ is akin to the flexibility of the structure and the characteristic frequency of the RUMs; more structure flexibility and smaller RUM frequency would result in a smaller value of T₀. It is also clear that T₀ will increase due to the anharmonic term measured by C.

T₁ corresponds to the temperature at which $B^{\prime }_0$ reaches its minimum. From equation (18) one can obtain

$$\begin{eqnarray}T_1 &\approx & aK_1 ^2 + bK_1 ^3 - C \left (aK_1^2 \right )^2 \nonumber \\ &=& \left (\frac {K_1}{k}\right )K_1 + \beta \left ({\frac {K_1}{k}}\right )^2 K_1 - C \left ({\frac {K_1}{k}}\right )^2 K_1^2\label {eq:T1} \end{eqnarray} \tag{ 21 }$$

which shows that T₁ is also determined by the flexibility of the structure and the characteristic frequency of the RUMs. By comparison with equation (20), we can see that the value of T₁ should be larger than two times T₀, given that the anharmonic coefficient C is a small quantity.

Note that, other than the anharmonic terms (involving C) in T₀ and T₁, the dominant terms in these two quantities are related to the value of $\beta$ which is material-dependent. This makes the value of T₀ or T₁ not such a good indicator for comparing anharmonicity between different materials. To cancel the dependence of the dominant term on $\beta$, the difference between T₁ and two times T₀ is calculated,

$$\begin{equation}\label {eq:T1T0} T_1 - 2T_0 = \left (\frac {K_1}{k}\right )K_1 - C \left (\frac {K_1}{k}\right )^2 K_1^2 - C\beta ^2 \left (\frac {K_1}{k}\right )^4K_1^2. \end{equation} \tag{ 22 }$$

This value is positive and shows how fast $B^{\prime }_0$ can reach its minimum on heating and provides a better measurement of the anharmonicity of the system, i.e. the larger the anharmonicity C, the smaller the value of T₁ − 2T₀, as shown in figure 3.

As a demonstration, we now apply the phenomenological expression of equation (1) to some silicious zeolite structures that show isotropic NTE [4], including FAU, KFI, LTA and TSC [6]. We have carried out molecular dynamics (MD) simulations for the structures in [4] using DL_POLY [7], with O–O, Si–O and Si–Si interactions described by Coulomb interactions and Buckingham potentials ϕ(r) of the form

$$\begin{equation} \phi _{ij}(r_{ij}) = A_{ij} \exp (-r_{ij}/\rho _{ij}) - C_{ij}r^{-6}_{ij}, \end{equation} \tag{ 23 }$$

where r_ij is the distance between two atoms of type i and j. The parameters A_ij, ρ_ij and C_ij for each atom pair type are taken from the force field of Tsuneyuki [8]. To obtain the values of $B_0^{\prime }$ at different temperatures, we used the third-order Birch–Murnaghan equation of state (EoS) [9] to fit to the isotherms of these structures from the MD, as shown in figure 4 for T = 1 K. Fitting to the isotherms at other temperatures from the simulations and testing the calculated results against the available experiments can be found in [4].

Zoom In Zoom Out Reset image size

According to the simple model in section 2, T₀ and T₁ in the phenomenological expression are akin to the (rotation) modes that lead to NTE of the system. To have an idea of how these modes exist in the zeolite structures we have studied, the density of states (DoS) and the Grüneisen parameters were calculated by harmonic lattice dynamics using GULP [10] with the same potential model as in the MD. As shown in figure 5, where the DoS of each structure is coloured according to the values of Grüneisen parameters, all the modes with significant negative Grüneisen parameters hence most relevant to NTE of the structure are below 3.0 THz (∼150 K). Thus, to reduce the error introduced by the deviation of the classical simulation from the real quantum picture at low temperatures [11], we will only use the data points of $B^{\prime }_0$ at temperature beyond 200 K when any mode with significant negative Grüneisen parameter has been well excited.

Zoom In Zoom Out Reset image size

The zero-temperature value of $B^{\prime }_0$, i.e. the value of $\beta$, was obtained from the harmonic lattice calculations. Figure 6 shows the fitted temperature dependence of $B^{\prime }_0$ for each zeolite using the phenomenological expression (equation (1)). The fitted values of variables T₀ and T₁ are listed in table 1. Although the data points below 200 K were not used in the fitting, they show consistent trend with the curves.

Zoom In Zoom Out Reset image size

Table 1.  Values of T₀ and T₁ for the zeolite structures from fitting to the $B^{\prime }_0$ data using the phenomenological expression of equation (1). The value of $B^{\prime }_0$ at zero temperature, i.e. $\beta$, for each structure was calculated using harmonic lattice dynamics.

Zeolite	$\beta$	T0 (K)	T1 (K)	T1 − 2T0 (K)
FAU	2.4	89(16)	493(98)	314(114)
KFI	2.6	84(12)	604(113)	437(125)
LTA	3.0	17(2)	260(25)	225(27)
TSC	0.5	12(2)	384(59)	361(61)

From figure 5, we found that, among all the modes that have negative Grüneisen parameters, those with frequencies no larger than 1.5 THz (indicated by the dotted line) account for 70% of the total number of modes and these modes also constituent 70% NTE of the structure in terms of the value of the overall Grüneisen parameter. It is clear that the fitted T₀ of the zeolites agrees with this characteristic frequency (1.5 THz ≈ 75 K). The values of T₁ − 2T₀ of the zeolites indicate that LTA should have larger anharmonicity than TSC and FAU should have larger anharmonicity compared to KFI.

Because $B^{\prime }_0$ is a second-order quantity, it is very sensitive to errors in isotherm data, which gives us a challenge to directly measure its temperature dependence accurately, especially for low-temperature points [3]. However, by using the proposed phenomenological expression of equation (1), one is able to have access to the temperature dependency of $B^{\prime }_0$ of NTE materials at low temperatures, with the zero-temperature value easily calculated from a reliable theoretic model and the data points at high temperatures where the experimental error becomes acceptable.