Anisotropic Thermal Conduction in Transition Metal Dichalcogenide Nanocomposites with Rough Interfaces

Abstract

We present a theory of thermal conduction in a transition metal dichalcogenide nanocomposite structure with rough interfaces that accounts for the anisotropic conductivities of the host, the insert and the interface regions. The host and insert conductivities are calculated using a semi ab-initio method. The effects of specularity in phonon interface scattering and the thermal boundary resistance is incorporated through linking a phonon wavevector dependent specular scattering parameter to the average height of surface inhomogeneities, and the conductivity of the composite is calculated by employing an extension of a modified effective medium approach. Our work for spherical inserts of WS${}_2$ in MoS${}_2$ predicts that the effects of specular scattering due to surface roughness is more pronounced for inserts smaller than 100 nm, even at volume fractions of the order of 0.05.

1. Introduction

The determination of key parameters of a nanocomposite structure that control its physical properties, such as thermal, electronic and optical characteristics, is essential for future technological advances. The calculation of the bulk properties of nanocomposites where a small fraction of insertions of one material are embedded within a matrix of a different material, perhaps with a boundary layer between the insertion and the matrix, is a long-standing problem. Effective medium approaches (EMAs) are a widely used method of calculating properties (such as the thermal conductivity [1] or the dielectric constant [2]) of composite systems using the properties of their individual constituents [2,3]. In thermal or electrical conductivity calculations, these constituents usually consist of inserts of one material embedded in a matrix of another, along with the effects of the boundary between them. Such a boundary is usually considered to be either an infinitely thin interface layer possesed of a Kapitza or interface resistance [4] (e.g., Ref. [1]) or a finite region with a conductivity of its own (e.g., Ref. [5]). The former approach can be adapted to model a wide variety of systems, for instance cases involving a bimodal distribution of insert particle sizes (e.g., Ref [6]) or systems containing pores (e.g., Ref [7]). However, this form of EMA is only valid for micrometer sized inserts, and requires further modification if it is to be valid for nanocomposite systems. For thermal conductivity calculations, Minnich and Chen [8] proposed a modified effective medum approach (mEMA) which is suitable for nanoscale systems such as the recently synthesised Bi${}_{100-x}$Sb${}_x$/Al${}_2$O${}_3$ nanocomposites [9]. In the mEMA the effects of phonon scattering from the boundaries of spherical inserts are included in the calculation of the matrix and insert thermal conductivities that are input into the closed-form expression for the effective thermal conductivity. This approach has been generalised to account for different insert shapes, where the effects of boundary scattering are not isotropic [10,11,12,13], and also to systems containing multiple particle sizes and orientations [14].

In composites formed of anisotropic materials (e.g., layered systems such as transition metal dichalcogenides) the directional dependence of the matrix thermal conductivity or dielectric constant must be accounted for [15,16,17]. If a Kapitza resistance exists in such a system, it is reasonable to assume that will also be directionally dependent and must be accounted for in a similar fashion (the principle of doing so within a Morika-Tanaki approach has been discussed in, for example, Ref. [18]). Within the context of the multiple scattering approach used in this study, an extension to Minnich and Chen’s [8] approach that accounts for anisotropic boundary resistance effects (emEMA) has been recently proposed in [19].

However, the originally formulated mEMA [8] and its recent extension [19] assume that the phonon scattering from insert boundaries is fully diffusive. In other words, these formulations treat the boundary as being perfectly rough. This is unlikely to be the case in real systems, since fully diffusive and specular scattering are idealised limiting cases. A simple phenomenological approach to the problem has been proposed [11], in which an empirically weighted average of the two limits is calculated. While this gives physically plausible results, there is no obvious connection between the value of the weighting parameter and the physical properties of the nanocomposite. In this work we suggest a scheme based upon Koh et al. [20] proposed correspondence between the relative proportions of different boundary scattering types and the average height of inhomogeneities on the boundary layer that is wave-vector dependent. In other words, our proposed extenstion improves on the approach adopted by Behrang et al. [11] for including both specular and diffuse scattering contributions to phonon boundary scattering rate as well as the thermal interface (Kapitza) scattering term.

Transition metal dichalcogenides (TMDCs) such as MoS${}_2$ and WS${}_2$ possess high electrical conductivities and low thermal conductivities. Because of this, they are thought to be good candidates for use in thermoelectric applications. Recent studies [21,22,23] demonstrate early attempts at fabricating TMDC-based nanostructures and composites. Through a judicious choice of insert and matrix material it will be possible to further reduce the thermal conductivity of these systems without strongly affecting the thermoelectric power factor, and so improve the thermoelectric figure of merit. In order to accurately predict the thermal conductivity reduction in a given nanocomposite, some means of accounting for the quality of the matrix-insert interface is required. Therefore in anticipation of future thermal conductivity measurements, we implement this scheme within the emEMA for a layered (hence anisotropic) bulk transition metal dichalcogenide nanocomposite consisting of spherical inserts of bulk 2H WS${}_2$ embedded within a matrix of bulk 2H MoS${}_2$, where semi-ab initio calculations [24] are used to obtain the required insert and matrix thermal conductivities.

2. Theoretical Method

An emEMA calculation for a system with anisotropic interface resistance is performed in two steps. Firstly we compute the effective thermal conductivity of the insert and the interface regions [19]:

$${\kappa }^{*}={\kappa }^i{\left[\mathbb{I}+\frac{3}{a}{\mathbb{N}}_s{\kappa }^i{\mathbb{R}}^K\right]}^{-1},$$

(1)

where ${\kappa }^i$ is the thermal conductivity of the insert in the absence of the boundary, ${\kappa }^{*}$ is the effective thermal conductivity of the insert and the boundary, a is the radius of the spherical insert, ${\mathbb{R}}^K$ is the thermal boundary or Kapitza resistance, and ${\mathbb{N}}_s$ is the surface depolarisation tensor (see Appendix A) accounting for the effects of anisotropic Kapitza resistance. The thermal conductivities, depolarisation tensor, Kapitza resistance and identity matrix $\mathbb{I}$ are $3\times 3$ matrices, and we assume that the system is oriented so that ${\kappa }^i$ and ${\mathbb{R}}^K$ are diagonal.

We then calculate the thermal conductivity of the effective medium ${\kappa }^E$ using the standard anisotropic EMA equation [15]:

$$\begin{array}{cc}\hfill {\kappa }^E& ={\kappa }^m+f({\kappa }^{*}-{\kappa }^m){\mathbb{J}}^{-1}{\kappa }^m,\hfill \\ \hfill \mathbb{J}& ={\kappa }^m+(1-f){\mathbb{N}}_m({\kappa }^{*}-{\kappa }^m).\hfill \end{array}$$

(2)

Here f is the volume fraction of inserts, ${\kappa }^m$ is the thermal conductivity of the matrix and ${\mathbb{N}}_m$ the matrix depolarisation tensor (see Appendix A) accounting for any anisotropy in ${\kappa }^m$. We assume the simplest possible case where ${\kappa }^m$ and ${\kappa }^{*}$ are both diagonal.

The required values of ${\kappa }^i$ and ${\kappa }^m$ may be computed via any reliable method for calculating the thermal conductivity. In this work we use a semi-ab initio method [24], which we summarise in Appendix B.

In order to include the effects of specularity, we define a momentum dependent parameter $s_{\mathit{q}}$ as the fraction of modes of momentum $\mathit{q}$ that undergo diffusive scattering from the insert interfaces. Generalising the simple expression of Koh et al. [20], we write:

$$s_{\mathit{q}}=1-l_{\mathit{q}},$$

(3)

where $l_{\mathit{q}}=exp(-{ϵ}^2/{\lambda }_{\mathit{q}}^2)$ with $ϵ$ being the average height of surface inhomogeneities and ${\lambda }_{\mathit{q}}$ being the wavelength corresponding to momentum $\mathit{q}$. We define $ϵ=\eta a_0$ where $a_0$ is the lattice spacing, and ${\lambda }_{\mathit{q}}=2\pi /Q$, where $Q=2\pi \left|\mathit{q}\right|/a_0$ and $\left|\mathit{q}\right|=\sqrt{q_x^2+q_y^2+q_z^2}$ is the norm of the lattice momentum vector in Cartesian co-ordinates. From these definitions it follows that:

$$s_{\mathit{q}}=(1-e^{-{\eta }^2{\left|\mathit{q}\right|}^2}).$$

(4)

The value of $\eta$ controls the extent of specular scattering. If $\eta \to 0$, then the surface becomes smooth and we obtain the specular limit $s_{\mathit{q}}=0$. If $\eta \to \infty$ then the surface becomes infinitely rough and we obtain the diffuse limit $s_{\mathit{q}}=1$. Physically meaningful values of $\eta$ will lie somewhere between these extremes, although in principle it is possible for a surface to be smooth enough or rough enough that all scattering is effectively specular or diffusive.

The effects of finite $\eta$ enter into our calculation on two levels: through the thermal boundary resistance at the EMA level, and through modifications of the effective boundary lengths of the matrix and insert at the mEMA level. We begin with the latter.

Scattering from the sample boundary itself is taken to be purely diffusive and unaffected by the specularity parameter. Minnich and Chen [8] include the effects of phonon scattering from insert interfaces by means of effective boundary lengths, and these must be modified so as to include the effects of specularity. In the case of the insert, any mode not scattered from the interface boundary due to specularity must be scattered from the sample boundary, and so the effective boundary length must be a weighted average of the sample length and the insert diameter. In the case of the matrix, only the interface density contribution to the effective boundary length is affected by the specularity. We may therefore write the the effective boundary length $L_B$ as follows:

$$L_B^{-1}=\left\{\begin{array}{cc}{s}_{\mathit{q}}{L}_{B,I}^{-1}+l_{\mathit{q}}{L}_{B,S}^{-1}\hfill & \mathrm{for}\mathrm{inserts}\hfill \\ L_{B,S}^{-1}+s_{\mathit{q}}{L}_{IS}^{-1}\hfill & \mathrm{for}\mathrm{the}\mathrm{matrix}\hfill \end{array}\right.,$$

(5)

where $L_{B,I}=2a$ is the insert boundary length, $L_{B,S}$ is the sample boundary length, $L_{IS}=4/\mathsf{\Phi }$ is the scattering length due to interface density $\mathsf{\Phi }=6f/L_{B,I}$ and $l_{\mathit{q}}=1-s_{\mathit{q}}$. In the case of the semi-ab initio thermal conductivity calculations performed in this study, the incorporation of these effective lengths into the overall scattering rate is discussed in Appendix B, specifically Equation (A10).

Next, we discuss the thermal boundary resistance. This is a challenging property to calculate with any accuracy, since the assumptions of the prevailing models are strictly valid only in the continuum limit and exclude the effects of inelastic scattering and other anharmonic processes (see Ref. [25] and citations within), and at best only give plausible bounds on the experimental value of the boundary resistance [4,25]. Nevertheless, we must qualitatively account for the effects of specularity on the thermal boundary resistance in some fashion, and so we follow in the spirit of previous theoretical work on mEMA-type calculations [11] by making use of an approximate model based on Chen’s expressions for the case of superlattice systems [26]. Our intent is to interpolate between the diffuse mismatch model (diffuse limit) and the acoustic mismatch model (specular limit) in a mode-dependent fashion reflecting the value of $s_{\mathit{q}}$, making use of the full band dispersions in our calculation. We use the equivalent equilibrium approach [26], and express $R_{TB}$ as the sum of weighted diffuse $R_{TB}^D$ and specular $R_{TB}^S$ contributions:

$$R_{TB}=R_{TB}^D+R_{TB}^S.$$

(6)

We begin with $R_{TB}^D$. Defining $\langle C_j{v}_j\rangle ={\sum }_{\mathit{q}s}{C}_{j,\mathit{q}s}{v}_{j,\mathit{q}s}$ and ${\langle C_j{v}_j\rangle }_W={\sum }_{\mathit{q}s}{C}_{j,\mathit{q}s}{v}_{j,\mathit{q}s}{s}_{\mathit{q}}$ where $j=i,m$ denotes the insert or matrix value respectively, $C_{j,\mathit{q}s}$ is the specific heat for the mode with momentum $\mathit{q}$ and band number s and $v_{j,\mathit{q}s}$ is the mode speed (which in the x or y direction is $\sqrt{v_x^2+v_y^2}$ and in the z direction is $|v_z|$), we may write:

$$R_{TB}^D\approx \frac{2({\langle C_i{v}_i\rangle }_W+{\langle C_m{v}_m\rangle }_W)}{\langle C_i{v}_i\rangle \langle C_m{v}_m\rangle }.$$

(7)

For $R_{TB}^S$, having defined the impedences $\langle Z_j\rangle ={\rho }_j{\sum }_{\mathit{q}s}{v}_{j,\mathit{q}s}{n}_{\mathit{q}s}/{\sum }_{\mathit{q}s}{n}_{\mathit{q}s}$ and ${\langle Z_j\rangle }_W={\rho }_j{\sum }_{\mathit{q}s}{v}_{j,\mathit{q}s}{l}_{\mathit{q}}{n}_{\mathit{q}s}/{\sum }_{\mathit{q}s}{n}_{\mathit{q}s}$ where ${\rho }_j$ is a density and $n_{\mathit{q}s}$ is the Bose-Einstein distribution, we can write the transmission coefficients:

$$\begin{array}{cc}\hfill {\mathcal{T}}_{mi}\left({\mu }_i\right)=& \frac{4\langle Z_i\rangle \langle Z_m\rangle }{{(\langle Z_i\rangle {\mu }_m+\langle Z_m\rangle {\mu }_m)}^2},\hfill \\ \hfill {\mathcal{T}}_{im}\left({\mu }_m\right)=& \frac{\langle C_m{v}_m^3\rangle }{\langle C_i{v}_i^3\rangle }{\mathcal{T}}_{mi}\left({\mu }_i\right),\hfill \\ \hfill {\mathcal{T}}_{mi}^W\left({\mu }_i\right)=& \frac{4\langle Z_i\rangle \langle Z_m\rangle }{{({\langle Z_i\rangle }_W{\mu }_i+{\langle Z_m\rangle }_W{\mu }_m)}^2},\hfill \end{array}$$

(8)

where ${\mu }_j=cos{\theta }_j$ where ${\theta }_j$ is the angle of incidence, and $\langle C_j{v}_j^3\rangle ={\sum }_{\mathit{q}s}{C}_{j,\mathit{q}s}{v}_{j,\mathit{q}s}^3$. Taking ${\mu }_{c,j}$ to be the value of ${\mu }_i$ corresponding to the critical angle above which all incident phonons are reflected, using the following integrals

$$\begin{array}{cc}\hfill I_{jk}=& {\int }_0^{{\mu }_{c,j}}{\mathcal{T}}_{jk}\left({\mu }_j\right){\mu }_{j}d{\mu }_j,\hfill \\ \hfill I_{mi}^W=& {\int }_0^{{\mu }_{c,m}}{\mathcal{T}}_{mi}^W\left({\mu }_m\right){\mu }_{m}d{\mu }_m,\hfill \end{array}$$

(9)

we may write:

$$R_{TB}^S\approx \frac{2}{\langle C_m{v}_m\rangle I_{mi}^W}(1-I_{mi}-I_{im}).$$

(10)

Note that $R_{TB}$ is not a linearly weighted average, since the weighting of the diffusive contribution enters its definition linearly whereas for the specular contribution it enters quadratically.

As an example calculation, we examine the effects of average inhomogeneity height $\eta$ on the effective thermal conductivity of bulk 2H WS${}_2$ embedded within a matrix of bulk 2H MoS${}_2$, where the relevent input thermal conductivities are obtained using a semi-ab initio [24] method described in Appendix B. We assume that the cross-plane and in-plane axes of both materials align, and in calculating the thermal boundary resistances we distinguish the in-plane values ${\mathbb{R}}_{xx}=R_{TB,x}={\mathbb{R}}_{yy}=R_{TB,y}$ from the cross-plane values ${\mathbb{R}}_{zz}=R_{TB,z}$ by setting $v=\sqrt{v_x^2+v_y^2}$ in the former case and $v=|v_z|$ in the latter. All other components of $\mathbb{R}$ are set to 0.

3. Results

Figure 1 displays the results varying $\eta$ in a system with $f=0.2$, $L_{B,S}=10$ $\mathsf{\mu }$, and $L_{B,I}=1$ $\mathsf{\mu }$m. Panel (a) displays $R_{TB}$ in the cross-plane (z) and in-plane (x) directions. In both directions it undergoes a rapid decrease as the temperature rises before saturating above 400 K. $R_{TB}$ at $\eta =5.0$ is close to that of the diffuse limit; in the range $0.1\ge \eta \ge 0.001$ we see little change in its value, suggesting that we have attained the specular limit. Any interesting physics regarding the boundary resistance must therefore occur within the range $5.0>\eta >0.1$. We find that $R_{TB,z}>R_{TB,x}$, which is consistent with ${\kappa }_{zz}$ being smaller than ${\kappa }_{xx}$ in the constituent materials. Note that the change in $R_{TB}$ as $\eta$ decreases is not monotonic: It decreases from the diffuse value, reaches a minimum in the vicinity of $\eta \approx 1.0$ to $\approx 0.5$ and increases towards the specular limiting value. This is a consequence of $s_{\mathit{q}}$ entering into our definition of the diffuse contribution as $s_{\mathit{q}}$ and into the specular contribution as ${(1-s_{\mathit{q}})}^2$ which would not be seen in a linear mixing approach such as that of Ref. [11].

Panel (b) of Figure 1 displays ${\kappa }^E$ in the cross-plane and in-plane directions. The cross-plane ${\kappa }_{zz}^E$ is, as we would expect from the constituents, smaller than ${\kappa }_{xx}^E$. The overall behaviour with decreasing $\eta$ in both directions is a monotonic increase from the value in the diffuse limit to saturation as the specular limit is approached for $\eta \le 0.1$. Note that while $\eta =5.0$ does not result in a hugely different value of $R_{TB}$ from that of the diffuse limit, for ${\kappa }^E$ the difference is far more apparent due to the effects of phonon scattering from the interface boundaries.

Results for ${\kappa }^m$ and ${\kappa }^i$ are displayed in panels (c) and (d) of Figure 1, respectively. Similarly to the ${\kappa }^E$ results we see that decreasing $\eta$ takes us monotonically from the diffuse limit (where scattering from interface boundaries is important) to the specular limit where there is no scattering from interface boundaries, only the sample edges. Interestingly, ${\kappa }^m$ appears to be slightly less sensitive to decreasing $\eta$ close to the specular limit, becoming saturated at $\eta =0.01$ rather than at $\eta =0.1$ as is seen for ${\kappa }^i$. This is a result of $L_{IS}$ (the contribution to the scattering length due to interface density) being larger in size than $L_{B,I}$; as $\eta$ tends towards zero the effective scattering length of the matrix will approach the value of the sample length much more quickly.

Figure 2 displays ${\kappa }_{xx}^E$ and ${\kappa }_{zz}^E$ for a number of insert sizes, volume fractions and $\eta$ values between 5.0 and 1.0, with sample size $L_{B,S}=10$ $\mathsf{\mu }$m. These are intended to capture the effects of physically plausible surface roughnesses, i.e., with deviations from the smooth surface of between 1 and 5 times the lattice constant in height. For all $\eta$ values, ${\kappa }^E$ in both directions increases with $L_{B,I}$ and decreases with f, as one would expect. The effect of decreasing $\eta$ is to reduce the effects of scattering from interfaces, and so we see enhancement of the peak in thermal conductivity below 100 K as the surface of the insert becomes more smooth. This enhancement is most noticable for 10 nm inserts (Figure 2a,b), which in the diffuse limit lack a pronounced peak. As the temperature is increased, the enhancement is decreased, and the thermal conductivities for all finite $\eta$ values considered tend towards the values seen in the diffuse limit. This is to be expected due to the dominance of anharmonic as opposed to length-based scattering at higher temperatures. The onset of anharmonic dominance generally occurs at lower temperatures for larger insert sizes, however, even at room temperature we can see from

Table 1. Effective thermal conductivities in units of W m${}^{-1}$ K${}^{-1}$ at different $\eta$, $L_{B,I}$, and volume fraction for $L_{B,S}=10$ $\mathsf{\mu }$m at 300 K.

	Volume Fraction		0.05				0.1				0.2
		$\mathit{\eta }$	Diff.	1.0	2.0	5.0	Diff.	1.0	2.0	5.0	Diff.	1.0	2.0	5.0
	${\mathit{L}}_{\mathit{B},\mathit{I}}$
${\mathit{\kappa }}_{\mathit{x}x}^E$	10 nm		47.9	73.2	59.3	51.7	28.9	52.5	39.8	31.7	15.5	33.7	23.4	16.9
	100 nm		110.2	121.9	114.7	112.5	87.8	104.4	94.1	90.8	61.9	80.6	69.2	61.2
	1000 nm		140.0	143.5	141.4	140.5	133.0	139.5	135.6	134.0	121.5	132.7	126.0	123.1
${\mathit{\kappa }}_{\mathit{z}z}^E$	10 nm		1.13	2.69	2.28	1.70	0.67	2.15	1.71	1.07	0.37	1.67	1.18	0.61
	100 nm		3.81	4.91	4.52	4.37	2.73	4.04	3.64	3.39	1.77	3.17	2.79	2.41
	1000 nm		6.41	6.76	6.59	6.56	5.67	6.22	5.97	5.93	4.67	5.44	5.13	5.05

that the effects of varying $\eta$ can still be quite significant for 10 nm and 100 nm inserts, particularly for $f=0.2$.

Why is the deviation from the diffuse result greater for smaller particles? Recall from Equation (5) that the effective boundary length $L_B$ for both the insert and the matrix is dependent on both the insert size $L_{B,I}$ and $s_{\mathit{q}}$. We may vary $s_{\mathit{q}}$ from 1 (diffuse limit) to 0 (specular limit) by decreasing $\eta$. A small $L_{B,I}$ means a smaller $L_B$ in the diffuse limit than would be the case with a $L_{B,I}$, and so for a small $L_B$ a wider range of values will be traversed as $s_{\mathit{q}}$ is varied than in the case where it is larger. The same decrease in $s_{\mathit{q}}$ will therefore cause a relatively greater change in $L_B$ for the case where $L_{B,I}$ is smaller.

For practical purposes, the above conclusion should be qualified slightly. We can estimate from the lattice constant of WS${}_2$ that for $L_{B,I}=10$ nm, there is a maximum meaningful value of $\eta \approx 16$, which is the point at which the average inhomogeneity height is equal to the radius of the insert. This means that it is entirely possible for $\eta >5.0$ and so for calculations considering only diffusive scattering to remain good models. However, if the quality of the interface-matrix boundary is sufficiently high, then the effects of specular scattering from the interfaces of small inserts must be included if the correct effective thermal conductivity is to be calculated.

In fact, comparing these results with those of Ref. [19] we find that the effects of specularity can significantly outweigh the effects of the corrections due to matrix and thermal boundary resistance anisotropy. Behrang et al’s calculation of the effects of specularity through simple mixing for spherical Si inserts in a Ge matrix [11] suggest that only a slight momentum independent specularity is needed to mimic Monte Carlo results. However, direct qualitative comparison between our results and theirs is difficult since they treat thermal boundary resistance differently from the standard approach, such that they predict smaller thermal conductivities than Minnich-Chen type [8] diffuse limit calculations for intermediate f and higher thermal conductivities for larger f, whereas we would expect a Minnich-Chen calculation (with appropriate modifications as in Ref. [19]) to provide the lower bound for our thermal conductivities.

4. Conclusions

We have used an extension of a modified effective medium approach (mEMA) to study thermal conduction in a transition metal dichalcogenide nanocomposite structure with rough interfaces by accounting for the anisotropic conductivities of the host, the insert and the interface regions. The role of specularity in phonon-interface scattering and in thermal boundary resistance is incorporated by implementing a simple phenomenological model relating momentum-dependent specular scattering and surface roughness. In general, it is found that the effect of specular scattering due to interface roughness is more pronounced for inserts smaller than 100 nm in the vicinity of 300 K. In particular, from comparison with calculations carried out in the diffuse limit, for spherical WS${}_2$ inserts in a matrix of MoS${}_2$, we predict that the effects of specular scattering due to surface roughness is more pronounced for smaller insert sizes, even at volume fractions of the order of 0.05. In other words, the surface roughness of the insert is a key parameter in controlling the thermal conductivity of nanocomposites with insert diameters of the order of tens of nanometers. For applications where the thermal conductivity should be as low as possible (e.g., thermoelectric applications) the surface of such inserts should be as rough as possible.