text:
# Surface scattering analysis of phonon transport in the quantum limit using an elastic model

## I Introduction

Landauer’s formulation of quantum transport showed that when elastic scattering dominates, the electrical conductance can be related to the transmission coefficient of the electron wavesL57 . In the ideal case of no scattering, this leads to a universal conductance that is quantized in units of $`e^2/h`$ at low temperatures, with an additional quantum of conductance added as each channel or mode of the conductance pathway opens up. The application of similar ideas to the phonon counterpart, namely thermal conductance, was recently derived by a number of authorsACR98 ; RK98 ; B99 , and is now recognizedBV00  to be related to earlier work on the entropy transport at low temperaturesP83 . Rego and Kirczenow have extended the concept of the universality of the thermal conductance to particles of arbitrary statistics (anyons)RK99 .

In the case of electrical resistance, the chemical potential or the number of conducting modes can be varied at very low temperatures, giving sharp jumps between various quantized values of the resistance. On the other hand, thermal transport by phonons necessarily requires nonzero temperatures to populate the modes of the conducting pathway, and the width of the Bose distribution function smears out the quantization of the conductance. Only at very low temperatures, where just the modes of the conducting pathway with zero frequency at long wavelengths contribute to the thermal conductance, the quantization of the ideal conductance becomes apparent in a universal thermal conductance $`N_0K_u`$ with $`K_u=(\pi ^2/3)k_B^2T/h`$ the universal conductance per mode, with $`k_B`$ Boltzmann’s constant and $`h`$ Planck’s constant, and $`N_0`$ is the number of modes with zero frequency at long wavelengths, which is four for a freely suspended elastic beam connecting the two thermal reservoirs. Note that this value of the low temperature conductance in the absence of scattering is independent of the dimensions and elastic properties of the thermal pathway.

A low temperature thermal conductance consistent with the predicted universal value was measured by Schwab et al.SHWR00  in experiments on a lithographically defined mesoscopic suspended beam (of dimensions about $`1\mu \mathrm{m}\times 200\mathrm{n}\mathrm{m}\times 60\mathrm{n}\mathrm{m}`$). Whilst their elegant experiment displays the universality of ballistic phonon transport, the experiment also showed a *decrease* in the thermal conductance below the universal value in the temperature range of $`0.08\mathrm{K}<T<0.4\mathrm{K}`$ that cannot be explained by the ballistic theory, since in this theory an increase in the thermal conductance is expected as the temperature is raised and more modes are excited. The decrease in thermal conductance is presumably associated with the scattering of the thermal phonons, and can be understood using the ideas of Landauer in terms of the scattering coefficient of the vibrational waves. This is the topic of the present paper.

In this paper we calculate the effect on the low temperature thermal conductance of the scattering of the thermal phonons by surface roughness, which is likely to be the major source of scattering in mesoscopic samples. The scattering of scalar waves, described by the simple wave equation, in waveguides with rough surfaces has been investigated by many workers, including ourselves, using both numerical and analytic methodsS97 ; KFFL99 ; SFMY99 ; B99 ; SC00 . However, for the low frequency modes of interest in the low temperature thermal conductance, the physical vibrational waves have quite different properties than the waves in the scalar model. For example the dispersion relations of the modes are different, with two of the four modes with zero long-wavelength frequency having a quadratic dispersion at small wave vectors, rather than the linear dependence given by the simple scalar theory. To understand the experimental results quantitatively, a more accurate treatment of the vibrational waves is needed. At low temperatures the wavelengths of the thermally excited modes are large compared with the atomic spacing, and so a treatment based on the equations of macroscopic elasticity theory is appropriate. BlencoweB95  has considered the scattering of elastic waves in a thin plate waveguide with rough surfaces, but prior to our work, the scattering of elastic waves confined in a beam-like wave guide with rough surfaces has not been considered.

Previously, we have investigated the effect of surface scattering on the low temperature thermal conductance using the scalar wave modelSC00 . In that paper we noted the apparent discrepancy between the results of scalar model with a simple assumption for the nature of the surface roughness and the data by Schwab et al. below a temperature of $`0.1\mathrm{K}`$: the data seemed to show a delay of the onset of scattering as the temperature increased that was not predicted by the model. However, since the scalar model does not properly account for the properties of the elastic waves, it was not clear whether this discrepancy is due to an inadequate modelling of the surface roughness, or the flaw in the description of the waves themselves. To resolve this matter, and obtain a more accurate account of the scattering of the waves by rough surfaces, we develop a theory of based on the full elasticity equations, and use this to calculate the thermal conductance at low temperatures. A short version of this work has been previously publishedSC01l .

In Section II, the scattering of elastic waves confined to a beam of rectangular cross-section with rough surfaces is calculated using the full three dimensional elasticity theory. We use a Green theorem approach, and calculate the scattering coefficient to quadratic order in the amplitude of the surface roughness. These results are quite general, but rather intractable for further progress, since the structure of the modes in an elastic beam cannot be determined in closed form. Thus in Section III we reduce the expressions to a thin plate limit to provide a closed form for the displacement fields, and to obtain analytical expressions for the scattering behavior. In Section IV the general behavior of the scattering and the effect on the thermal conductance is analyzed in detail, using a simple description of the surface roughness, to investigate the physical consequences of the novel features of the elastic waves. In Section V we use our theory to attempt to fit the data of Schwab et alSHWR00  using more realistic descriptions of the surface roughness. A number of the more difficult issues that arise in the elasticity theory are described in appendices.

Although our main interest is the scattering of thermally excited vibrational waves in mesoscopic systems at low temperatures, the formulation of the surface scattering is quite general, and can be applied to other situations, such as the scattering of mechanically excited modes in macroscopic samples for example.

## II General Formalism

### II.1 The model

The main focus of this paper is the effect of surface roughness on the low temperature thermal conductance of mesoscopic structures. The geometry we consider is a freely suspended elastic beam, which we call the bridge, connecting two thermal reservoirs. We will consider a beam of rectangular cross-section of dimensions width $`W`$ (in the $`y`$ direction) and depth $`d`$ (in the $`z`$ direction). Mesoscopic structures are often produced lithographically from epitaxially grown material. We choose a convention that the depth is the dimension in the growth direction, and the width in the lithographically defined transverse direction. We define the length of the rectangular beam of nominally uniform cross section as $`L`$. In practice the bridge may be joined to the reservoirs smoothly, by a portion of continuously growing width, to eliminate or reduce the scattering of the vibration modes off a sharp junction. We will suppose that the scattering by roughness is important only in some narrower portion of length $`L`$.

The thermal conductance is given by the expressionACR98 ; RK98 ; B99

$$K=\frac{\mathrm{}^2}{k_BT^2}\underset{m}{}\frac{1}{2\pi }_{\omega _m}^{\mathrm{}}𝒯_m(\omega )\frac{\omega ^2e^{\beta \mathrm{}\omega }}{(e^{\beta \mathrm{}\omega }1)^2}𝑑\omega ,$$
(1)
where $`\omega _m`$ is the cutoff frequency of the $`m`$ th mode, $`\beta =1/(k_BT)`$, and $`T`$ is the temperature. The integration is over the frequency $`\omega `$ of the modes $`m`$ that propagate in the structure. The transmission coefficient is unity for the ideal case. Any scattering reduces the thermal conductance, and scattering of the lowest modes can reduce the conductance below the universal value at low temperatures.

To actually perform the scattering calculation we imbed the rough beam of length $`L`$ in an infinite beam of the same cross section but with smooth surfaces outside of the region of length $`L`$, Fig. 1. Thus the mathematical calculation is the scattering of a wave incident from $`x=\mathrm{}`$ on a rough portion of the beam with surfaces at $`y=\pm W/2\pm f_1(x,z)`$ and at $`z=\pm d/2\pm f_2(x,y)`$, with the height functions $`f_{1,2}`$, defining the roughness, nonzero only in a finite region $`0<x<L`$. Forward scattering is evaluated from the intensity of waves as $`x+\mathrm{}`$, and backward scattering from the intensity of waves as $`x\mathrm{}`$.

To calculate the scattering amplitude, we take a Green function approach similar to our previous work on the scalar modelSC00 .

The displacement field $`𝐮`$ away from any sources satisfies the wave equation:

$$\rho _t^2u_i=_jT_{ij}$$
(2)
where $`\rho `$ is the mass density, and

$$T_{ij}=C_{ijkl}_ku_l$$
(3)
is the stress tensor field with $`C_{ijkl}`$ the elastic modulus tensor. The subscript $`i`$ runs over the $`3`$ Cartesian coordinates, we use the symbol $`_x`$ to denote the derivative $`/x`$ etc., and repeated indices are to be summed over. The displacement field satisfies stress free boundary conditions at the surfaces

$$T_{ij}n_j|_S=0$$
(4)
where $`S`$ denotes the surface boundaries and $`n_j`$ is normal to the surface. Assuming harmonic time dependence at frequency $`\omega `$, Eq. (2) becomes

$$\rho \omega ^2u_i+C_{ijkl}_j_ku_l=0.$$
(5)
We approximate the material of the system as an isotropic solid. Then the elastic modulus tensor is

$$C_{ijkl}=\lambda \delta _{ij}\delta _{kl}+\mu \left(\delta _{ik}\delta _{ji}+\delta _{il}\delta _{kj}\right)$$
(6)
where $`\lambda `$ and $`\mu `$ are Lamé constants ($`\mu `$ is also the shear modulus)

$$\lambda =E\sigma /(1+\sigma )(12\sigma ),\mu =E/2\left(1+\sigma \right)$$
(7)
with $`E`$ Young’s modulus and $`\sigma `$ the Poisson ratio.

Even in a rectangular beam geometry the displacement fields in the propagating modes yielded by these equations are complicated, and cannot be found analytically. The modes can be grouped into four classes according to their signature under the parity operations $`yy`$ and $`zz`$. Some modes show regions of anomalous dispersion where the group velocity $`d\omega /dk`$ is negative: these regions require a careful examination of the notions of “forward” and “backward” scattering for the waves. The lowest frequency mode of each class has a frequency that tends to zero at small wave number. These four modes are the only ones excited at low enough temperature, and are the ones contributing to the universal thermal conductance. The structure of these modes at small wave numbers is simple and can be calculated using familiar macroscopic arguments of elasticity theory: they are the compression, torsion, and (two orthogonal) bending modes.

We define a Green function $`G_{iq}(𝐱,𝐱^{};t,t^{})`$ to satisfy the wave equation with a source term $`\delta _{iq}\delta (𝐱𝐱^{})\delta (tt^{})`$, and $`\mathrm{\Gamma }_{ijq}`$ to be the corresponding stress

$$\mathrm{\Gamma }_{ijq}C_{ijkl}_kG_{lq}.$$
(8)
It is convenient to introduce the frequency space version of the Green function

$$G_{iq}(𝐱;𝐱^{};t,t^{})=\frac{d\omega }{2\pi }G_{iq}(𝐱;𝐱^{};\omega )e^{i\omega \left(tt^{}\right)},$$
(9)
with a similar expression defining $`\mathrm{\Gamma }_{ijq}(𝐱,𝐱^{};\omega )`$. Inserting $`G`$, $`\mathrm{\Gamma }`$, and the source term into Eq. (5) gives

$$\rho \omega ^2G_{iq}(𝐱,𝐱^{};\omega )+_j\mathrm{\Gamma }_{ijq}(𝐱,𝐱^{};\omega )=\delta _{iq}\delta \left(𝐱𝐱^{}\right)$$
(10)
where $`𝐱`$ is the observation coordinate and $`𝐱^{}`$ is the source coordinate.

Equations (5) and (10) lead to Green’s theorem expressing the displacement field at frequency $`\omega `$ in terms of a surface integral

$$u_q(𝐱)=_S^{}\left[n_j^{}T_{ij}\left(𝐱^{}\right)G_{iq}(𝐱^{},𝐱;\omega )n_j^{}u_i\left(𝐱^{}\right)\mathrm{\Gamma }_{ijq}(𝐱^{},𝐱;\omega )\right]𝑑S^{}.$$
(11)
We are free to choose any closed integration surface $`S^{}`$. One choice is to use the physical rough surface thereby eliminating the first term in Eq. (11) due to the boundary condition Eq. (4). However, the resulting integration over the rough surface is not easy. Instead, we integrate over the smoothed surfaces at $`y=\pm W/2`$ and $`z=\pm d/2`$ and impose the boundary conditions on the Green function to be stress free on these smoothed surfaces

$$\mathrm{\Gamma }_{ijq}n_j|_S=0,$$
(12)
together with cross sections at $`x^{}\pm \mathrm{}`$ to close the surface.

The total field $`𝐮`$ can be written as the sum of incident and scattered waves

$$𝐮=𝐮^{\mathrm{in}}+𝐮^{\mathrm{sc}}\text{.}$$
(13)
It can be shown (see Appendix A) that the integration over the sections at $`x^{}\pm \mathrm{}`$ on the right hand side of Eq. (11) just gives $`u_q^{\mathrm{in}}`$. In the integration over the smoothed surfaces at $`y=\pm W/2`$ and $`z=\pm d/2`$ the second term in the integrand vanishes due to Eq. (12). Thus we find the expression for the scattered field

$$u_q^{\mathrm{sc}}(𝐱)=_S\left[n_j^{}T_{ij}\left(𝐱^{}\right)G_{iq}(𝐱^{},𝐱;\omega )\right]𝑑S^{},$$
(14)
with the surface $`S`$ the smoothed surfaces $`y=\pm W/2`$ and $`z=\pm d/2`$. The stress field $`T_{ij}`$ on the smoothed surface is evaluated by expanding about its value on the *rough* surfaces, where Eq. (4) applies.

The rest of the section goes as follows: firstly, we find an explicit expression for the Green function with stress free boundary conditions; then we apply the boundary perturbation method to project the stress at the rough surfaces onto the smooth surfaces by expanding the stress free boundary terms around the smooth surfaces using the small roughness as the expansion parameter; and finally we evaluate the strength of the scattered waves to give the scattering coefficient.

### II.2 Green function

We evaluate $`G_{iq}(𝐱,𝐱^{};\omega )`$ as an expansion in the complete orthonormal set of normal modes $`𝐮^{(k,m)}\left(𝐱\right)`$ in the ideal geometry, which satisfy Eq. (5) and stress free boundaries at the smooth surfaces. Here $`k`$ is the wave number in the $`x`$ direction, and $`m`$ labels the branch of the dispersion curve. We define $`\omega _m(k)`$ as the frequency of the mode $`m`$ at wave number $`k`$ in the ideal geometry. The modes satisfy the completeness relation

$$\underset{m}{}\frac{dk}{2\pi }u_i^{(k,m)}(𝐱^{})^{}u_j^{(k,m)}(𝐱)=\delta _{ij}\delta (𝐱𝐱^{}).$$
(15)
Substituting this expression on the right hand side of Eq. (10) leads to the expression for the Green function

$$G_{iq}(𝐱^{},𝐱;\omega )=\underset{m}{}\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑k\frac{\varphi _i^{(k,m)}(y^{},z^{})^{}\varphi _q^{(k,m)}(y,z)}{\rho \left[\left(\omega +iϵ\right)^2\omega _m^2(k)\right]}e^{ik\left(xx^{}\right)},$$
(16)
where we write

$$u_i^{(k,m)}(𝐱)=\varphi _i^{(k,m)}(y,z)e^{ikx}$$
(17)
with$`\varphi _i^{(k,m)}`$ giving the transverse dependence of the displacement field. In Eq. (16) $`ϵ`$ is a positive infinitesimal number to incorporate causality, $`G_{iq}(𝐱,𝐱^{};t,t^{})=0`$ for $`t<t^{}`$.

Equation (16) can now be evaluated by contour integration. The integrand has poles labelled by an index $`n`$ near values $`k=k_n`$ on the real axis which are given by solutions to the dispersion relation $`\omega _m(k_n)=\omega `$ for all branches $`m`$. (We take an incident wave with $`\omega >0`$.) Note that for branches with regions of anomalous dispersion there may be more than one solution to this equation for some $`\omega `$, so that the index $`n`$ is not identical to the branch index $`m`$. The poles are shifted slightly off the real axis by the infinitesimal $`\epsilon `$ in Eq. (16), and are given by expanding about $`k_n`$

$$k=k_n+\frac{iϵ}{\omega _nv_g^{\left(n\right)}},$$

with $`v_g^{\left(n\right)}`$ the group velocity at the $`n`$th pole $`d\omega _m/dk|_{k=k_n}`$. Notice the poles are in the upper half plane for $`v_g^{(n)}>0`$, and in the lower half plane for $`v_g^{\left(n\right)}<0`$.

Now we can perform the $`k`$ integration by complex integration. Consider first the case, $`x>x^{}`$. The contour must be closed in the upper half plane so that the contribution from the semicircle at large $`|k|`$ vanishes. The contour integration then picks up contributions from the poles in the upper half plane, i.e. wave numbers with $`v_g^{\left(n\right)}>0`$. On the other hand, for $`x<x^{}`$ the contour must be closed in the lower half plane and it is poles at wave numbers with $`v_g^{\left(n\right)}<0`$ that give nonzero residue. Forward scattering or backscattering is thus seen to be determined by the sign of the group velocity $`v_g^{\left(n\right)}`$ rather than by the sign of $`k_n`$, as indeed would be expected physically.

Evaluating the residues gives the expression for the Green function:

$$G_{iq}(𝐱^{},𝐱;\omega )=i\underset{n}{}^{}\frac{u_i^{(n)}(𝐱^{})^{}u_q^{(n)}(𝐱)}{2\rho \omega _nv_g^{\left(n\right)}},$$
(18)
where $`u_i^{(n)}(𝐱)`$ is written for $`u_i^{(k,m)}(𝐱)`$ at the value of the wave number $`k=k_n`$ satisfying $`\omega _m(k_n)=\omega `$. The prime on the sum is used to denote the fact that the sum runs over intersections $`n`$ with $`v_g^{(n)}>0`$ for $`x>x^{}`$, and over $`n`$ with $`v_g^{(n)}<0`$ for $`x<x^{}`$.

The group velocity $`v_g^{\left(n\right)}`$ does not have an analytical expression for a rectangular beam, and is obtained by differentiating the dispersion curve which must be found numerically. Alternatively, to avoid numerical differentiation, we can rewrite $`v_g^{\left(n\right)}`$ in terms of the average power flow in mode $`n`$. Since $`u_i^{\left(n\right)}`$ is normalized, the power $`P_n`$ in mode $`n`$ can be written as

$$P_n=\frac{1}{2}\mathrm{Re}\left(i\omega T_{ix}^{\left(n\right)}u_i^{\left(n\right)}\right)𝑑y𝑑z=\frac{1}{2}\rho \omega ^2v_g^{\left(n\right)},$$
(19)
the first expression of the equality expressing the energy flux in terms of the rate of work done across a section, and the second in terms of the group velocity and the average energy density evaluated as twice the average kinetic energy. Then $`v_g^{\left(n\right)}`$ can be evaluated in terms of $`P_n`$ as

$$v_g^{\left(n\right)}=2P_n/\rho \omega ^2$$
(20)
and $`P_n`$ has an expression directly in terms of displacement field given by the first equality in Eq. (19)

$$P_n=\frac{1}{2}\mathrm{Re}\left(i\omega T_{ix}^{\left(n\right)}u_i^{\left(n\right)}\right)𝑑y𝑑z.$$
(21)
This expression for $`v_g^{\left(n\right)}`$ can also be derived directly from the equations of motion A .

### II.3 Boundary perturbation

In this section we show the boundary perturbation technique for the rough surfaces on the sides (i.e. the $`xz`$ boundary planes). We work out the scattering coefficient explicitly for the surface near $`y=W/2`$. The surface near $`y=W/2`$ will give a similar contribution and, assuming uncorrelated roughness on the two surfaces, is accounted for by multiplying the single-surface scattering rate by $`2`$ at the end of the calculation. The results for the top and bottom surfaces can be obtained by interchanging $`y`$ and $`z`$ whenever they occur in the indices in the displacement fields and stress tensors in the calculation below.

In order to calculate the stress on the smooth surface appearing in Eq. (14), we expand the stress $`T_{ij}`$ in a Taylor series about the flat surface, and impose stress free boundary conditions at the rough surface which is the small distance $`f_1`$ away. We also assume $`f_1`$ is differentiable.

The unit vector $`\widehat{n}`$ normal to the rough boundaries to first order in $`f_1`$ is

$$\widehat{n}\widehat{y}_xf_1(x,z)\widehat{x}_zf_1(x,z)\widehat{z}.$$
(22)
Then the stress free surface boundary conditions Eq. (4) can be written

$$\left[T_{iy}_xf_1(x,z)T_{ix}_zf_1(x,z)T_{iz}\right]_{y=\frac{W}{2}+f}=0$$
(23)
Now expanding Eq. (23) in the neighborhood of $`y=W/2`$ and taking only the lowest order in $`f_1`$ and $`f_1^{}`$, we obtain

$`T_{xy}|_{y=\frac{W}{2}}`$ $`\left(_xf_1(x,z)T_{xx}+_zf_1(x,z)T_{xz}f_1(x,z)_yT_{xy}\right)|_{y=\frac{W}{2}}`$ (24)
$`T_{zy}|_{y=\frac{W}{2}}`$ $`\left(_xf_1(x,z)T_{zx}+_zf_1(x,z)T_{zz}f_1(x,z)_yT_{zy}\right)|_{y=\frac{W}{2}}`$ (25)
$`T_{yy}|_{y=\frac{W}{2}}`$ $`f_1(x,z)_yT_{yy}|_{y=\frac{W}{2}}`$ (26)
where the first two expressions for $`T_{xy}`$ and $`T_{zy}`$ have been used to simplify $`T_{yy}`$. Since the terms on the right hand side of Eqs. (24-26) are explicitly first order in the small parameter $`f_1`$, the stress field $`T_{ij}`$ on the right hand side can be evaluated at zeroth order, i.e. for ideal smooth surfaces. These results are used in Eq. (14).

### II.4 Scattering coefficient

We now evaluate the expression for the scattered field given by an integration over the beam surfaces Eq. (14). To calculate the scattering coefficient, we consider an incident wave of unit amplitude in a single mode $`m`$. Again in this section we will outline the calculation for the scattering by the single surface at $`y=W/2`$, and will include the effects of the other surfaces at the end. We therefore have

$$u_q^{\text{sc}}\left(𝐱\right)=\left[T_{iy}\left(𝐱^{}\right)G_{iq}(𝐱^{},𝐱;\omega )\right]_{y^{}=\frac{W}{2}}𝑑x^{}𝑑z^{}.$$
(27)
We can now evaluate the forward and backscattering amplitudes by using Eq. (18) for the Green function in Eq. (27), and evaluating the scattered wave at large positive and negative $`x`$

$$u_q^{\mathrm{sc}}\left(x\mathrm{},y,z\right)_{\mathrm{}}^x𝑑x^{}_{\frac{d}{2}}^{\frac{d}{2}}𝑑z^{}\underset{n,v_g^{(n)}>0}{}\frac{i}{2\rho \omega v_g^{\left(n\right)}}\left[T_{iy}\left(𝐱^{}\right)u_i^{\left(n\right)}\left(𝐱^{}\right)^{}\right]_{y^{}=\frac{W}{2}}u_q^{\left(n\right)}\left(𝐱\right),$$
(28)
$$u_q^{\mathrm{sc}}\left(x\mathrm{},y,z\right)_x^{\mathrm{}}𝑑x^{}_{\frac{d}{2}}^{\frac{d}{2}}𝑑z^{}\underset{n,v_g^{(n)}<0}{}\frac{i}{2\rho \omega v_g^{\left(n\right)}}\left[T_{iy}\left(𝐱^{}\right)u_i^{\left(n\right)}\left(𝐱^{}\right)^{}\right]_{y^{}=\frac{W}{2}}u_q^{\left(n\right)}\left(𝐱\right).$$
(29)
The stress tensor $`T_{ij}`$ corresponding to the full displacement field of the wave is evaluated from Eqs. (24-26). Since these expressions explicitly include the small roughness amplitude $`f_1`$ on the right hand side, to calculate the scattering at lowest order in the roughness amplitude it is sufficient to replace all $`T_{ij}`$ on the right hand side by the value $`T_{ij}^{(m)}`$ in the incident mode $`m`$. From Eqs. (28), (29) we see that $`𝐮^{\text{sc}}(𝐱)`$ is expressed as a sum over modes $`𝐮^{\left(n\right)}(𝐱)`$, and the coefficient of each mode is then the scattering amplitude $`t_{n,m}`$ from incident mode $`m`$ into mode $`n`$, so that

$$t_{n,m}=_{\mathrm{}}^{\mathrm{}}𝑑x_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\frac{i}{2\rho \omega v_g^{\left(n\right)}}\left[T_{iy}^{(m)}\left(𝐱\right)u_i^{\left(n\right)}\left(𝐱\right)^{}\right]_{y=W/2},$$
(30)
where we can now extend the integration limit to $`\pm \mathrm{}`$ since $`f_1`$, and so the integrand, is zero outside the domain of roughness $`0<x<L`$. Again mode indices $`n`$ for which $`v_g^{(n)}>0`$ represent the forward-scattered waves and those with $`v_g^{(n)}<0`$ the backward-scattered waves.

Now use the expression for the stress tensor on the smooth surfaces obtained in the previous section Eqs. (24-26) and integrate the resulting expressions by parts with respect to $`x`$ or $`z`$ to rewrite the terms in $`_xf_1`$ and $`_zf_1`$ as integrations over $`f_1`$. After these manipulations we find $`t_{n,m}`$ can be written

$$t_{n,m}=\frac{i}{2\rho \omega v_g^{\left(n\right)}}_{\mathrm{}}^{\mathrm{}}𝑑x_{\frac{d}{2}}^{\frac{d}{2}}𝑑zf_1(x,z)\mathrm{\Gamma }^{(m,n)}(x,z)$$
(31)
where

$`\mathrm{\Gamma }^{(m,n)}(x,z)`$ $`=[(_xT_{xx}^{\left(m\right)}+_yT_{xy}^{\left(m\right)}+_zT_{xz}^{\left(m\right)})u_x^{\left(n\right)}+(_xT_{zx}^{\left(m\right)}+_yT_{zy}^{\left(m\right)}+_zT_{zz}^{\left(m\right)})u_z^{\left(n\right)}`$
$`+_yT_{yy}^{\left(m\right)}u_y^{\left(n\right)}+T_{xx}^{\left(m\right)}_xu_x^{\left(n\right)}+T_{zz}^{\left(m\right)}_zu_z^{\left(n\right)}+T_{zx}^{\left(m\right)}(_xu_z^{\left(n\right)}+_zu_x^{\left(n\right)})]_{y=W/2}.`$ (32)
Applying the equations of motion Eq. (5) and remembering$`T_{iy}^{(m)}|_{y=W/2}=0`$ for all $`i`$ and for all $`x,z`$ leads to the somewhat simpler expression

$`\mathrm{\Gamma }^{(m,n)}(x,z)`$ $`=[\rho \omega ^2(u_x^{(m)}u_x^{\left(n\right)}+u_y^{(m)}u_y^{\left(n\right)}+u_z^{(m)}u_z^{\left(n\right)})`$
$`+T_{xx}^{\left(m\right)}_xu_x^{\left(n\right)}+T_{zz}^{\left(m\right)}_zu_z^{\left(n\right)}+T_{xz}^{\left(m\right)}(_zu_x^{\left(n\right)}+_xu_z^{\left(n\right)})]_{y=W/2}.`$ (33)
Notice that the scattering separates into a kinetic term (the first line) and a stress term (the second line).

The above form for $`\mathrm{\Gamma }^{(m,n)}`$ is still neither instructive nor practical for numerical evaluation. It can be further simplified using the expressions Eqs. (3) and (7) for the stress tensor in terms of the displacements. First we use the boundary condition $`T_{yy}^{(m)}=0`$ for the y-stress to give at $`y=W/2`$

$$_yu_y^{(m)}=\frac{\sigma }{\left(1\sigma \right)}(_xu_x^{(m)}+_zu_z^{\left(m\right)}).$$
(34)
This can be used to simplify the expressions for the other components of the stress tensors at $`y=W/2`$

$`T_{xx}^{(m)}`$ $`={\displaystyle \frac{E}{\left(1\sigma ^2\right)}}(_xu_x^{(m)}+\sigma _zu_z^{\left(m\right)}),`$ (35)
$`T_{zz}^{(m)}`$ $`={\displaystyle \frac{E}{\left(1\sigma ^2\right)}}(\sigma _xu_x^{(m)}+_zu_z^{\left(m\right)}),`$ (36)
$`T_{xz}^{(m)}`$ $`={\displaystyle \frac{E}{2\left(1+\sigma \right)}}(_xu_z^{(m)}+_zu_x^{(m)}).`$ (37)
Inverting these gives at $`y=W/2`$

$`_xu_x^{\left(m\right)}`$ $`={\displaystyle \frac{1}{E}}(T_{xx}^{\left(m\right)}\sigma T_{zz}^{\left(m\right)}),`$ (38)
$`_zu_z^{\left(m\right)}`$ $`={\displaystyle \frac{1}{E}}(T_{zz}^{\left(m\right)}\sigma T_{xx}^{\left(m\right)}),`$ (39)
$`_xu_z^{(m)}+_zu_x^{(m)}`$ $`={\displaystyle \frac{2(1+\sigma )}{E}}T_{xz}^{\left(m\right)}.`$ (40)
We emphasize that Eqs. (34-40) are only true on the stress free boundaries, and are not generally true in the bulk of the material.

Using these results we get

$$t_{n,m}=\frac{i}{2\rho \omega v_g^{\left(n\right)}}_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\stackrel{~}{f}_1(k_mk_n,z)\overline{\mathrm{\Gamma }}^{(m,n)}(z)$$
(41)
with

$`\overline{\mathrm{\Gamma }}^{(m,n)}`$ $`=\{\rho \omega ^2(\varphi _x^{(m)}\varphi _x^{\left(n\right)}+\varphi _y^{(m)}\varphi _y^{\left(n\right)}+\varphi _z^{(m)}\varphi _z^{\left(n\right)})`$ (42)
$`{\displaystyle \frac{1}{E}}[(\overline{T}_{xx}^{\left(m\right)}\overline{T}_{zz}^{\left(n\right)}+\overline{T}_{zz}^{\left(m\right)}\overline{T}_{zz}^{\left(n\right)})(\sigma \overline{T}_{zz}^{\left(m\right)}\overline{T}_{xx}^{\left(n\right)}+\overline{T}_{xx}^{\left(m\right)}\overline{T}_{zz}^{\left(n\right)})]{\displaystyle \frac{1}{\mu }}\overline{T}_{xz}^{\left(m\right)}\overline{T}_{zx}^{\left(n\right)}\}_{y=W/2}`$
where we have introduced the explicit $`x`$ dependence of $`u_i^{(n)}(𝐱)`$ as in Eq. (17) and the stress tensor

$$T_{ij}^{(n)}(𝐱)=\overline{T}_{ij}(y,z)e^{ik_nx},$$
(43)
so that the $`x^{}`$ integration is just the Fourier transform $`\stackrel{~}{f}`$ of the roughness function, and $`\overline{\mathrm{\Gamma }}`$ is a function of the $`z`$ coordinate only.

Alternatively, using Eqs. (35-37) we can derive an expression explicitly in the displacement fields, which is useful for numerical evaluation,

$`\overline{\mathrm{\Gamma }}^{(m,n)}`$ $`=\{\rho \omega ^2(\varphi _x^{(m)}\varphi _x^{\left(n\right)}+\varphi _y^{(m)}\varphi _y^{\left(n\right)}+\varphi _z^{(m)}\varphi _z^{\left(n\right)})`$
$`{\displaystyle \frac{2\mu }{\left(1\sigma \right)}}\left[(k_mk_n\varphi _x^{(m)}u_x^{\left(n\right)}+_z\varphi _z^{\left(m\right)}_z\varphi _z^{\left(n\right)})+\sigma (ik_m\varphi _x^{(m)}_z\varphi _z^{\left(n\right)}ik_n_z\varphi _z^{\left(m\right)}\varphi _x^{\left(n\right)})\right]`$
$`\mu (ik_m\varphi _z^{(m)}_z\varphi _x^{\left(n\right)}+k_mk_n\varphi _z^{(m)}\varphi _z^{\left(n\right)}+_z\varphi _x^{(m)}_z\varphi _x^{\left(n\right)}ik_n_z\varphi _x^{(m)}\varphi _z^{\left(n\right)})\}_{y=W/2}.`$ (44)
The scattering rate is given by multiplying $`\left|t_{n,m}\right|^2`$ by the ratio of the group velocities in the scattered and incident waves<sup>1</sup><sup>1</sup>1Note that the scattering amplitude normalized by the energy flux $`\overline{t}_{n,m}=\sqrt{\left|v_g^{\left(n\right)}\right|/\left|v_g^{\left(m\right)}\right|}t_{n,m}`$ can be seen from Eq. (41) to explicitly satisfy the reciprocity relation $`\overline{t}_{n,m}=\overline{t}_{m,n}^{}`$. For systems in which energy is conserved, as in our case, reciprocity is equivalent to time-reversal invarianceP98 .. We also treat the roughness of the surface statistically, and take an ensemble average (denoted by angular brackets) to give the final expression for the scattering rate $`\gamma _{n,m}`$ from mode $`m`$ to mode $`n`$ by the per unit length of single rough surface at $`y=W/2`$ given by

$$\gamma _{n,m}L=\frac{v_g^{\left(n\right)}}{v_g^{\left(m\right)}}\left|t_{n,m}\right|^2=\frac{1}{4\rho ^2\omega ^2v_g^{\left(m\right)}v_g^{\left(n\right)}}\left|_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\stackrel{~}{f}_1(k_mk_n,z)\overline{\mathrm{\Gamma }}^{(m,n)}(z)\right|^2.$$
(45)
We are interested in the reduction of the phonon heat transport due to rough surfaces. Only the backscattered waves (those with $`v_g^{\left(n\right)}<0`$) reduce the amount of heat transmitted. Thus we define $`\gamma _m`$, the thermal attenuation coefficient of mode $`m`$ per unit length, to be the sum of the scattering rates from the incident mode $`m`$ to *all* possible backscattered modes, per unit length of rough surface. This can be written for scattering off the single rough surface considered so far

$$\gamma _mL=\underset{\begin{array}{c}n\\ v_g^{(n)}<0\end{array}}{}\gamma _{n,m}L=\underset{\begin{array}{c}n\\ v_g^{(n)}<0\end{array}}{}\frac{1}{4\rho ^2\omega ^2v_g^{\left(m\right)}v_g^{\left(n\right)}}\left|_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\stackrel{~}{f}_1(k_mk_n,z)\overline{\mathrm{\Gamma }}^{(m,n)}(z)\right|^2.$$
(46)
To include the second rough side surface, assuming uncorrelated roughness, we simply have to multiply the expression for $`\gamma _m`$ by a factor of $`2`$. The expression for scattering off the top and bottom surfaces, if these are rough too, can be derived in a similar manner, and the result may be obtained by exchanging $`y`$ and $`z`$ in Eq. (46). The total scattering rate is the sum of the scattering off all the surfaces.

We have assumed that the amplitude of the surface roughness is small, allowing us to use perturbation theory to derive the above expressions. In this case, the attenuation coefficient gives an exponential decay of the energy flux, so that the energy flux transmission coefficient is

$$𝒯_m=\mathrm{exp}\left[\gamma _mL\right].$$
(47)
## III Thin Plate Limit

Although the expression in the previous section is general and applicable to any rectangular waveguide with rough surfaces, there are no closed-form expressions for the displacement fields in general, and so a direct evaluation of the scattering would have to be done completely numerically. Here, we instead use the *thin plate approximation* $`dW`$ LF7 ; CL00 , which yields closed form expressions for the displacement fields of the modes (in terms of a dispersion curves $`\omega _m(k)`$ given by numerical solution of a simple transcendental equation). The thin plate limit captures the important properties of the elastic modes, for example the quadratic dispersion of the bending modes at small wave numbers, and regions of anomalous dispersion, as well as providing analytical expressions enabling us to do further analysis of the scattering. The thin plate theory is applicable where the thickness of the sample is much less than the width and the wavelengths are much greater than the thickness, which is the case for many mesoscopic systems at low temperatures.

The use of the thin plate limit for mesoscopic structures was proposed in reference CL00 , where the calculation of the structure of the modes is described in more detail. It is found that the modes can be separated into two classes: *in-plane modes*, where the polarization of the displacement is largely in the $`xy`$ plane (together with small strains in the $`z`$-direction given by the Poisson effect) and the displacement field is completely specified by giving the *vertically averaged horizontal displacement components* $`\overline{u}_x(x,y)`$ and $`\overline{u}_y(x,y)`$; and *flexural modes*, where the displacement is primarily in the $`z`$ direction and is specified by a vertical displacement field $`\overline{u}_z(x,y)`$. Within each class we can further distinguish the modes by their parity under $`yy`$. For the in-plane modes we define the mode as even if $`\overline{u}_x(x,y)=\overline{u}_x(x,y)`$ and odd if $`\overline{u}_x(x,y)=\overline{u}_x(x,y)`$. Similarly, the even flexural modes have $`\overline{u}_z(x,y)=\overline{u}_z(x,y)`$ and the odd modes have $`\overline{u}_z(x,y)=\overline{u}_z(x,y)`$. As in the general case, there are four branches of the dispersion curves that tend to zero frequency as the wave number goes to zero, corresponding to one mode from each of these classes. The low frequency, even in plane mode corresponds to the compression mode, and the odd mode to a bending mode. The low frequency even, flexural mode corresponds to the second bending mode, and the low frequency odd flexural mode is the torsion mode.

Explicit expressions for the displacement fields can be obtained using the method described in reference CL00 . For the in-plane modes we find, up to a normalization factor $`A_1`$ that is common to both even and odd parity waves, the even modes

$$\overline{u}_x(x,y)=ikA_1\left[\frac{k^2\chi _1^2}{2k^2}\mathrm{cos}\left(\frac{\chi _2W}{2}\right)\mathrm{cos}\left(\chi _1y\right)\mathrm{cos}\left(\chi _2y\right)\mathrm{cos}\left(\frac{\chi _1W}{2}\right)\right]e^{ikx}$$
(48)
$$\overline{u}_y(x,y)=A_1\left[\frac{k^2\chi _1^2}{2\chi _1}\mathrm{cos}\left(\frac{\chi _2W}{2}\right)\mathrm{sin}\left(\chi _1y\right)+\chi _2\mathrm{cos}\left(\frac{\chi _1W}{2}\right)\mathrm{sin}\left(\chi _2y\right)\right]e^{ikx},$$
(49)
and the odd modes

$$\overline{u}_x(x,y)=ikA_1\left[\frac{k^2\chi _1^2}{2k^2}\mathrm{sin}\left(\chi _1y\right)\mathrm{sin}\left(\frac{\chi _2W}{2}\right)\mathrm{sin}\left(\frac{\chi _1W}{2}\right)\mathrm{sin}\left(\chi _2y\right)\right]e^{ikx}$$
(50)
$$\overline{u}_y(x,y)=A_1\left[\frac{k^2\chi _1^2}{2\chi _1}\mathrm{cos}\left(\chi _1y\right)\mathrm{sin}\left(\frac{\chi _2W}{2}\right)+\chi _2\mathrm{sin}\left(\frac{\chi _1W}{2}\right)\mathrm{cos}\left(\chi _2y\right)\right]e^{ikx},$$
(51)
where $`\chi _1=\left(\omega ^2/c_t^2k^2\right)^{1/2}`$ and $`\chi _2=\left(\omega ^2/c_l^2k^2\right)^{1/2}`$, with $`c_t`$ the transverse sound velocity and $`c_l`$ the longitudinal velocity *in a large thin plate*

$$c_t=\sqrt{\frac{E}{2\rho \left(1+\sigma \right)}},c_l=\sqrt{\frac{E}{\rho \left(1\sigma ^2\right)}},$$
(52)
and $`\omega `$ and $`k`$ related by the dispersion curve which must be found numerically. In the thin plate limit it is sufficient to take for the in-plane modes

$`u_x(x,y,z)`$ $`\overline{u}_x(x,y),`$ (53)
$`u_y(x,y,z)`$ $`\overline{u}_y(x,y),`$ (54)
$`u_z(x,y,z)`$ $`0.`$ (55)
Similarly, the vertical displacement field for the even flexural modes is

$$\overline{u}_z(x,y)=A_2\left[\mathrm{cosh}\left(\frac{\chi _{}W}{2}\right)\mathrm{cosh}\left(\chi _+y\right)\frac{k^2\sigma \chi _+^2}{k^2\sigma \chi _{}^2}\mathrm{cosh}\left(\frac{\chi _+W}{2}\right)\mathrm{cosh}\left(\chi _{}y\right)\right]e^{ikx},$$
(56)
and for the odd flexural modes

$$\overline{u}_z(x,y)=A_2\left[\mathrm{sinh}\left(\frac{\chi _{}W}{2}\right)\mathrm{sinh}\left(\chi _+y\right)\frac{k^2\sigma \chi _+^2}{k^2\sigma \chi _{}^2}\mathrm{sinh}\left(\frac{\chi _+W}{2}\right)\mathrm{sinh}\left(\chi _{}y\right)\right]e^{ikx},$$
(57)
where $`\chi _+=\left(k^2+\sqrt{\rho d/D}\omega \right)^2`$ and $`\chi _{}=\left(k^2\sqrt{\rho d/D}\omega \right)^2`$, with $`D=Ed^3/12\left(1\sigma ^2\right)`$ the flexural rigidify, and again $`\omega `$ and $`k`$ are related by the appropriate dispersion curve. In the classical thin plate theory, the displacement fields are given in terms of $`\overline{u}_z`$ by the expressions

$`u_x(x,y,z)`$ $`z_x\overline{u}_z(x,y),`$ (58)
$`u_y(x,y,z)`$ $`z_y\overline{u}_z(x,y),`$ (59)
$`u_z(x,y,z)`$ $`\overline{u}_z(x,y).`$ (60)
This approximation is adequate for evaluating the surface stress integrals in Eq. (46) but turns out not to be sufficiently accurate to evaluate the energy flux expression for the group velocity Eq. (21). We discuss this case in section III.3.2 below.

### III.1 Ideal thermal conductance

Since our quantitative calculation of the scattering coefficient relies on the analytic expressions for the elastic modes available only in the thin plate limit, it is essential to estimate the temperature range where the thin plate limit is applicable for a given experimental structure. On the other hand, as the wavelength becomes much smaller than the dimensions of the structure, we should to be able to treat the waves in terms of separate longitudinal and transverse waves in the bulk of the material, without worrying too much about the complicated standing wave transverse mode structure important for the long wavelength modes. In this regime, which we refer to as the bulk mode limit, the counting of the modes is insensitive to the details of the boundary conditions, and is the same as for a scalar wave approximation. The ideal thermal conductance depends only on cut-off frequency of the modes (see Eq. (1)), and we can assess the applicability of these simple limiting approximations by comparing the mode cutoff frequencies with results from a numerical calculation of the full elasticity theory. For the full elastic theory, we use the “xyz algorithm”NAW97 .

For the bulk mode calculation, there are three polarizations (one longitudinal and two transverse) with propagation velocities $`c_{3l}`$ and $`c_t`$, respectively with $`c_t`$ as before and

$$c_{3l}=\sqrt{\frac{E\left(1\sigma \right)}{\rho \left(1+\sigma \right)\left(12\sigma \right)}}.$$
(61)
The precise details of the boundary conditions are unimportant in the mode counting for large mode numbers. If we assume standing waves in the transverse direction corresponding to zero normal derivative boundary conditions on the wave functions, the cutoff frequencies are

$$\omega _{t,mn}=c_t\sqrt{\left(\frac{m\pi }{W}\right)^2+\left(\frac{n\pi }{d}\right)^2},\text{ (two fold degenerate)}$$
(62)
for the transverse waves, and

$$\omega _{l,mn}=c_{3l}\sqrt{\left(\frac{m\pi }{W}\right)^2+\left(\frac{n\pi }{d}\right)^2},\text{ (non degenerate)}$$
(63)
for the longitudinal waves, with $`m,n=0,1,2\mathrm{}`$. For large $`m,n`$ we can use the continuous form for the frequency $`\omega _N`$ of the $`N`$th mode

$$N=\frac{dW}{4\pi }\omega _N^2\left(\frac{2}{c_t^2}+\frac{1}{c_{3l}^2}\right).$$
(64)
Figure 2 shows the cutoff frequencies as a function of mode number for the thickness to width ratio $`d/W=0.38`$. The thin plate theory gives a good approximation at lower frequencies. The accuracy of thin plate theory becomes better as $`d/W`$ gets smaller. For example, in the case of $`d/W=0.1`$ (not shown), the error in the cutoff frequencies of the first $`13`$ modes is less than $`3\%`$, whilst the error is as large as $`5\%`$ for the first $`7`$ modes for the case $`d/W=0.38`$ shown in the figure. In terms of the ideal (no-scattering) thermal conductance (Eq. (1) with the transmission coefficient set to unity), we find that for $`d/W=0.38`$ the error in the thermal conductance is less than $`4\%`$ up to $`T0.4`$ K. Thus the thin plate limit is adequate to examine the scattering effects in this temperature range. At large frequencies $`\omega W/c_t>30`$ the elasticity theory results approach closely the continuum bulk mode calculations, Eq. (64). The thin plate approximation clearly fails in this limit, since it predicts $`N\omega `$ corresponding to a $`2d`$ structure.

### III.2 Attenuation coefficient in the thin plate limit

The thin plate approximation is implemented by noticing that the stress free boundary conditions imply that the stress components $`T_{iz}`$ are zero on the top and bottom surfaces. For small thickness this implies that the components $`T_{iz}`$ for any $`i`$ are small everywhere. In most situations these components can be approximated as zeroLF7 . This simplifies many of the terms appearing in Eq. (41). Also, at low temperatures, only modes with no strong dependence on the $`z`$ coordinate will be excited, so that the mode sum extends over modes with increasing numbers of nodes in the $`y`$ direction only.

In this section we calculate the scattering of the elastic waves by surface roughness for a thin plate. We assume that the roughness is confined to the sides, since in the experiments theses are prepared lithographically, whereas the top and bottom surfaces are produced by the epitaxial growth process.

For simplicity, we assume the roughness function $`f_1`$ has no $`z`$ dependence. This is probably a reasonable description of the roughness produced by a typical lithographic process of anisotropic chemical etch<sup>2</sup><sup>2</sup>2Without this assumption, we would find slightly different $`z`$ averages of the roughness function $`\stackrel{~}{f}_1`$ involved for the scattering of the in-plane modes and of the flexural modes—in fact a direct average for the inplane modes and an average weighted by $`z^2`$ for the flexural modes. In addition there would now be scattering from in-plane to flexural modes, and *vice versa*.. Then the Fourier transformed roughness function $`\stackrel{~}{f}_1\left(k_mk_n\right)`$ can be pulled outside of the $`z`$ integral in Eq. (46) and the statistical average over the roughness can be performed to give

$$\left|\stackrel{~}{f}_1\left(k\right)\right|^2=\stackrel{~}{g}\left(k\right)L,$$
(65)
where $`\stackrel{~}{g}\left(k\right)`$ is the Fourier transform of the roughness correlation function

$$\stackrel{~}{g}\left(k\right)=𝑑xe^{ikx}f_1\left(x\right)f_1\left(0\right).$$

Equation (46) leads to the back-scattering rate from mode $`m`$ to mode $`n`$

$$\gamma _{n,m}=\frac{\stackrel{~}{g}\left(k_mk_n\right)}{2\rho ^2\omega ^2v_g^{\left(m\right)}v_g^{\left(n\right)}}\left|_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\overline{\mathrm{\Gamma }}^{(m,n)}(z)\right|^2.$$
(66)
where Eq. (46) is multiplied by a factor of $`2`$ to account for the two surfaces at $`y=\pm W/2`$.

With the closed forms of the displacement fields at hand, we can obtain the analytical expression for the attenuation coefficient. We first evaluate $`\overline{\mathrm{\Gamma }}^{(m,n)}`$ from Eq. (42). Since $`T_{iz}^{(m)}0`$, the expression for $`\overline{\mathrm{\Gamma }}`$ reduces to

$$\overline{\mathrm{\Gamma }}^{(m,n)}\left[\rho \omega ^2\left(\varphi _x^{(m)}\varphi _x^{\left(n\right)}+\varphi _y^{(m)}\varphi _y^{\left(n\right)}+\varphi _z^{(m)}\varphi _z^{\left(n\right)}\right)\frac{1}{E}\left(\overline{T}_{xx}^{\left(m\right)}\overline{T}_{xx}^{\left(n\right)}\right)\right]_{y=W/2}.$$
(67)
In addition, putting $`T_{zz}^{(m)}`$ in Eq. (37) at the stress free boundary to zero gives

$$_zu_z^{\left(m\right)}=\sigma _xu_x^{(m)}$$
(68)
so that $`T_{xx}^{(m)}`$ from Eq. (35) simplifies to

$$T_{xx}^{(m)}=E_xu_x^{(m)}.$$
(69)
Now Eq. (66) can be written as

$$\gamma _{n,m}=\frac{\stackrel{~}{g}\left(k_mk_n\right)}{2\rho ^2\omega ^2v_g^{\left(m\right)}v_g^{\left(n\right)}}\left|_{\frac{d}{2}}^{\frac{d}{2}}𝑑z\left[\rho \omega ^2\varphi _i^{(m)}\varphi _i^{\left(n\right)}+Ek_sk_n\varphi _x^{(m)}\varphi _x^{(n)}\right]_{y=\frac{W}{2}}\right|^2,$$
(70)
where the index $`i`$ is summed over $`x,y,z`$. The scattering in the thin plate limit is seen to have two components: the kinetic term, the first term in the $`[]`$ in Eq. (70), which involves all components of the displacement; and the stress term, the second term, which just depends on the longitudinal displacement.

To see how the scattering rate scales with the parameters it is useful to rewrite Eq. (70) as

$$\begin{array}{c}\gamma _{n,m}L=\frac{\stackrel{~}{g}\left(k_mk_n\right)L}{2W^4}\times \frac{W^2\omega ^2}{v_g^{\left(m\right)}v_g^{\left(n\right)}}\hfill \\ \hfill \times \frac{\left|_{d/2}^{d/2}𝑑z\left[\varphi _i^{(m)}\varphi _i^{\left(n\right)}+\frac{Ek_sk_n}{\rho \omega ^2}\varphi _x^{(m)}\varphi _x^{(n)}\right]_{y=\frac{W}{2}}\right|^2}{\left(_{d/2}^{d/2}𝑑z_{W/2}^{W/2}\frac{dy}{W}\varphi _i^{(m)}\varphi _i^{(m)}\right)^{1/2}\left(_{d/2}^{d/2}𝑑z_{W/2}^{W/2}\frac{dy}{W}\varphi _i^{(n)}\varphi _i^{(n)}\right)^{1/2}}.\end{array}$$
(71)
The first factor is a dimensionless measure of the strength of the roughness; the second factor is a dimensionless ratio that depends, through the dispersion relation, only on the geometric ratio $`d/W`$ and the Poisson ratio $`\sigma `$; and the final factor involves integrals over the displacement fields, where we have introduced the explicit normalization factors in the denominator so that we may evaluate the ratio using convenient unnormalized expressions for the displacements.

### III.3 Evaluating the group velocity

As we have seen in Eq. (20), we can avoid evaluating the group velocity appearing in Eq. (70) via numerical differentiating the dispersion curve by instead relating the group velocity to the energy flux in the mode, which in turn can be written as an explicit integral Eq. (21). Thus we need to evaluate the expression (we suppress the mode index in this section )

$$P=\frac{1}{2}\mathrm{Re}\left[i\omega \left(T_{xx}u_x^{}+T_{yx}u_y^{}+T_{zx}u_z^{}\right)𝑑y𝑑z\right].$$
(72)
involving the displacement fields and their derivatives.

In the thin plate limit the $`z`$ components of the stress are small. If we approximate $`T_{zz}=0`$ then expressions Eqs. (3) and (7) can be used to evaluate the $`z`$-component of the strain

$$_zu_z=\frac{\sigma }{(1\sigma )}\left(_xu_x+_yu_y\right).$$
(73)
This can be then used to simplify the in-plane components of the stress

$`T_{xx}`$ $`={\displaystyle \frac{E}{(1\sigma ^2)}}(_xu_x+\sigma _yu_y),`$ (74)
$`T_{yy}`$ $`={\displaystyle \frac{E}{(1\sigma ^2)}}(\sigma _xu_x+_yu_y),`$ (75)
$`T_{yx}`$ $`={\displaystyle \frac{E}{2(1+\sigma )}}(_xu_y+_yu_x).`$ (76)
These expression are used to evaluate the first two terms in the integrand in Eq. (72). The evaluation of the last term in the integrand turns out to depend on whether we are looking at the in-plane or flexural modes, and we now consider each case in turn.

#### III.3.1 In-plane modes

For the in-plane modes in the thin plate limit it is sufficiently accurate to approximate $`T_{zx}0`$, and we can evaluate the remaining terms in $`P`$ with the approximations $`u_x\overline{u}_x`$, $`u_y\overline{u}_y`$ independent of $`z`$. This yields

$$P=\mathrm{Re}\left\{\frac{i\omega Ed}{4\left(1\sigma ^2\right)}𝑑y\left[2(_x\overline{u}_x+\sigma _y\overline{u}_y)\overline{u}_x^{}+(1\sigma )\left(_x\overline{u}_y+_y\overline{u}_x\right)\overline{u}_y^{}\right]\right\}.$$
(77)
#### III.3.2 Flexural modes

For the flexural mode the approximations $`T_{zx}0`$ and $`u_z(x,y,z)\overline{u}_z(x,y)`$ independent of $`z`$ lead to the expressions for the horizontal displacements

$`u_x(x,y,z)`$ $`z\overline{u}_z(x,y),`$ (78)
$`u_y(x,y,z)`$ $`z\overline{u}_z(x,y).`$ (79)
Using these expressions with Eqs. (74-75) shows that the first two terms in Eq. (72) are of order $`d^3`$, i.e. *third order* in the expansion parameter of thin plate theory $`d/W`$. It turns out that to this order, we *cannot* neglect the last term in $`T_{zx}`$, even though all $`z`$-components in the stress tensor are nominally “small”. Indeed comparing the group velocity evaluated from Eq. (72) neglecting the term in $`T_{zx}`$ with those given by numerically differentiating the dispersion curve shows a clear discrepancy. This same problem comes up in deriving the wave equation for the flexural waves

$$\rho d\omega ^2\overline{u}_z=D_{}^4\overline{u}_z.$$
(80)
The term on the left hand side is the mass per unit area times the vertical acceleration, which is given by the integral over the depth of $`_xT_{zx}+_yT_{zy}`$. Clearly the components of $`T_{zi}`$ cannot be neglected completely. Their “smallness” is what leads to the unusual fourth order derivative appearing in this wave equation, with a coefficient again proportional to $`d^3`$.

We have used two methods to arrive at the correct calculation of the energy flux integral for the flexural waves, which is then used to calculate the group velocity for the these waves. The first is to use an improved approximation to the expressions for the in-plane displacements Eqs. (74,75) and a nonzero $`T_{zx}`$ following the approach of Timoshenko T . The second evaluates the energy flux in terms of the vertical displacement and an effective vertical force, and in addition the rotational displacement $`\theta `$ and corresponding torque $`M`$, as is used in the macroscopic derivationLF7  of the wave equation (80). Either of these methods leads to the expression for the energy flux

$$\begin{array}{c}P\frac{1}{2}\omega D\mathrm{Re}\{dy[2k^3\overline{u}_z\overline{u}_z^{}+k(1\sigma )(_y\overline{u}_z)_y\overline{u}_z^{}k(1+\sigma )(_y^2\overline{u}_z)\overline{u}_z^{}]\hfill \\ \hfill +Dk\left[(1\sigma )(_y\overline{u}_z)\overline{u}_z^{}\right]_{y=\frac{W}{2}}Dk\left[(1\sigma )(_y\overline{u}_z)\overline{u}_z^{}\right]_{y=\frac{W}{2}}\}\end{array}$$
(81)
The derivations are displayed in Appendix B. The comparison of the group velocity derived from Eq. (81) and from numerically differentiating the dispersion curve now shows agreement to high accuracy.

## IV Scattering Analysis

The thermal attenuation is calculated from Eq. (70) for normalized mode displacement fields or (71) in general. The group velocity for each mode can be accurately evaluated numerically from the equality $`v_g=2P/\rho w^2`$, with the energy flux $`P`$ given by Eq. (77) for the in plane modes and Eq. (81) for the flexural modes (both expressions are for normalized displacement fields). These are all explicit results in terms of the mode displacements, which are given by Eqs. (48-51) for the in-plane modes, and Eqs. (56,57) for the flexural modes.

Before analyzing the scattering behavior, we first need to have a good understanding of the dispersion relation of the modes, since the scattering rates are strongly dependent on this.

### IV.1 Dispersion relation and group velocity

The dispersion relations for a representative case are shown in Fig. (3). For this example we have used a Poisson ratio of $`0.24`$, and a depth to width ratio of $`d/W=0.375`$, values corresponding to the experimental work of Schwab et al. SHWR00 . As we have discussed, the modes fall into four classes, depending on their parity signatures. We label the lowest mode from each class, the one with zero frequency as the wave number goes to zero, as mode $`0`$, and the modes with successively higher cutoff frequencies in each class as mode $`1`$, mode $`2`$, etc., in that class.

Notice that one of the curves in the figure, the one for the in-plane mode with cutoff frequency $`\omega W/c_t5`$, shows anomalous dispersion with the frequency *decreasing* as the wave number increases up to about $`3W^1`$. (This is actually an even mode, and some higher even and odd modes also show anomalous dispersion.) The dispersion curves for all modes $`n>0`$ have zero slope, and so zero group velocity, at onset. As we will see later, this results in a diverging scattering rate at each mode onset. For the $`n=0`$ modes, as $`\omega 0`$ two of the modes (the compression and torsion modes) have linear dispersion, whilst the other two lowest modes (in-plane and flexural bending modes) exhibit quadratic dispersion. Figure (4) shows the group velocities $`v_g`$ for the four lowest in-plane modes. The group velocity of the bending mode approaches zero as $`\omega 0`$ whilst that of the compression mode becomes constant. The group velocity of the compression mode suddenly drops to $`0.5c_t`$ around $`\omega W/c_t4.6,`$ then gradually recovers and approaches $`0.9c_t`$. These features of the dispersion curve will be reflected in the behavior of the scattering of the waves.

### IV.2 Scattering behavior

We first consider the scattering and reduction of the thermal transport by white noise roughness $`\stackrel{~}{g}(k)=\stackrel{~}{g}(0)`$. This allows us to focus on the role of geometry and the unusual mode structure of the elastic waves in the physics.

In the low frequency limit the dispersion curve and the spatial dependence of the modes take on the simple analytic forms shown in Table 1, allowing us to make analytic predictions for the scattering at low frequencies, and then the thermal conductance at low temperatures. Since only small wave vector scattering is involved in these calculations, the results are true for a general roughness correlation function, providing $`\stackrel{~}{g}(0)`$ is nonzero. The mode structure in Table 1 may be calculated from Eqs. (48-57) taking $`k0`$ or from arguments of macroscopic elasticity theory.

The contributions to the thermal attenuation coefficient in the low frequency limit ($`\omega W/c_t1`$) from the various scattering processes are shown in the Table 2 <sup>3</sup><sup>3</sup>3A more accurate expression for the scattering between the in-plane compression and bending modes is giving by keeping the next order term in $`\overline{\omega }`$ which is $`O(\overline{\omega }^{1/2})`$ for this scattering process. To include this correction multiply the expression in the table by $`(1+\sqrt{\overline{\omega }}/12^{1/4})^2`$.. The expressions take on their simplest form if we introduce the frequency scaled with the velocity of the long wavelength compression mode $`\overline{\omega }=\omega c_E/W`$ with $`c_E=\sqrt{E/\rho }=\sqrt{2(1+\sigma )}c_t`$. The power laws can largely be understood from the prefactor in Eq. (71), $`\gamma _{n,m}\omega ^2/v_g^{\left(m\right)}v_g^{\left(n\right)}`$. The group velocity $`v_g`$ becomes a constant at small frequencies for the compression and torsion modes. Thus the torsion-torsion and compression-compression scattering shows the $`\omega ^2`$ dependence corresponding to Rayleigh scattering in one dimension, and as was found for scalar waves with linear dispersion. On the other hand for the bending modes $`v_g\omega ^{1/2}`$. This has the important consequence that the in-plane bend-bend scattering increases more rapidly at low frequencies proportional to $`\omega `$, and the torsion-bend and compression-bend scattering have an $`\omega ^{3/2}`$ frequency dependence. For the flexural bend-bend scattering the two terms in the braces in Eq. (70) cancel to leading order resulting in smaller scattering $`O(\omega ^3)`$ than given by the prefactor alone. Note that the expressions for the flexural modes involve additional factors of $`W/d`$, so that these modes will be scattered more strongly at a given $`\omega `$ in the thin plate limit. This is because these modes are softer, so that the scattering wave vectors are larger for the same frequency.

Numerical results for the attenuation coefficient $`\gamma _m`$ of the four lowest modes are shown in Fig. (5) for the in-plane and Fig. (6) for the flexural modes. The plot for the in-plane modes in particular shows interesting structure deriving from the complicated dispersion curves of Fig. 3. Much of this structure can be understood from the product of group velocities in the denominator of Eq. (70). In particular there is a square root divergence in $`\gamma _m`$ at the onset frequency of each mode where the group velocity is zero. In addition, the large scattering around $`\omega W/c_t=5`$ derives from the region of anomalous dispersion, since the group velocity is small in this frequency range. The insert to Fig. (5) shows an expanded view of the low frequency behavior, using the results from Table 2 together with the next order correction for the compression-bend scattering. The agreement for the compression mode is very good even up to $`\omega W/c_t3`$, whereas for the bend mode the correspondence is only good for $`\omega W/c_t0.5`$. The scattering for the flexural modes shows generally similar results, Fig. (6) although the behavior is simpler corresponding to the rather featureless dispersion curves. At low frequencies, (insert to Fig. (6)), the scattering of the flexural-bend mode is small, since the intramode scattering is reduced by the cancellation discussed above.

Figure (7) shows the total scattering $`_m\gamma _m`$ for the in-plane modes on a log-log plot, again with white noise roughness. At very low frequencies the scattering varies proportional to $`\omega `$ corresponding to the dominant intramode scattering of the compression mode at low frequencies (Table 2). For frequencies up to $`\omega W/c_t3.5`$, the first nonzero onset frequency of an in-plane mode, the analytic low-frequency expression given by summing the in-scattering expressions from Table 2 (cc, cb, bc, and bb), shown as the dotted line in Fig. 7, gives a good approximation to the full results. At higher frequencies the total scattering increases rapidly, following a general trend proportional to $`\omega ^4`$ (dashed line) together with divergent scattering at each mode onset frequency. The $`\omega ^4`$ power law can be understood as the combination of the explicit $`\omega ^2`$ dependence of Eq. (70), together with two powers of $`\omega `$ coming from the number of modes available for scattering from and to.

### IV.3 Change in the thermal conductance

In the weak scattering limit the change in thermal conductance at low temperatures can be derived directly from the expressions for the scattering at low frequencies. If we write the thermal attenuation coefficient of mode $`m`$ as $`\gamma _mL=A(\omega /\omega _0)^p`$, where $`p`$ is the power law obtained in the low frequency limit and $`\omega _0`$ some characteristic frequency, then the corresponding contribution of the suppression of the thermal conductance from this mode is

$$\delta K_m/K_u=AI_p(T/T_0)^p$$
(82)
with $`T_0=\mathrm{}\omega _0/k_B`$ the corresponding characteristic temperature and $`K_u=\pi ^2k_B^2T/3h`$ the universal thermal conductance. The constant $`I_p`$ can be obtained evaluating the integral

$$I_p=\frac{3}{\pi ^2}_0^{\mathrm{}}𝑑y\frac{y^{p+2}e^y}{(e^y1)^2}.$$
(83)
Thus the power law for the temperature dependence of the depression of the thermal conductivity is the same as the one for the low frequency behavior of the scattering coefficient.

Figures (8) and (9) show the thermal conductance depression scaled with the universal value $`K_u`$ as a function of the appropriate scaled temperature for the lowest in-plane and flexural modes, showing the deviation from the low temperature power laws as the temperature is raised. For the in plane modes we use the characteristic temperature $`T_E=\mathrm{}c_E/k_BW`$ and for the flexural modes $`T_F=\mathrm{}c_Ed/k_BW^2`$. The individual plots are then independent of the geometry. To combine the contributions from the in-plane and flexural modes the ratio $`d/W`$ is needed to relate the two temperature scale factors. In the thin plate limit $`T_F=(d/W)T_ET_E`$.

## V Comparison with Experiment

### V.1 Experimental geometry

Based on the SEM micrograph of the experimental structureSsp , we set the dimensions of the structure in the following way. In the experimental structure of Schwab et al. the thermal pathway was constructed with the shape function $`W(x)=W\mathrm{cosh}(Ax)`$ so that the beam width becomes large and joins smoothly to the thermal reservoirs at the ends, reducing the scattering due to the geometric imperfection at these junctions. Unfortunately this makes the calculation of the behavior of the elastic waves in the beams much harder. However both with and without the scattering off surface roughness, we expect the narrow portion of the beam to dominate the behavior. Thus we simplify the structure and model it as an elastic beam with rectangular cross section of width $`W`$, depth $`d`$, and effective length $`L`$. We estimate the width as the narrowest width of the structure, $`W160nm`$, and $`L=1`$ $`\mu `$m as the length over which the width is approximately constant. The thickness of the material was $`d=60`$ nm. The accuracy of the length estimation is not very critical, since the only length dependence in the scattering rate $`\gamma `$ appears in the combination $`\delta ^2L`$ where $`\delta `$ is the rms roughness which is a parameter of the model, so that any error in the assignment of $`L`$ will just change the value assigned to $`\delta `$. The width $`W`$ on the other hand plays a crucial role, for example determining the frequency cutoffs of the various modes, and so the temperature dependence of the thermal conductivity.

### V.2 Roughness correlation function

Since the nature of the surface roughness on the experimental structure is not known, to fit the experimental data we need a sensible parameterization of the roughness. As a starting point we choose a Gaussian correlation function for the roughness, leading to the spectral density

$$\stackrel{~}{g}\left(k\right)=\sqrt{\pi }a\delta ^2\mathrm{exp}\left[\frac{a^2}{4}k^2\right].$$
(84)
This parameterization of the roughness contains two parameters: $`\delta `$ the rms roughness and $`a`$ the correlation length.

To analyze the data, we first quantify the amount of scattering by subtracting the data of Schwab et al. from the ideal thermal conductance obtained numerically using the “xyz” algorithmNAW97 . Then we attempt to fit the data by adjusting the two parameters $`a`$ and $`\delta ^2L`$.

The inadequacy of Eq. (84) in fitting the experimental data is shown by the low temperature fits in Fig. (10). At these low temperatures only small wave number modes are excited, so that the exponential term in Eq. (84) can be approximated as unity and $`\stackrel{~}{g}\left(k\right)\stackrel{~}{g}\left(0\right)=\sqrt{\pi }a\delta ^2`$. Thus the roughness parameters only appear in the combination $`a\delta ^2`$, and this quantity can be varied to attempt to fit the low temperature region. As seen from the figure, increasing $`a\delta ^2`$ causes scattering that is systematically larger than the experimental data at the low temperatures, while decreasing $`a\delta ^2`$ does not provide enough scattering in the range $`0.1<T<0.2`$ K.

Although there is considerable scatter in the data over the range of the fit, the systematic differences between the predictions and the data lead us to propose a modified form of the roughness correlation that reduces the scattering at small wave numbers

$$\stackrel{~}{g}\left(k\right)=\sqrt{\pi }a\delta ^2\mathrm{exp}\left[\frac{a^2\left(kk_0\right)^2}{4}\right].$$
(85)
A nonzero value of the parameter $`k_0`$ leads to a roughness correlation function that is maximum at a length scale of order $`k_0^1`$, and serves to reduce the scattering at long wavelengths. As mentioned in the introduction section, the same discrepancy (i.e., the overestimation of the scattering at long wavelengths in the theory compared with experiment) was found using the scalar model of the elastic wavesSC00 . The full elasticity theory considered here actually makes the discrepancy worse, since the scattering at small frequencies now is predicted to increase more rapidly at small frequencies than the $`\omega ^2`$ found in the scalar theory, varying as $`\omega ^p`$ with $`p<2`$ for most of the scattering processes, see Table 2.

To fit the data of Schwab et al, we need to determine three parameters: $`k_0`$, $`a`$, and $`\delta `$. As well as looking at the fit by eye, we evaluate the quality of the fit by calculating the mean square deviation of the data from the theory curve over the temperature range up to $`0.4K`$. At higher temperatures many modes becoming excited, and the scattering of individual modes becomes strong, so that our theory is less reliable. Since the onset frequency of the scattering at low frequencies and the initial decrease in thermal conductance with increasing temperature near the onset is mainly determined by $`k_0`$, this parameter is the easiest to determine. We find the value $`k_0W=4.9`$ , rather insensitive to the values of $`a`$ and $`\delta `$. We have also investigated the fit with a delta function noise correlation function (i.e. $`a0`$ in Eq. (85)). With this roughness, to get a maximum in the decrease in the thermal conductance at about $`T0.2K`$ requires $`k_0W`$ values between $`4.5`$ and $`5`$, which agrees with the previous estimate.

The remaining fit parameters $`a`$ and $`\delta `$ are not well determined. Visual inspection of the fit and a plot of the error as a function of $`(a,\delta )`$ shows a reasonable fit over the range $`a/W=5,\delta /W=0.17`$ to $`a/W=6.7,\delta /W=0.32`$. A choice of the roughness amplitude towards the lower end of this range seems most likely physically (e.g., $`\delta /W=0.3150`$ nm is unlikely for roughness generated by chemical etch). We therefore use $`\delta /W=0.2`$ and $`a/W=5.5`$ to generate the curves.

A difficulty of fitting the data is the lack of data points at very low temperatures: it is in this range where only a few modes are involved that we have a very good understanding of the scattering. At higher temperatures many more modes become involved, and the scattering of individual modes becomes strong, so that the second order approximation used in calculating the scattering will not be good. A full test of the theory explaining the reduction in the thermal conductance in terms of the scattering off surface roughness requires more data below a temperature of about $`0.08K`$ for the type of geometry used by Schwab et al., or systems with smaller geometries where the effects can be measured at higher temperatures.

### V.3 Individual mode contribution to the thermal conductance

It is interesting to investigate the contribution to the total thermal conductance of the individual modes with the roughness parameters used to fit the experimental data. This is shown in Fig. (13). The flex-bending mode shows a much smaller contribution to the reduction in $`K`$ at low temperatures for the reason we have already discussed. The modes with nonzero onset frequencies start to contribute significantly above about $`T0.2K`$, and this is the predominant cause for the increase in thermal conductivity above this temperature, since the recovery of the thermal conductance for the lowest mode occurs very slowly.

## VI Conclusion

We have investigated the effect of surface roughness on the scattering of elastic waves in a rectangular beam or waveguide, and the resulting depression of the thermal conductance in the low temperature quantized limit, using full elasticity theory. Our formulation is quite general, but to obtain concrete results we have specialized to the thin-plate limit, which should be a reasonable approximation for many mesoscopic experiments where the depth of the structures is fixed by the epitaxial growth, whilst the width is determined lithographically. The thin plate limit preserves the peculiar features of the elastic waves in the full elastic theory, namely a quadratic dispersion at long wavelengths for two of the low frequency modes, and regions of negative dispersion in the spectra. A robust result is that the low frequency asymptotic dependence of the scattering by unstructured roughness of the modes that propagate at low frequencies (the ones that are important in the low temperature universal thermal conductance) depends on the structure of the modes and the dispersion relation, and is *not* the simple $`\omega ^2`$ dependence of Rayleigh scattering as found in the scalar approximation to the modes. We find different power laws for the various mode scattering processes that can be understood largely from the dispersion relations at: $`\omega `$ for intramode scattering for the in-plane bend mode (the flex-bend intramode scattering is anomalous because of a cancellation between leading order terms, and varies as $`\omega ^3`$); $`\omega ^{3/2}`$ for scattering between the bend modes and the modes with linear dispersion (torsion and compression modes); and the usual $`\omega ^2`$ for the intramode scattering of the modes with linear dispersion. The current experimental data on the suppression of the low temperature thermal conductance below the universal value does not extend to low enough temperatures to provide a good test of these predictions. To investigate this prediction further, it would be interesting to extend the experiments to lower temperatures, or to smaller devices such as carbon nanotubes, where the characteristic temperature scales (when a typical thermally excited phonon has a wavelength comparable with the device dimensions) are higher.

We have used our results to understand the data the data of Schwab et al., who observed a depression of the thermal conductance below the universal value in the temperature range of $`0.1K`$ to $`0.4K`$. Although the scatter in the data is considerable at these low temperatures, the observations seem to show a delay in the onset of the depression scattering as the temperature is raised, beyond what can be fitted with our predictions for unstructured surface roughness. We tentatively resolve this delay by supposing that the surface roughness has a maximum amplitude at some nonzero length scale, which we parameterize by a shifted Gaussian correlation function. Due to the lack of data at low temperatures, a precise determination of the roughness parameters is not possible. However, we do obtain a fit to the data with parameters that do not look unreasonable when compared with electron micrographs of the actual devices.

Our results are based on second order perturbation theory, and the thermal conductance is evaluated assuming the scattering over the length of the device is small. This is a good approximation at low temperatures, but the scattering becomes strong at higher temperatures, particularly for the new modes excited as the temperature is raised, which have a diverging scattering at onset due to the flat dispersion relation here. At higher temperatures multiple scattering and perhaps phonon localization will therefore become important. Kambili et al.KFFL99  and Sanchez-Gil et al. SFMY99  have numerically investigated the these effects in the simplified scalar wave approximation. It would be interesting in the future to extend their work to the full elasticity model.

###### Acknowledgements.

This work was supported by NSF grant no. DMR-9873573. We thank Ruben Krasnopolsky for help on the numerical codes.

## Appendix A Incident and Scattered Fields

Using Green’s theorem we have expressed the displacement field at frequency $`\omega `$ in terms of the surface integral

$$u_q(𝐱)=_S^{}\left[n_j^{}T_{ij}\left(𝐱^{}\right)G_{iq}(𝐱^{},𝐱)n_j^{}u_i\left(𝐱^{}\right)\mathrm{\Gamma }_{ijq}(𝐱^{},𝐱)\right]𝑑S^{}$$
(86)
Eq. (86) involves the integration over a closed surface $`S^{}`$, which we have chosen to be the smooth boundaries together with the cross sections at $`x^{}\pm \mathrm{}`$. In this appendix we show that the integration over the sections at $`\pm \mathrm{}`$ simply yields the incident field $`u_q^{\mathrm{in}}`$, and this allows us to deduce the expression for the scattered field as an integration over the side surfaces. To deduce this result, we first need to derive what are known as reciprocity relations for the elastic modesA .

Let $`𝐮^{(r)}`$ and $`𝐮^{(s)}`$ be the displacement fields for modes $`r`$ and $`s`$ in the ideal beam, and $`𝐓^{(r)}`$, $`𝐓^{(s)}`$ the corresponding stress tensor fields. The modes satisfy the wave equation at frequency $`\omega `$, so that

$$\begin{array}{c}\rho \omega ^2u_i^{(r)}+_jT_{ij}^{(r)}=0\\ \rho \omega ^2u_i^{(s)}+_jT_{ij}^{(s)}=0\end{array}$$
(87)
Multiply the first equation by $`u_i^{(s)}`$ and the complex conjugate of the second by $`u_i^{(r)}`$, subtract the two equations, integrate over a volume of the beam between $`x=x_1`$ and $`x=x_2`$, and finally use the divergence theorem to find

$$_S\left[u_i^{(s)}T_{ij}^{(r)}u_i^{(r)}T_{ij}^{(s)}\right]\widehat{n}_j𝑑S=0,$$
(88)
where the integral is over the surface bounding the volume, consisting of the sides of the beam between $`x_1`$ and $`x_2`$, and the sections at $`x_1`$ and $`x_2`$. The integrations over the sides of the beam are zero by the stress free boundary conditions. For the integration over the sections introduce the explicit $`x`$-dependence $`𝐮^{(r)}=\varphi (y,z)e^{ik_rx}`$ and $`𝐓^{(r)}=\overline{𝐓}^{(r)}(y,z)e^{ik_rx}`$ with $`k_r`$ the wave number of mode $`r`$ at frequency $`\omega `$ etc. Then Eq. (88) reduces to

$$\left(1e^{i\left(k_rk_s\right)\left(x_1x_2\right)}\right)\left[\varphi _i^{(s)}\overline{T}_{ix}^{(r)}\varphi _i^{(r)}\overline{T}_{ix}^{(s)}\right]𝑑y𝑑z=0$$
(89)
and the integral is independent of $`x`$. Unless the prefactor is zero, this shows us that the integral over the section must be zero, and so

$$\left[u_i^{(s)}T_{ix}^{(r)}u_i^{(r)}T_{ix}^{(s)}\right]𝑑y𝑑z=0,k_rk_s.$$
(90)
This is one version of the reciprocity relations.

For our purposes it is more convenient to express the condition for the reciprocity integral to be zero in terms of the group velocity rather than the wave number. To do so, we need to consider the dispersion curves. The condition for the reciprocity integral to be nonzero, $`k_r=k_s`$ for modes $`r,s`$ at the same frequency $`\omega `$, actually implies $`r`$ and $`s`$ are the *same* mode, so that in fact $`v_g^{(r)}=v_g^{(s)}`$. The only other possibility is that $`r`$ and $`s`$ are modes with dispersion curves that cross at frequency $`\omega `$, $`k=k_r=k_s`$. However only modes of different $`y,z`$ parity signatures can cross, and then the integration over the section for these different modes in Eq. (90) is again zero. Thus we can rewrite the reciprocity relation as

$$\left[u_i^{(s)}T_{ix}^{(r)}u_i^{(r)}T_{ix}^{(s)}\right]𝑑y𝑑z=0,v_g^{(r)}v_g^{(s)}.$$
(91)
If $`r`$ and $`s`$ are the same mode, the integral is related to the energy flux and hence to the group velocity (see Eq. (19))

$$𝑑y𝑑z\left(u_i^{(r)}T_{ij}^{(r)}u_i^{(r)}T_{ij}^{(r)}\right)=2i\rho \omega v_g^{(r)}.$$
(92)
We now use Eqs. (91,92) to evaluate the contributions to Eq. (86) from the integrations over the sections at $`x^{}\pm \mathrm{}`$.

Let us first consider $`x^{}\mathrm{}`$. According to Eq. (18) the $`x^{}`$ dependence of the Green’s function pair $`𝐆,𝚪`$ consist of modes $`𝐮_s(x^{})^{}`$ with $`v_g^{(s)}<0`$ since here $`x^{}>x`$ for any finite $`x`$. On the other hand the field pair $`𝐮,𝐓`$ are made up of the incident wave, and waves scattered from the roughness at finite $`x`$, and so consist of modes $`𝐮_r(x^{})`$ with $`v_g^{(r)}>0`$. The integral in Eq. (86) over the section at $`x^{}\mathrm{}`$ is therefore the sum of terms involving $`\left[u_i^{(s)}T_{ix}^{(r)}u_r^{(r)}T_{ix}^{(s)}\right]𝑑y𝑑z`$ with $`v_g^{(r)}`$ and $`v_g^{(s)}`$ of opposite sign. All these terms are zero by Eq. (91), and so there is no contribution from the section at $`x^{}\mathrm{}`$.

Similar arguments apply to the section at $`x^{}\mathrm{}`$. The Green function is made up of modes with $`v_g>0`$. The scattered component of the field $`𝐮`$ consists of modes with $`v_g<0`$, and there is no contribution to the integral over the section from these modes. On the other hand the incident wave $`𝐮^{\mathrm{in}}`$ is mode $`𝐮_m`$ with $`v_g^{(m)}>0`$, and there is the single term with $`v_g^{(n)}=v_g^{(m)}`$ surviving in the sum over modes in the Green function. Using Eq. (92) the integral just gives $`u_q^{(m)}(x)`$. Writing $`𝐮=𝐮^{\mathrm{in}}+𝐮^{\mathrm{sc}}`$ then leads to Eq. (14) in the text.

## Appendix B Energy Flux for Flexural Modes

The classical thin plate approximation of setting $`T_{zi}=0`$ is not sufficient to calculate the energy flux of the flexural modes using the integral Eq. (72). In this appendix we evaluate the correct expression for the energy flux by two different methods, first using the extended thin-plate theory of TimoshenkoT  (see also GraffG ), and then using a method in terms of the energy of plate deformationsLF7  that avoids these difficulties.

In the extended thin plate approximation of Timoshenko the $`z`$-dependence of the in-plane displacements is still approximated as linear

$`u_x(x,y,z)`$ $`z\psi _x(x,y),`$ (93)
$`u_y(x,y,z)`$ $`z\psi _y(x,y).`$ (94)
However, the $`x,y`$ dependence is no longer assumed to be given by the gradient of the mean vertical displacement $`\overline{u}_z(x,y)`$, but by the more general expression

$$\psi =_{}\overline{u}_z+_{}S+_{}\times (\zeta \widehat{z})$$
(95)
introducing the scalar potential $`S(x,y)`$ and vector potential $`\zeta (x,y)`$ defining the corrections to the in-plane strain and rotation. Here $`_{}=(_x,_y)`$ is the horizontal gradient. In addition the vertically averaged stress $`T_{zx}`$ is taken to be

$$T_{zx}\kappa ^2\frac{E}{2\left(1+\sigma \right)}(_x\overline{u}_z+_zu_x)$$
(96)
(with a similar expression for $`T_{zy}`$ given by replacing the subscript $`x`$ with $`y`$ everywhere). Here the “shear correction factor” $`\kappa `$, a number of order unity, is introduced to take into account deviations of the in-plane displacements from the assumed linear dependence on $`z`$T . In the usual thin plate approximation $`T_{zi}`$ are set to zero and $`\psi =_{}w`$, so that $`(u_x,u_y)=z_{}w`$.

With the Timoshenko approximations, the equations of motion for the three components of displacement are now investigated.

The equations of motion for the horizontal displacement lead to an equation relating $`\psi `$ to $`\overline{u}_z`$G

$$\frac{D}{2}\left\{(1\sigma )^2\psi +(1+\sigma )_{}_{}\psi \right\}\kappa ^2\mu d\left(\psi +_{}\overline{u}_z\right)=0$$
(97)
(remember $`D=Ed^3/12(1\sigma ^2)`$, with $`E`$ Young’s modulus, and $`\mu `$ the shear modulus). The inertial terms $`_t^2\psi `$ turn out to be negligible in this equation. Using Eq. (95), Eq. (97) becomes

$$D_{}_{}^2(Sw)\kappa ^2\mu d_{}S+\frac{D}{2}(1\sigma )_{}\times (_{}^2\zeta \widehat{z})\kappa ^2\mu d_{}\times (\zeta \widehat{z})=0.$$
(98)
Taking the vertical curl of Eq. (98) gives

$$\frac{D}{2}(1\sigma )_{}^2\mathrm{\Omega }\kappa ^2\mu d\mathrm{\Omega }=0$$
(99)
with $`\mathrm{\Omega }=\widehat{z}_{}\times \psi =_{}^2\zeta `$ the rotation. For a wave disturbance $`e^{ikx}`$, this gives an exponential dependence on $`y`$, $`e^{\pm \lambda y}`$ with

$$\lambda ^2\frac{2\kappa ^2\mu d}{D(1\sigma )}d^2.$$
(100)
Since $`\lambda ^1dW`$, the rotation will be large only over a boundary layer region with width of order $`d`$ near the edges $`y=\pm W/2`$, where the solution takes the form

$$\mathrm{\Omega }(x,y\pm W/2)\mathrm{\Omega }(\pm W/2)e^{ikx}e^{\lambda |yW/2|}.$$
(101)
The vector potential $`\zeta `$ has a similar solution, so that last two terms in Eq. (98) cancel. This leaves for the scalar potential $`S`$

$$_{}\left(D_{}^2(S\overline{u}_z)\kappa ^2\mu dS\right)=0$$
(102)
which immediately gives

$$D_{}^2S\kappa ^2\mu dS=D_{}^2\overline{u}_z.$$
(103)
(We are only interested in $`_{}S`$ and so do not need to keep track of the arbitrary gradient-free function that could be added to this equation.)

The equation of motion for the vertical displacement isG

$$\kappa ^2\mu d_{}^2S=\rho d\omega ^2\overline{u}_z.$$
(104)
Together Eqs. (103,104) give

$$\rho d\omega ^2(\overline{u}_z\frac{D}{\kappa ^2\mu d}_{}^2\overline{u}_z)=D_{}^4\overline{u}_z.$$
(105)
This is the usual fourth order wave equation, with a small correction term of order $`(d/W)^2`$ (the second term in the brackets on the left hand side). Note that solutions to this equation vary on the long scale of order $`k^1`$ or $`W`$, and not the small scale $`\lambda ^1d`$, so that to a good approximation we have

$`\rho d\omega ^2\overline{u}_z`$ $`=D_{}^4\overline{u}_z,`$ (106a)
$`S`$ $`={\displaystyle \frac{D}{\kappa ^2\mu d}}_{}^2\overline{u}_z.`$ (106b)
The first equation is now the standard fourth order wave equation. The second equation for $`S`$ shows it to be small compared with $`\overline{u}_z`$ by of order $`(d/W)^2`$.
The boundary conditions at the edges are that all stresses are zero, so that in particular at $`y=\pm W/2`$

$$𝑑zT_{zy}=\kappa ^2\mu d(_y\overline{u}_z+\psi _y)=0.$$
(107)
Substituting Eq. (95) into this gives

$$_yS_x\zeta =0.$$
(108a)
Equation (108a) together with Eq. (106b) tells us the size of the $`\zeta `$ correction, which at $`y=\pm W/2`$ takes the value
$$\zeta (x,y=\pm W/2)=\frac{D}{\kappa ^2\mu d}\frac{1}{ik}\left(_y_{}^2\overline{u}_z\right)|_{y=\pm W/2}.$$
(109)
This expression can be simplified using the boundary condition $`T_{yy}=0`$ at $`y=\pm W/2`$, which from Eq. (75) and Eqs. (78,79) gives at $`y=\pm W/2`$

$$_y^2\overline{u}_z=\sigma _x^2\overline{u}_z=\sigma k^2\overline{u}_z,$$
(110)
so that

$$\zeta (x,y=\pm W/2)=\frac{ikD(1\sigma )}{\kappa ^2\mu d}\left(_y\overline{u}_z\right)|_{y=\pm W/2}.$$
(111)
The potential $`\zeta `$ is only large in the boundary layers near the edges where it takes the form

$$\zeta (x,y\pm W/2)=\frac{ikD(1\sigma )}{\kappa ^2\mu d}(_y\overline{u}_z)|_{y=\pm W/2}e^{\lambda |yW/2|}$$
(112)
Thus finally we have expressions for the horizontal displacement field, Eqs. (95) and (93) together with Eqs. (106b) and (112) defining $`S`$ and $`\zeta `$, and Eq. (96). These can be used to calculate the additional contribution to the energy flux coming from the $`T_{zx}`$ term in Eq. (72). (The corrections to $`u_x`$ and $`u_y`$ derived here do not change the contributions from the first two terms in Eq. (72) to the order we require, since these terms are already third order in the small parameter $`d/W`$.)

We therefore need to evaluate

$$T_{zx}u_z^{}𝑑y𝑑z\kappa ^2\mu d𝑑y(_xS+_y\zeta )\overline{u}_z^{}.$$
(113)
Both terms in the integral give contributions at the same order. The first term, coming from the correction to the in-plane strain Eq. (106b), is

$$\kappa ^2\mu d𝑑y(_xS)\overline{u}_z^{}=D(_x_{}^2\overline{u}_z)\overline{u}_z^{}𝑑y.$$
(114)
The second term in the integrand is only large in the boundary layer region near the edges and from Eq. (111) evaluates to the edge contributions

$$\kappa ^2\mu d𝑑y(_y\zeta )\overline{u}_z^{}=ikD(1\sigma )[(_y\overline{u}_z)\overline{u}_z^{}]|_{y=W/2}^{y=W/2}.$$
(115)
Combining these expressions for Eq. (113) with Eqs. (77) together with (78) and (79) yields the final expression

$$\begin{array}{c}P\frac{\omega kD}{2}\mathrm{Re}\{dy[2k^2\overline{u}_z\overline{u}_z^{}+(1\sigma )\left(_y\overline{u}_z\right)\left(_y\overline{u}_z\right)^{}(1+\sigma )\left(_y^2\overline{u}_z\right)\overline{u}_z^{}]\hfill \\ \hfill +\left[(1\sigma )\left(_y\overline{u}_z\right)\overline{u}_z^{}\right]_{y=W/2}^{y=W/2}\}.\end{array}$$
(116)
which is identical to Eq. (81).

An alternative approach to calculate the energy flux is to use the expression for the energy of distortions of the plate evaluated using the lowest order expressions Eqs. (74,75) and (78,79)LF7

$$F=\frac{1}{2}D\left[\left(_{}^2\overline{u}_z\right)^2+2(1\sigma )\left\{\left(\frac{^2\overline{u}_z}{xy}\right)^2\frac{^2\overline{u}_z}{x^2}\frac{^2\overline{u}_z}{y^2}\right\}\right]𝑑x𝑑y.$$
(117)
It turns out that the higher order corrections discussed above are not needed in this expression, and so we can derive the energy flux without these difficulties. The functional derivative of $`F`$ with respect to $`\overline{u}_z`$ yields the vertical force per unit area in the interior of the plate, which can be used to derive the fourth order wave equation, as well as expressions for the energy flux into the plate across the boundaries. The latter expressions give us the result for the energy flux along the beam

$$P=\frac{1}{2}\mathrm{Re}\left\{i\omega \left[M_x\theta _x^{}+V\overline{u}_z^{}dy+(F_c\overline{u}_z^{}|_{y=W/2}+F_c\overline{u}_z^{}|_{y=W/2})\right]\right\}$$
(118)
where

$$V=D_x[_x^2\overline{u}_z+\left(2\sigma \right)_y^2\overline{u}_z]$$
(119)
is the effective vertical force that couples to the vertical displacement $`\overline{u}_z`$,

$$M_x=D\left(_x^2\overline{u}_z+\sigma _y^2\overline{u}_z\right)$$
(120)
is the torque that couples to the angular displacement $`\theta _x=\overline{u}_z/x`$, and

$$F_c(y=\pm W/2)=\pm 2D(1\sigma )_{xy}^2\overline{u}_z|_{y=\pm W/2}$$
(121)
is a vertical force localized at the edges of the plate.

Substituting Eqs. (119-121) into Eq. (118) gives

$$\begin{array}{c}P=\frac{1}{2}\mathrm{Re}\{\frac{\left(i\omega D\right)}{2}[dy(_x^2\overline{u}_z+\sigma _y^2\overline{u}_z)(_x\overline{u}_z)^{}+dy[_x^3\overline{u}_z+(2\sigma )_x_y^2\overline{u}_z]\overline{u}_z^{}\hfill \\ \hfill 2(1\sigma )\left(_x_y\overline{u}_z\right)\overline{u}_z^{}|_{y=W/2}+2(1\sigma )\left(_x_y\overline{u}_z\right)\overline{u}_z^{}|_{y=W/2}]\}\end{array}$$
(122)
Evaluating $`_x=ik`$, and using integration by parts, we again get Eq. (81).