Large thermal anisotropy in monoclinic niobium trisulfide: A thermal expansion tensor study

TL;DR

This paper introduces a method based on the Gruneisen formalism to calculate the thermal expansion tensor for low-symmetry crystals, and applies it to reveal high anisotropy in monoclinic niobium trisulfide, challenging common beliefs about negative TEC components.

Contribution

The study develops a first-principles approach to compute the TEC tensor in monoclinic crystals, enabling analysis of anisotropic thermal expansion in low-symmetry systems.

Findings

Unveiled high TEC anisotropies in monoclinic NbS3

Negative TEC component attributed to mechanical properties, not anharmonicity

Method can be generalized to triclinic systems with multiple off-diagonal elements

Abstract

We present a method based on the Gruneisen formalism to calculate the thermal expansion coefficient (TEC) tensor that is applicable to any crystal system, where the number of phonon calculations associated with different deformations scales linearly with the number of lattice parameters. Compared to simple high-symmetry systems such as cubic or hexagonal systems, a proper consideration of low-symmetry systems such as monoclinic or triclinic crystals demands a clear distinction between the TEC tensor and the lattice-parameter TECs along the crystallographic directions. The latter is more complicated and it involves integrating the equations of motion for the primitive lattice vectors, with input from the TEC tensor. A first-principles study of the TEC is carried out for the first time on a monoclinc crystal, where we unveil high TEC anisotropies in a recently reported monoclinic phase of niobium trisulfide (NbS3) crystal with a relatively large primitive cell (32 atoms per cell) using density-functional theory. We find the occurrence of a negative TEC tensor component is largely due to the mechanical property rather than the anharmonic effect, contrary to the common belief. Our theoretical treatment of the monoclinic system with a single off-diagonal tensor element could be routinely generalized to any crystal system, including the lowest-symmetry triclinic system with three off-diagonal tensor elements.

Figures

Tables

Table 1. Table 1: Thirteen deformation types for the calculation of all thirteen elastic constant for a monoclinic crystal. The strain vector ϵ bold-italic-ϵ {\boldsymbol{\epsilon}} defines a deformation applied to the crystal in equilibrium according to Eq.  18 . The third column describes the increase of the total energy Δ ​ E Δ 𝐸 \Delta E as a quadratic dependence on ϵ italic-ϵ \epsilon .

Deformation	Strain	$\Delta E/\Omega$ to $O\left({ϵ}^2\right)$
type	vector $ϵ$
1	$(ϵ,0,0,0,0,0)$	$\frac{1}{2}{c}_1{ϵ}^2$
2	$(0,ϵ,0,0,0,0)$	$\frac{1}{2}{c}_5{ϵ}^2$
3	$(0,0,ϵ,0,0,0)$	$\frac{1}{2}{c}_8{ϵ}^2$
4	$(0,0,0,ϵ,0,0)$	$\frac{1}{2}{c}_{10}{ϵ}^2$
5	$(0,0,0,0,ϵ,0)$	$\frac{1}{2}{c}_{12}{ϵ}^2$
6	$(0,0,0,0,0,ϵ)$	$\frac{1}{2}{c}_{13}{ϵ}^2$
7	$(ϵ,ϵ,0,0,0,0)$	$\frac{1}{2}\left(c_1+2c_2+c_5\right){ϵ}^2$
8	$(ϵ,0,ϵ,0,0,0)$	$\frac{1}{2}\left(c_1+2c_3+c_8\right){ϵ}^2$
9	$(0,ϵ,ϵ,0,0,0)$	$\frac{1}{2}\left(c_5+2c_6+c_8\right){ϵ}^2$
10	$(ϵ,0,0,0,ϵ,0)$	$\frac{1}{2}\left(c_1+2c_4+c_{12}\right){ϵ}^2$
11	$(0,ϵ,0,0,ϵ,0)$	$\frac{1}{2}\left(c_5+2c_7+c_{12}\right){ϵ}^2$
12	$(0,0,ϵ,0,ϵ,0)$	$\frac{1}{2}\left(c_8+2c_9+c_{12}\right){ϵ}^2$
13	$(0,0,0,ϵ,0,ϵ)$	$\frac{1}{2}\left(c_{10}+2c_{11}+c_{13}\right){ϵ}^2$

Equations

$$\frac{d{ϵ}_{ij}}{dT}={\alpha }_{ij}(T)$$

$$\left(\begin{array}{ccc}{\alpha }_{11}& 0& {\alpha }_{31}\\ 0& {\alpha }_{22}& 0\\ {\alpha }_{31}& 0& {\alpha }_{33}\end{array}\right)$$

$${\alpha }_{\ell }(T)=\frac{1}{\ell (T)}\frac{d\ell (T)}{dT}$$

$$\mathbf{a}=a\mathbf{i},$$

$$\mathbf{b}=b\mathbf{j},$$

$$\mathbf{c}=c\cos \beta \mathbf{i}+c\sin \beta \mathbf{k}.$$

$${\mathbf{a}}^{\prime }=\mathbf{a}=a\mathbf{i},$$

$${\mathbf{b}}^{\prime }=\mathbf{b}=b\mathbf{j},$$

$${\mathbf{c}}^{\prime }=\mathbf{c}-\mathbf{a}=c^{\prime }\cos {\beta }^{\prime }\mathbf{i}+c^{\prime }\sin {\beta }^{\prime }\mathbf{k},$$

$$\alpha_{v}=\left|\begin{array}[]{ccc}\ 1+\alpha_{11}&0&\alpha_{31}\\ 0&1+\alpha_{22}&0\\ \alpha_{31}&0&1+\alpha_{33}\end{array}\right|-1$$

$${\alpha }_v={\alpha }_a+{\alpha }_b+{\alpha }_c+\frac{\beta \cos \beta }{\sin \beta }{\alpha }_{\beta }$$

$${\alpha }_{ij}(T)=\frac{1}{\Omega }\sum _{k=1}^3\sum _{l=1}^3{S}_{ijkl}{I}_{kl}(T)$$

$$I_{ij}(T)=\frac{\Omega }{{\left(2\pi \right)}^3}\sum _{\lambda }{\int }_{\mathrm{BZ}}d\mathbf{q}{\gamma }_{ij,\lambda \mathbf{q}}c({\nu }_{\lambda \mathbf{q}},T)$$

$${\gamma }_{ij,\lambda \mathbf{q}}=-\frac{1}{{\nu }_{\lambda \mathbf{q}}}\frac{\partial {\nu }_{\lambda \mathbf{q}}}{\partial {ϵ}_{ij}}$$

$$I_{ij}(T)={\int }_{{\nu }_{min}}^{{\nu }_{max}}d\nu {\Gamma }_{ij}(\nu )c(\nu ,T)$$

$${\Gamma }_{ij}(\nu )=\frac{\Omega }{{\left(2\pi \right)}^3}\sum _{\lambda }{\int }_{\mathrm{BZ}}d\mathbf{q}{\gamma }_{ij,\lambda \mathbf{q}}\delta (\nu -{\nu }_{\lambda \mathbf{q}})$$

$${\alpha }_i(T)=\frac{1}{\Omega }{S}_{ij}{I}_j(T)$$

$$\left(\begin{array}{ccc}{\alpha }_{11}(T)& {\alpha }_{12}(T)& {\alpha }_{31}(T)\\ {\alpha }_{12}(T)& {\alpha }_{22}(T)& {\alpha }_{23}(T)\\ {\alpha }_{31}(T)& {\alpha }_{23}(T)& {\alpha }_{33}(T)\end{array}\right)\to \left(\begin{array}{ccc}{\alpha }_1(T)& \frac{1}{2}{\alpha }_6(T)& \frac{1}{2}{\alpha }_5(T)\\ \frac{1}{2}{\alpha }_6(T)& {\alpha }_2(T)& \frac{1}{2}{\alpha }_4(T)\\ \frac{1}{2}{\alpha }_5(T)& \frac{1}{2}{\alpha }_4(T)& {\alpha }_3(T)\end{array}\right)$$

$$\left(\begin{array}{ccc}{I}_{11}(T)& I_{12}(T)& I_{31}(T)\\ I_{12}(T)& I_{22}(T)& I_{23}(T)\\ I_{31}(T)& I_{23}(T)& I_{33}(T)\end{array}\right)\to \left(\begin{array}{ccc}{I}_1(T)& I_6(T)& I_5(T)\\ I_6(T)& I_2(T)& I_4(T)\\ I_5(T)& I_4(T)& I_3(T)\end{array}\right)$$

$$C=\left(\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& C_{15}& 0\\ C_{12}& C_{22}& C_{23}& 0& C_{25}& 0\\ C_{13}& C_{32}& C_{33}& 0& C_{35}& 0\\ 0& 0& 0& C_{44}& 0& C_{46}\\ C_{15}& C_{25}& C_{35}& 0& C_{56}& 0\\ 0& 0& 0& C_{46}& 0& C_{66}\end{array}\right)$$

$$A=\left(\begin{array}{ccc}a& 0& c\cos \beta \\ 0& b& 0\\ 0& 0& c\sin \beta \end{array}\right)$$

$$E=\left(\begin{array}{ccc}{ϵ}_{11}& {ϵ}_{12}& {ϵ}_{31}\\ {ϵ}_{12}& {ϵ}_{22}& {ϵ}_{23}\\ {ϵ}_{31}& {ϵ}_{23}& {ϵ}_{33}\end{array}\right)=\left(\begin{array}{ccc}{ϵ}_1& \frac{1}{2}{ϵ}_6& \frac{1}{2}{ϵ}_5\\ \frac{1}{2}{ϵ}_6& {ϵ}_2& \frac{1}{2}{ϵ}_4\\ \frac{1}{2}{ϵ}_5& \frac{1}{2}{ϵ}_4& {ϵ}_3\end{array}\right)$$

$${\gamma }_{31,\lambda \mathbf{q}}=-\frac{1}{{\nu }_{\lambda \mathbf{q}}}\frac{\partial {\nu }_{\lambda \mathbf{q}}}{\partial {ϵ}_{31}}=-\frac{2}{{\nu }_{\lambda \mathbf{q}}}\frac{\partial {\nu }_{\lambda \mathbf{q}}}{\partial {ϵ}_5}$$

$${\gamma }_{5,\lambda \mathbf{q}}\equiv -\frac{2}{{\nu }_{\lambda \mathbf{q}}}\frac{\partial {\nu }_{\lambda \mathbf{q}}}{\partial {ϵ}_5}$$

$${\gamma }_{5,\lambda \mathbf{q}}^b=-\frac{1}{{\nu }_{\lambda \mathbf{q}}}\frac{\partial {\nu }_{\lambda \mathbf{q}}}{\partial {ϵ}_5}$$

$${\Gamma }_5(\nu )={\Gamma }_{31}(\nu )=\frac{\Omega }{{\left(2\pi \right)}^3}\sum _{\lambda }{\int }_{\mathrm{BZ}}d\mathbf{q}{\gamma }_{5,\lambda \mathbf{q}}\delta (\nu -{\nu }_{\lambda \mathbf{q}})$$

$$I_5(T)=I_{31}(T)=\frac{\Omega }{{\left(2\pi \right)}^3}\sum _{\lambda }{\int }_{\mathrm{BZ}}d\mathbf{q}{\gamma }_{5,\lambda \mathbf{q}}c({\nu }_{\lambda \mathbf{q}},T)$$

$$C=\left(\begin{array}{cccccc}{c}_1& c_2& c_3& 0& c_4& 0\\ c_2& c_5& c_6& 0& c_7& 0\\ c_3& c_6& c_8& 0& c_9& 0\\ 0& 0& 0& c_{10}& 0& c_{11}\\ c_4& c_7& c_9& 0& c_{12}& 0\\ 0& 0& 0& c_{11}& 0& c_{13}\end{array}\right)$$

$$\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\\ {\alpha }_5\end{array}\right)=\frac{1}{\Omega }{\left(\begin{array}{cccc}{c}_1& c_2& c_3& c_4\\ c_2& c_5& c_6& c_7\\ c_3& c_6& c_8& c_9\\ c_4& c_7& c_9& c_{12}\end{array}\right)}^{-1}\left(\begin{array}{c}{I}_1\\ I_2\\ I_3\\ I_5\end{array}\right)$$

$$\tilde{\alpha }(T)=\left(\begin{array}{ccc}{\alpha }_1(T)& \frac{1}{2}{\alpha }_6(T)& \frac{1}{2}{\alpha }_5(T)\\ \frac{1}{2}{\alpha }_6(T)& {\alpha }_2(T)& \frac{1}{2}{\alpha }_4(T)\\ \frac{1}{2}{\alpha }_5(T)& \frac{1}{2}{\alpha }_4(T)& {\alpha }_3(T)\end{array}\right)$$

$$\left(\begin{array}{ccc}{a}_1(T^{\prime })& 0& c_1(T^{\prime })\\ 0& b_2(T^{\prime })& 0\\ a_3(T^{\prime })& 0& c_3(T^{\prime })\end{array}\right)=D\left(\begin{array}{ccc}{a}_1(T)& 0& c_1(T)\\ 0& b_2(T)& 0\\ a_3(T)& 0& c_3(T)\end{array}\right)$$

$$D=I+\tilde{\alpha }(T)\Delta T$$

Full text

Thermal expansion, anharmonic effects, monoclinic crystals, elastic constants, density-functional calculations, phonon calculations, niobium trisulfide, Grüneisen parameter

I Introduction

Thermal expansion of materials is a direct manifestation of the anharmonic effects inherent in a solid. A good understanding of the phonon dynamics affected by temperature-induced crystal deformations is often required for improved device applicationsGoodwin et al. (2008); Lin et al. (2011); Kim et al. (2012). So far ab initio studies of anharmonic effects in high-symmetry crystals such as cubic and hexagonal systems have been carried out via the quasiharmonic approximation (QHA) Mounet and Marzari (2005); Stoffel et al. (2015). The application of QHA on low-symmetry crystals is extremely challenging because many sets of time-consuming phonon calculations are required to search for the free energy minimum in a much larger lattice-parameter space, which may explain why, to the best of our knowledge, no studies have been carried out on, say, a monoclinic system. Even if a direct free energy minimization via QHA could be carried out, it may be still difficult to understand the underlying physics without interpreting the phonon frequency shifts due to a crystal deformation, phonon density of states, mechanical properties such as elastic constants and elastic compliances, and heat capacities. Since many technologically important materials are also found in low crystal symmetries, e.g., hybrid perovskitesCalabrese et al. (1991); Du et al. (2017); Giovanni et al. (2018) and ceramicsZhao and Vanderbilt (2002); Haggerty et al. (2014), it is desirable to be able to address the anharmonicity in these materials through other means. Recently, the Grüneisen approach has between applied to high-symmetry crystal systems and good agreements between theoretical predictions and experimental results Ding and Xiao (2015); Liu et al. (2017); Gan and Lee (2018). In this work, we extend the utility of the Grüneisen formalismGrüneisen (1926); Schelling and Keblinski (2003) to handle low-symmetry crystals. As will be shown later in this paper, dealing with low-symmetry crystals requires a clear distinction to be made between the thermal expansion coefficient (TEC) tensor and the lattice-parameter TECs. These two concepts, unfortunately, coincide with each other in high-symmetry crystals such as the cubic crystal. Other methods to calculate TECs are the vibrational self-consistent-field methodMonserrat et al. (2013) and the self-consistent quasiharmonic approximationHuang et al. (2016). We should note that the knowledge of the third-order interatomic force constantsBarron and Klein (1974); Shiomi et al. (2011) can be used directly to calculate the Grüneisen parameters for the evaluation of TECs. Another recent work Lee and Gan (2017) outlines a reverse process where numerically computed Grüneisen parameters could be used to extract the third-order interatomic force constants for the calculations of thermal conductivity.

The monoclinic NbS3-IVBloodgood et al. (2018) (with four lattice parameters) belongs to a class of semiconducting transition metal chalcogenides. Transition metal chalcogenides have received attention due to their potential to replace traditional semiconductor technologies. They have been investigated for their electronic and optical properties, and their potential use in electronic devicesRibeiro et al. (2018). The extensive interest in metal chalcogenidesZhao et al. (2011); Chong et al. (2014) arises from the low-dimensionality of the structure. As a prototypical example of metal trichalcogenides (MX3), niobium trisulfide (NbS3) comprises of [NbS6] units, where the sulfur atoms are arranged in a triangular prism around a central Nb atomvan Smaalen (1988). Neighboring trigonal prisms share a triangle face, thus fusing individual units to form an infinite chain along the $b$ axis. The chains are held together by additional covalent bonds to form bilayer sheets which interact via van der Waals forces. Experimentally, NbS3 has been found to crystallize in six polymorphs: I, II, III, IV, V and a high-pressure phase. The first polymorph of NbS3, NbS3-I, first identified using single crystal x-ray diffractionJellinek (1960); Rijnsdorp and Jellinek (1978), belongs to space group $P\bar{1}$. Using electron diffraction studies, NbS3-II was the second polymorph to be identified, although the structure remains unknown to date. Both NbS3-II and NbS3-III exhibit charge density wave (CDW) transitionsZettl et al. (1982); Zybtsev et al. (2017). Recently, Bloodgood et al.Bloodgood et al. (2018) reported the synthesis of new monoclinic polymorphs of NbS3, where structural information was obtained via single crystal x-ray diffraction, electron microscopy, and Raman spectroscopy studiesBloodgood et al. (2018). The new polymorphs, NbS3-IV and NbS3-V, were found to crystallize in the space groups $P12_{1}/c1$ and $P12_{1}/m1$, respectively. Presently, first-principles studies of the anharmonic effects that account for thermal expansion are lacking for these compounds.

The outline of this paper is as follows: Sec. II.1 discusses the concepts of a TEC tensor and lattice-parameter TECs and the distinction between them. Sec. II.2 presents the Grüneisen formalism for the computation of the TEC tensors. For efficient implementation and a good physical interpretation, Sec. II.3 shows how the tensorial expression can be turned into a simple matrix expression. Sec. II.4 discusses the deformations that are to be applied to the crystal for the calculations of elastic constants and Grüneisen parameters. The procedure for elastic constant calculations is outlined in Sec. II.5 where we observe that all required elastic constants for the TEC calculations can be obtained from symmetry-preserving deformations. In Sec. II.6, we describe how the lattice-parameter TECs may be obtained from a TEC tensor by solving a set of coupled differential equations. The reverse of this process is presented in Appendix A. The results are presented in Sec. III. Finally Sec. IV contains the conclusions.

II Methodology

II.1 TEC tensor versus lattice-parameter TECs

Thermal expansion of a crystal can generally be described by the thermal expansion coefficient (TEC) tensor with components $\alpha_{ij}(T)$ that depend on temperature $T$. When the temperature changes, the crystal will change its shape according to a strain tensor $\epsilon_{ij}$ induced by the TEC tensor where

[TABLE]

In three dimensions, the most general $\alpha_{ij}(T)$ is a $3\times 3$ symmetric tensor with six independent parameters. For a cubic crystal, it turns out that $\alpha_{11}(T)=\alpha_{22}(T)=\alpha_{33}(T)$ while all crossed terms $\alpha_{ij}(T)=0$, $i\neq j$, which means that only one independent parameter is needed to describe the TEC. Also it does not matter how the cubic crystal is oriented in space since the expansion is isotropic. Specifically, the TEC tensor components $\alpha_{ii}(T)$ are all equal for $i=1,2,3$, and they coincide with the meaning of the $a$ lattice-parameter TEC, denoted by $\alpha_{a}(T)=\frac{1}{a(T)}{\frac{da(T)}{dT}}$ for the cubic crystal where $a$ is the lattice constant. In another example, the hexagonal crystal is described by two independent TEC tensor components $\alpha_{a}(T)=\alpha_{11}(T)=\alpha_{22}(T)$, and $\alpha_{c}(T)=\alpha_{33}(T)$ where we must tacitly assume that the $c$ axis is pointing in the $z$ Cartesian direction. However, for a low-symmetry monoclinic crystal (with lattice parameters $a\neq b\neq c$ and $\beta\neq 90^{\circ}$), the matter becomes slightly more complicated. It is usually convenient to choose a crystal orientation where the $b$ lattice parameter is parallel to the $y$ Cartesian direction (the so-called ‘unique axis $b$’ in Ref. [Hahn, Ed.]). Symmetry considerations and Neumann’s principle dictates that the TEC tensor of a monoclinic crystal is of the following formSchlenker et al. (1975); Nye (1985)

[TABLE]

Here, there are four independent components of the TEC tensor. Similar to the basic definition of the lattice-parameter TEC $\alpha_{a}(T)$ for a cubic crystal, we define a lattice-parameter TEC of a monoclinic crystal according to

[TABLE]

where $\ell$ equals to one of the monoclinic lattice parameters $a$, $b$, $c$, and $\beta$. Since the $b$ axis is unique by choice, the $\alpha_{22}(T)$ component is decoupled from the other three components and coincides with the $b$ lattice-parameter TEC, where $\alpha_{b}(T)=\alpha_{22}(T)$. With a ‘unique axis b’ choice, it is still necessary to clearly state the crystal orientation before reporting the other three TEC tensor components $\alpha_{11}(T)$, $\alpha_{33}(T)$, and $\alpha_{31}(T)$. Unless the $a$ lattice parameter happens to be aligned along the $x$ Cartesian axis, $\alpha_{11}(T)$ is not same as the lattice-parameter TEC for $a$. It is then apparent that the temperature dependence of $a$, $c$, and $\beta$ collectively depend on $\alpha_{11}(T)$, $\alpha_{33}(T)$, and $\alpha_{31}(T)$, and vice versa. In Sec. II.6 and Appendix A we derive the relations between these two TEC descriptions. Throughout this work, we assume the orientation of a monoclinic crystal in equilibrium at $T=0$ K is as follows:

[TABLE]

The Cartesian unit vectors are ${\boldsymbol{i}}$, ${\boldsymbol{j}}$, and ${\boldsymbol{k}}$. It is important to note that for $T>0$ K, ${\boldsymbol{a}}$ may have a nonzero component along ${\boldsymbol{k}}$ due to presence of the off-diagonal $\alpha_{31}(T)$ component.

We shall now comment on the nonuniqueness of the choice of $a$, $c$, and $\beta$ of a monoclinic primitive cell, even within a ‘unique axis $b$’ choice. With the specific choice adopted in Eq. 4, we may derive from it a second set of primitive lattice vectors such as

[TABLE]

This will in general have the consequence that $\alpha_{c^{\prime}}\neq\alpha_{c}$ and $\alpha_{\beta^{\prime}}\neq\alpha_{\beta}$. We shall see the effect of such a choice on the results in Section III.

We may calculate the volumetric TEC $\alpha_{v}$ (with suppressed $T$ dependence for clarity) from

[TABLE]

The same $\alpha_{v}$ can also be calculated from the lattice-parameter TECs where

[TABLE]

This can be derived from the expression for the volume of the primitive cell $\Omega$ where $\Omega=abc\sin\beta$.

II.2 Grüneisen formalism

The TEC tensor components $\alpha_{ij}(T)$ for any crystal at temperature $T$ may be calculated from first-principles via a Grüneisen formalismGrüneisen (1926); Schelling and Keblinski (2003) expressed in a tensorial form:

[TABLE]

The elastic compliance tensor and equilibrium volume of the primitive cell are denoted by $S_{ijkl}$ and $\Omega$, respectively.

Central to the Grüneisen formalism are the integrated quantities $I_{ij}(T)$ given by

[TABLE]

Here $\nu_{\lambda{\boldsymbol{q}}}$ denotes the frequency of a phonon mode with mode index $\lambda$ and wave-vector ${\boldsymbol{q}}$. The heat capacity of a phonon mode of frequency $\nu$ is given by $c(\nu,T)=k_{{\scriptscriptstyle B}}(r/\sinh r)^{2}$, with $r=h\nu/2k_{{\scriptscriptstyle B}}T$, where $h$ and $k_{{\scriptscriptstyle B}}$ are the Planck and Boltzmann constants, respectively. The bandstructure-like deformation-dependent Grüneisen parameter is defined by

[TABLE]

which measures the fractional change in phonon mode frequency with respect to a small strain $\epsilon_{ij}$ applied to the crystal. Phonon frequencies may be calculated accurately either from a density-functional perturbation theory (DFPT) approachBaroni et al. (2001) or a supercell force-constant methodKresse et al. (1995); Ackland et al. (1997); Alfè (2009). The Grüneisen parameters are calculated based on Eq. 10 from phonon frequency shifts using a central-difference scheme with the help of changes in the dynamical matrices and a perturbation theoryGan et al. (2015). In this work, we apply strains of $\epsilon=\pm 0.01$ for all deformations to determine the Grüneisen parameters.

The integrated quantity $I_{ij}(T)$ in Eq. 9 is seen as an integral over the Brillouin zone (BZ) of the product of the mode-dependent heat capacities and mode-dependent Grüneisen parameters. It can also be written as an integral over the phonon frequencies

[TABLE]

where $\nu_{min}$ ($\nu_{max}$) is the minimum (maximum) frequency in the phonon spectrum. Here we have introduced a quantity called the density of the phonon states weighted by the Grüneisen parameters

[TABLE]

With the definition in Eq. 12, we turn the direct summation in Eq. 9 into an integral given in Eq. 11 involving the two quantities $c(\nu,T)$, which has a form that is easily understood, and $\Gamma_{ij}(\nu)$, which captures the net effect of anharmonicities of all phonon modes with the same frequency $\nu$.

II.3 Compact matrix expression

We find it convenient to use the analogy of the strain-stress relationshipNye (1985) $\epsilon_{i}=S_{ij}\sigma_{j}$ to cast Eq. 8 into

[TABLE]

with the prescription

[TABLE]

and

[TABLE]

Now, we need to discuss how to obtain $S_{ij}$ from the elastic constants. In general, the $6\times 6$ symmetric elastic constant matrix $C$ with matrix elements $C_{ij}$ consists of $21$ independent matrix elementsNye (1985).

For a monoclinic crystal with ‘a unique axis $b$’ choice, the matrix $C$ with $13$ independent matrix elements is

[TABLE]

We will give more details on the calculation of elastic constants in Sec. II.5. The elastic compliance matrix $S$ with matrix elements $S_{ij}$ is simply the inverse of $C$, i.e., $S=C^{-1}$.

II.4 Crystal deformation

We describe how to perform a deformation to the monoclinic crystal with lattice parameters $a$, $b$, $c$, and $\beta$. Deformations are required for the elastic constants and the mode-dependent Grüneisen parameters according to Eq. 10. At [math] K, we use a matrix $A$ obtained from the primitive lattice vectors (see Eq. 4)

[TABLE]

Next we let the deformation matrixNye (1985) be

[TABLE]

so that the primitive lattice vectors ${\boldsymbol{a}}$, ${\boldsymbol{b}}$, and ${\boldsymbol{c}}$ after a deformation are described by $A^{\prime}=(I+E)A$, where $I$ is the identity matrix. The full $3\times 3$ matrix $E$ can be equivalently characterized by a strain vector ${\boldsymbol{\epsilon}}=(\epsilon_{1},\epsilon_{2},\epsilon_{3},\epsilon_{4},\epsilon_{5},\epsilon_{6})$. In this notation ${\boldsymbol{\epsilon}}=(\epsilon_{1},0,0,0,0,0)$ describes a uniaxial strain of $\epsilon_{1}$ in the $x$ Cartesian direction, while ${\boldsymbol{\epsilon}}=(0,0,0,0,e_{5},0)$ describes a shear strain of $\epsilon_{31}={\frac{1}{2}}\epsilon_{5}$ in the $x$ and $z$ Cartesian directions. Since we use the strain vector ${\boldsymbol{\epsilon}}$ with single-index entries to describe a deformation, care should be exercised when dealing with the off-diagonal term

[TABLE]

We find it is necessary to define

[TABLE]

so to make $\gamma_{5,\lambda{\boldsymbol{q}}}=\gamma_{31,\lambda{\boldsymbol{q}}}$. Note that this is different from the ‘bare’ Grüneisen parameter (which holds for the diagonal elements $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$) where

[TABLE]

and hence $\gamma_{5,\lambda{\boldsymbol{q}}}=2\gamma^{b}_{5,\lambda{\boldsymbol{q}}}$. We then have from Eq. 12

[TABLE]

and from Eq. 15

[TABLE]

II.5 Elastic constant calculation

In this subsection, we describe how to determine the $13$ elastic constants $C_{ij}$ for a monoclinic crystal. We follow the approach in Refs. Beckstein et al., 2001; Dal Corso, 2016 due to its simplicity.

To further simplify notations we express the elastic constants in Eq. 16 as

[TABLE]

These constants can be obtained by performing thirteen independent deformation types according to Table 1. For each deformation type, a set of strains $\epsilon$ ranging from $-0.01$ to $0.01$ is used to obtain an energy $\Delta E/\Omega$ versus strain $\epsilon$ curve. A quadratic least-square fit is used to obtain the curvature $k$ for the parabola. Each curvature is in general a linear combination of a few elastic constants, e.g., for the deformation type $7$, since $\frac{\Delta E}{\Omega}={\frac{1}{2}}(c_{1}+2c_{2}+c_{5})\epsilon^{2}$ (see Table 1), $k_{7}={\frac{1}{2}}(c_{1}+2c_{2}+c_{5})$. The values of $c_{i}$, $i=1,2,\ldots,13$ can be obtained from $k_{i}$ by a matrix inversion.

We note that of all the thirteen deformations, ten of them preserve the space group of the crystal after a deformation is applied. These are the deformation types $1$, $2$, $3$, $5$, $7$, $8$, $9$, $10$, $11$, and $12$. Since the effect of temperature does not change the crystal type (i.e., the crystal remains monoclinic upon heating or cooling), we expect that $I_{4}(T)=I_{6}(T)=0$, and Eq. 13 is thus reduced to

[TABLE]

For a monoclinic crystal, we observe that all ten deformations needed for the TEC calculation can be chosen to be strictly space-group preserving. This is important for an accurate determination of elastic constants because we can maintain the same number of ${\boldsymbol{k}}$ points throughout all self-consistent density-functional theory calculations, which minimizes possible total-energy fluctuation. The existence of a minimum yet sufficient number of symmetry-preserving deformation types to calculate just all relevant elastic constants for the sole purpose of TEC determination is also observed in other crystal types: cubic with one uniform deformation, orthorhombicGan et al. (2015) with six deformations, hexagonal and trigonal with three deformationsGan and Liu (2016); Gan and Lee (2018), and of course, the triclinic system with the full 21 deformations for all 21 elastic constants (the space group of a triclinic crystal is either $P1$ or $P{\overline{1}}$).

II.6 Calculation of the lattice-parameter TECs

Having obtained the TEC tensor $\alpha_{i}(T)$, we still need a recipe to find the lattice-parameter TECs as defined in Eq. 3. Since $b$ is not coupled to $a$, $c$, and $\beta$, we have immediately $\alpha_{b}(T)=\alpha_{2}(T)$. To find the other three lattice-parameter TECs $\alpha_{a}$, $\alpha_{c}$, and $\alpha_{\beta}$, we need to know how the lattice parameters $a,c,\beta$ vary with $T$ and then use a simple numerical difference scheme (such as forward difference) to extract $\alpha_{\ell}$ according to Eq. 3. Here we need to derive a set of differential equations for numerical integration. We assume the primitive lattice vectors at temperature $T$ are ${\boldsymbol{a}}(T)=a_{1}(T){\boldsymbol{i}}+a_{3}(T){\boldsymbol{k}}$, ${\boldsymbol{b}}(T)=b_{2}(T){\boldsymbol{j}}$, and ${\boldsymbol{c}}(T)=c_{1}(T){\boldsymbol{i}}+c_{3}(T){\boldsymbol{k}}$. At $T=0$ K, we have the initial conditions ${\boldsymbol{a}}(0)=(a,0,0)$, ${\boldsymbol{b}}(0)=(0,b,0)$, ${\boldsymbol{c}}(0)=(c\cos\beta,0,c\sin\beta)$. From Eq. 14 we let

[TABLE]

When $T$ changes to $T^{\prime}=T+\Delta T$, we have

[TABLE]

where the deformation matrix $D$ is given by

[TABLE]

As $\Delta T$ goes to zero, we obtain five differential equations:

[TABLE]

Since $b(T)=|{\boldsymbol{b}}(T)|=b_{2}(T)$, we have immediately $\alpha_{b}(T)=\alpha_{2}(T)$ from Eq. 31, which is expected since the primitive lattice vector ${\boldsymbol{b}}$ is not coupled to ${\boldsymbol{a}}$ and ${\boldsymbol{c}}$. With the initial conditions $a_{1}(0)=a$, $a_{3}(0)=0$, $c_{1}(0)=c\cos\beta$, and $c_{3}(0)=c\sin\beta$, we use a robust fourth-order Runge-Kutta integration scheme to numerically integrate the differential equations. From $a(T)=\sqrt{a_{1}(T)^{2}+a_{3}(T)^{2}}$, and $c(T)=\sqrt{c_{1}(T)^{2}+c_{3}(T)^{2}}$, we obtain the lattice-parameter TECs $\alpha_{a}$ and $\beta_{c}$, and $\alpha_{\beta}$. Now we have outlined a procedure to obtain the lattice-parameter TECs from the TEC tensor. In Appendix A, we outline a procedure to reverse the process, which may be more applicable in experimental settings since lattice-parameter TECs are usually measured first.

III Results and discussions

We carry out density-functional theory (DFT) calculations within the local-density approximation, with projector augmented-wave (PAW) pseudopotentials as implemented in the Vienna Ab initio Simulation Package (VASP)Kresse and Furthmüller (1996) on NbS3-IV. The initial structure of NbS3-IV is taken from Ref. [Bloodgood et al., 2018] where $\beta=90.12^{\circ}$, $a=6.752$, $b=4.974$, and $c=18.132$ Å. This corresponds to the $P12_{1}/c1$ setting of the space group #14. The primitive cell consists of $32$ atoms, $8$ of which are inequivalent. A large cutoff energy of $400$ eV is used throughout this work. A $12\times 12\times 4$ Monkhorst-Pack mesh is used for the ${\boldsymbol{k}}$-point sampling for the relaxation of the lattice parameters of the primitive cell and atomic positions as well as for the calculations of elastic constants. Atomic relaxation is stopped when the maximum force on all atoms is less than $10^{-3}$ eV/Å. To fully optimize the monoclinic lattice parameters, we successively carry out deformations of types $1$, $2$, $3$, and $5$ in Table 1. For each type of deformation, we vary the strain $\epsilon$ from $-0.01$ to $0.01$ to locate the minimum in the energy $E$ versus $\epsilon$ curve. The predicted minimum in the $E$-$\epsilon$ curve allows us to adjust the lattice parameters before a deformation of another type is carried out. At the end of the procedure, we obtain $\beta=89.98^{\circ}$, $a=6.673$, $b=4.870$, and $c=17.837$ Å, which compare favorably with the experimental resultsBloodgood et al. (2018).

We then perform the calculations for the elastic constants as outlined in Sec. II.5. For each type of deformation, we also vary $\epsilon$ from $-0.01$ to $0.01$. The elastic constants $C_{ij}$ of NbS3-IV are $183.64$, $20.49$, $26.32$, $-0.58$, $158.08$, $7.99$, $0.81$, $47.79$, $1.21$, $7.20$, $0.13$, $23.47$, and $33.38$ GPa, for $ij=11$, $12$, $13$, $15$, $22$, $23$, $25$, $33$, $35$, $44$, $46$, $55$, and $66$, respectively. Subsequently, we obtain the elastic compliance matrix $S$ (in $10^{-3}$/GPa) where

[TABLE]

Next we perform the phonon calculations using the supercell force-constant methodGan et al. (2010); Liu et al. (2014). We use a $3\times 3\times 1$ supercell (with $288$ atoms) with a $4\times 4\times 4$ Monkhorst-Pack mesh for the ${\boldsymbol{k}}$-point sampling. A small $0.015$ Å atomic displacement is used to induce forces on all atoms in the supercell. Due to a relatively low symmetry of the space group #14, a total of $48$ independent self-consistent DFT calculations are to be performed for a complete set of phonon calculations.

For the equilibrium structure, we plot the phonon dispersions along a few high-symmetry directions as shown in Fig. 1(a). The phonon spectrum has a highest frequency of $577$ cm*-1* and a phonon gap between $407$ and $519$ cm*-1*. The phonon density of states $\rho(\nu)=\frac{\Omega}{(2\pi)^{3}}\sum_{\lambda}\int_{\rm BZ}d{\boldsymbol{q}}\ \delta(\nu-\nu_{\lambda{\boldsymbol{q}}})$ shown in Fig. 1(b) does not exhibit imaginary modes, implying that NbS3-IV is stable.

To obtain all the necessary Grüneisen parameters for calculation of the TEC tensor, we apply deformations of types $1$, $2$, $3$, and $5$ in Table 1 with $\epsilon=\pm 0.01$. This entails a total of eight sets of phonon calculations, with each set involves $48$ independent self-consistent DFT calculations on a large $3\times 3\times 1$ supercell. From the frequency changes, we obtain the Grüneisen parameters as shown in Fig. 2(a)-(d). We see that Grüneisen parameters with large magnitudes tend to be found near $\Gamma$, where low-frequency phonons are found. We calculate the density of Grüneisen parameters according to $g_{i}(\gamma)=\frac{\Omega}{(2\pi)^{3}}\sum_{\lambda}\int_{\rm BZ}d{\boldsymbol{q}}\ \delta(\gamma-\gamma_{i,\lambda{\boldsymbol{q}}})$ for each $i$, and the right panels of Fig. 2 indicate that most Grüneisen parameters are concentrated in the $-2$ to $2$ range.

We then evaluate $\Gamma_{i}(\nu)$, the phonon density of states weighted by the Grüneisen parameters, for $i=1,2,3,5$. A $30\times 30\times 10$ mesh is used for the ${\boldsymbol{q}}$-point sampling. Fig. 3(b) shows very high anisotropies of the Grüneisen parameters associated with different deformations. There is a gap between $407$ to $519$ cm*-1* in all $\Gamma_{i}(\nu)$, which is similar to that of the phonon density of states shown in Fig. 1(b). The results show that high-frequency modes (around $550$ cm*-1*) are dominated by negative Grüneisen parameters. The dependency of the heat capacity on temperature shown in Fig. 3(a) tells us that these modes contribute to $I_{i}(T)$ (via Eq. 11) only at high temperatures. The deformation due to an $xz$ shear strain gives $\Gamma_{5}(\nu)$ which is only slightly anharmonic relative to that of the diagonal deformations (i.e., $\Gamma_{1}(\nu)$, $\Gamma_{2}(\nu)$, and $\Gamma_{3}(\nu)$). These diagonal deformations are mostly positive, except for $\Gamma_{2}(\nu)$ in the range of $50$ to $75$ cm*-1*.

The integrated quantities $I_{i}(T)$ for $i=1$, $2$, $3$, and $5$ are shown in Fig. 4 may be obtained through either a direct summation in the ${\boldsymbol{q}}$-space according to Eq. 9 or through an integral over the frequencies according to Eq. 11. It is seen that $I_{1}(T)$ is larger than $I_{2}(T)$ and $I_{3}(T)$ at high temperatures above $100$ K, which is consistent with the data in Fig. 3(b). It is seen that $I_{5}(T)$ is indeed small compared to the diagonal quantities $I_{1}(T)$, $I_{2}(T)$, and $I_{3}(T)$ over the entire temperature range, which is again consistent with the small values of $\Gamma_{5}(\nu)$ shown in Fig. 3(b).

The TEC tensor components $\alpha_{i}$, $i=1,2,3,5$ are displayed as solid lines in Fig. 5. It is somewhat surprising that $\alpha_{1}(T)$ is the smallest among all diagonal TEC tensor components, despite the fact that $I_{1}(T)$ is the dominant diagonal component, at least for $T>100$ K as shown in Fig. 4. We can rationalize this result as follows. According to Eq. 13, $\alpha_{1}(T)$ is a linear combination of $I_{i}(T)$, weighted by coefficients taken from the first row of the elastic compliance matrix $S$ in Eq. 34. Since $S_{13}=-3.20\times 10^{-3}$/GPa is negative and large in magnitude compared to $S_{11}=5.97\times 10^{-3}$/GPa, and $I_{1}(T)$ and $I_{3}(T)$ are comparable in magnitudes, this makes $\alpha_{1}(T)$ the smallest among all diagonal $\alpha_{i}(T)$. It should be noted that $\alpha_{1}(T)$ is negative between [math] and $50$ K (shown in the inset of Fig. 5) partly because $I_{3}(T)$ is dominant for temperatures below $100$ K. On the other hand, the dominant value of $S_{33}=22.84\times 10^{-3}$/GPa over $S_{31}$ and $S_{32}$ makes $\alpha_{3}(T)$ the dominant diagonal TEC tensor component. Thus, we conclude that the value of a TEC depends on not only on the anharmonic effects measured by the Grüneisen parameters, as commonly believed, but also on the mechanical properties such as the elastic constants (or equivalently, the elastic compliance matrix). The volumetric TEC of NbS3-IV is $\alpha_{v}=34\times 10^{-6}$ K*-1* at high temperatures.

Next, we turn to the discussion of the lattice-parameter TECs for the crystal in a $P12_{1}/c1$ setting with $\beta=89.98^{\circ}$. Firstly, $\alpha_{b}(T)=\alpha_{2}(T)$ by the choice of the unique axis. From Fig. 5, we find that $\alpha_{a}(T)$ is indistinguishable from $\alpha_{1}(T)$ since the primitive lattice vector ${\boldsymbol{a}}$ is parallel to the $x$ axis at [math] K and it is always almost parallel to the $x$ axis even at higher $T$ due to small values of $\alpha_{5}(T)$. For this setting in which $\beta$ is very close to $90^{\circ}$, the primitive lattice vector ${\boldsymbol{c}}$ remains almost parallel to the $z$-axis at all temperatures, causing $\alpha_{c}(T)$ to be indistinguishable from $\alpha_{3}(T)$. We find that $\alpha_{\beta}(T)$ can be approximated very well (not shown) by $-\frac{2}{\pi}\alpha_{5}(T)$, since $\beta(T)$ is always very close to ${\frac{1}{2}}\pi$ and that $\beta={\frac{1}{2}}\pi-\epsilon_{5}$ (the physical meaning of $\epsilon_{5}$ is the decrease in angle between two elements in the $x$ and $z$ directions, see, e.g., Ref. Nye, 1985), and hence $\alpha_{\beta}=\frac{1}{\beta}{\frac{d\beta}{dT}}=-\frac{1}{\beta}\alpha_{5}(T)\approx-\frac{2}{\pi}\alpha_{5}(T)$.

We discuss next the lattice-parameter TECs of a new primitive cell derived from the $P12_{1}/c1$ setting with $\beta=89.98^{\circ}$, using the transformation rule as shown in Eq. 5. It turns out that this corresponds to a $P12_{1}/n1$ settingNespolo and Aroyo (2016) with $\beta^{\prime}=110.50^{\circ}$, $a^{\prime}=a=6.673$, $b^{\prime}=b=4.870$, $c^{\prime}=19.042$ Å. Since this setting also describes the same unrotated crystal as the $P12_{1}/c1$ setting, the TEC tensor components remain the same. For this new $P12_{1}/n1$ setting, we evolve the monoclinic lattice parameters as previously performed for the $P12_{1}/c1$ setting. It is not surprising to find that $\alpha_{b^{\prime}}(T)=\alpha_{b}(T)=\alpha_{2}(T)$. In addition, $\alpha_{a^{\prime}}(T)=\alpha_{a}(T)$, both of which are essentially the same as $\alpha_{1}(T)$. However, we do find that $\alpha_{c^{\prime}}(T)$ is somewhat different from $\alpha_{3}(T)$ as a result of the primitive lattice vector ${\boldsymbol{c}}^{\prime}$ being oriented away from $z$ axis. Finally, $\alpha_{\beta^{\prime}}(T)$ has an opposite sign compared to $\alpha_{\beta}(T)$ and agrees fortuitously with $\alpha_{5}(T)$.

IV Summary

With an extension of the Grüneisen formalism, we have, for the first time, studied the thermal expansion properties of a monoclinic phase of niobium trisulfide (NbS3) using density-functional theory calculations. Large anisotropies in the thermal expansion coefficient (TEC) tensor components in the $x$, $y$, $z$ directions have been found. A negative TEC tensor component $\alpha_{11}$ is found at low temperatures, the occurence of which is mainly due to the mechanical property characterized by the elastic compliance. This is somewhat surprising since a negative TEC is usually associated with the occurence of many phonon modes with negative Grüneisen parameters. We expect that the Grüneisen parameters generated in this work could be used to provide a good estimate of the thermal conductivity based on a Slack’s expressionToher et al. (2014); Slack (1979); Blanco et al. (2004). The obtained Grüneisen parameters can be further used to extract third-order interatomic force constants for the determination of thermal conductivity by solving the Boltzmann transport equationLi et al. (2014); Lindsay (2016). Our developed method can be applied to any crystal classes, which opens the door to studying many important materials at the level of accurate density-functional theory.

V Acknowledgments

We thank the National Supercomputing Center, Singapore (NSCC) and A*STAR Computational Resource Center, Singapore (ACRC) for computing resources. This work is supported in part by RIE2020 Advanced Manufacturing and Engineering (AME) Programmatic Grant No A1898b0043.

Appendix A Obtaining the TEC tensor from the lattice-parameter TECs

Given the lattice-parameter TECs $\alpha_{a}(T)$ $\alpha_{b}(T)$, $\alpha_{c}(T)$, and $\alpha_{\beta}(T)$, and temperature dependence of the primitive lattice vectors ${\boldsymbol{a}}$, ${\boldsymbol{b}}$, and ${\boldsymbol{c}}$ where ${\boldsymbol{a}}(T)=a_{1}(T){\boldsymbol{i}}+a_{3}(T){\boldsymbol{k}}$, ${\boldsymbol{b}}(T)=b_{2}(T){\boldsymbol{j}}$, and ${\boldsymbol{c}}=c_{1}(T){\boldsymbol{i}}+c_{3}(T){\boldsymbol{k}}$, we seek to obtain the TEC tensor components $\alpha_{i}(T)$, for $i=1,2,3,5$. As mentioned before, we have immediately $\alpha_{2}(T)=\alpha_{b}(T)$. In the following we present a method to calculate $\alpha_{i}(T)$ for $i=1,3,5$.

Consider $a(T)=\sqrt{a_{1}(T)^{2}+a_{3}(T)^{2}}$, we have

[TABLE]

where $y^{\prime}$ means the derivative of $y$ with respect to $T$. We have suppressed the explicit dependence of all quantities on temperature in Eq. 35 for clarity. Replacing $a_{1}^{\prime}$ and $a_{3}^{\prime}$ using Eqs. 29 and 30 we have

[TABLE]

The three unknowns $\alpha_{1}$, $\alpha_{3}$, and $\alpha_{5}$ are seen to appear in a linear equation. Therefore we need two other linear equations to enable a matrix inversion. A second relation is obtained by considering $c(T)=\sqrt{c_{1}(T)^{2}+c_{3}(T)^{2}}$ that gives

[TABLE]

The third relation is obtained by differentiating ${\boldsymbol{a}}\cdot{\boldsymbol{c}}=a_{1}c_{1}+a_{3}c_{3}=ac\cos\beta$ with respect to $T$. Collecting all three relations, we have

[TABLE]

where

[TABLE]

An inversion of Eq. 38 leads to

[TABLE]

Recall all quantities on the right hand side of Eq. 40 depend on $T$, we then have explicit (albeit complicated) expressions for $\alpha_{1}(T)$, $\alpha_{3}(T)$, and $\alpha_{5}(T)$. It may be pointed out that these expressions for TEC tensor components are somewhat different from that found in Ref. 31, presumably due to the fact that they always orient the $c$ lattice parameter so that it is pointing in the $z$ Cartesian direction. Specifically they have $\alpha_{2}(T)=\alpha_{b}(T)$ and $\alpha_{3}(T)=\alpha_{c}(T)$. However, we do not impose such condition but to allow the crystal to freely change its shape under the effect of temperature while keeping the Cartesian coordinates fixed.