On the link between mechanics and thermal properties: mechanothermics

We report on the theoretical derivation of macroscopic thermal properties (specific heat, thermal conductivity) of an electrically insulating rod connected to two reservoirs, from the linear superposition of its mechanical mode Brownian motions. The calculation is performed for a weak thermal gradient, in the classical limit (high temperature). The development is kept basic as far as geometry and experimental conditions are concerned, enabling an almost fully analytic treatment. In the modeling, each of the modes is subject to a specific Langevin force, which enables to produce the required temperature profile along the rod. The theory is predictive: the temperature gradient (and therefore energy transport) is linked to motion amplitude cross-correlations between nearby mechanical modes. This arises because energy transport is actually mediated by mixing between the modal waves, and not by the modes themselves. This result can be tested on experiments, and shall extend the concepts underlying equipartition and fluctuation-dissipation theorems. The theory links intimately the macroscopic size of the clamping region where the mixing occurs to the microscopic lengthscale of the problem at hand: the phonon mean-free-path. This clamping region, which is key, has received recently a renewed attention in the field of nanomechanics with topical works on"phonon shields"and"soft clamping". We believe that our work should impact the domain of thermal transport in nanostructures, with future developments of the theory toward the quantum regime.


I. INTRODUCTION
Context -Solid matter is constituted of atoms located at specific places in space. These positions are defined by the interactions among these atoms, whatever they might be (e.g. covalent, ionic, van der Waals). While these bonding forces are actually rather difficult to describe from first principles (especially at defects like e.g. grain boundaries in crystals, surfaces, voids) [1], they result in an effective potential V ( r i ) experienced by each atom (located at r i = {x i , y i , z i }) that keeps the piece of matter together: which means ∂V /∂x i = 0 and similarly for y i , z i at the rest positions r i,0 [2].
At equilibrium, atoms oscillate weakly around their rest positions: this is Brownian motion, understood as the response of the solid to the energy supplied randomly by the thermal baths to which it is inevitably connected (let it be the surrounding air and mostly the support that holds the object, at temperature T ) [3]. This agitation is conveniently described in terms of phonons, the collective mode excitations arising from the coupling between atoms (the proper eigenstates of the system). At lowest order, this coupling can be seen as a collection of restoring forces with effective spring constants ∂ 2 V /∂x i ∂x j (with all permutations x i ←→ y i , z i and x j ←→ y j , z j ); at higher orders, the anharmonicity of the potential V leads to the finite lifetime of these excitations [2,3]. Modern numerical methods which take into account interactions over a large number of atoms, namely Molecular Dynamics Simulations (MDS), can be used to calculate thermal properties of widely used materials, like e.g. crystalline Silicon [4]. Considering solids with no free electrons (insulators or intrinsic semiconductors at low enough T ) and with no extra degrees of freedom (like e.g. a magnetic moment), phonons lead at the macroscopic scale to both the description of thermal equilibrium, with the definition of the specific heat C v (at constant volume), and to the description of heat transport with the definition  The microscopic atomistic theory enables to built on one side elementary excitations (phonons) and on the other one elasticity theory (stress vs. strain relations, i.e. continuum mechanics). The former leads at the macroscopic scale to the bulk thermal properties: specific heat and thermal conductivity. It should therefore be possible to derive these properties from the latter. This is the subject of the manuscript (the " ? " marked arrow above).
Exerting a small force onto a solid results into a small reversible deformation of this body, of amplitude proportional to the force (the linear limit). Again, this can be seen as a consequence of the effective springs that link all the atoms together. One can then derive elasticity theory by averaging together these interactions over a small volume δτ , and transforming sums over atoms into integrals over space. The procedure is exact, provided this volume is small enough from a macroscopic point of view (which requires the inter-atomic interactions to be short-ranged), but still large enough to host a decent number of atoms; besides, deformations have then to occur on lengthscales larger than this volume, thus much larger than the inter-atomic distance (i.e. long wavelengths). This leads to the definition of the stiffness tensor [c ijkl ] that relates stresses to strains (with each index i, j, k, l referring to the 3 directions in space x, y, z) [2]. As for thermal properties, MDS can be used to evaluate the stiffness coefficients c ijkl of relevant materials, e.g. Silicon [5]. On the macroscopic scale, elasticity theory enables to describe how an object will distort when driven by an arbitrary DC or AC force applied onto a region of it. In the framework of the linear theory discussed here, this is performed by defining the normal modes of the structure, i.e. the standing wave motional states that the object can sustain, and then performing the proper linear superposition [2].
The atomistic description of solids thus leads both to the definition of macroscopic thermal and mechanical properties, Fig. 1; these are just two facets of the same underlying physics (i.e. the inter-atomic potential V ). And indeed, thermodynamics can be probed at the level of an individual mechanical normal mode (as opposed to the whole object). Brownian motion of micromechanical beams has been reported using various systems and techniques [6][7][8][9], from room temperature down to millikekvin temperatures. It is detected as peaks in the motion spectrum at the frequencies of the mechanical modes ω n , with widths Γ n the mechanical damping rates, and areas proportional to the stored elastic energies 1 2 k n < x 2 n >= 1 2 k B T , with x n (t) the amplitude of motion of mode n having an effective spring constant k n . k B is Boltzmann's constant and < · · · > the ensemble average. This result is known as the equipartition theorem, and applies when in-equilibrium [10]. The mathematical formalism describes Brownian motion as due to a stochastic force F n , the so-called Langevin force, that acts upon mode n. It is due to the thermodynamic bath (at temperature T ) that generates the damping Γ n ; there is therefore a link between F n and Γ n , which is called the fluctuation-dissipation theorem [10].
Experiments on single modes are also performed out-ofequilibrium [11,12]. In this case, measured properties (local temperature along the beam, mode spectrum) do depend on the heat flow imposed in the structure. The situation is thus much more complex than for the in-equilibrium case (where everything is defined through the single parameter T ), and heat flows are usually extremely large in the experimental realizations, which adds another degree of complexity [13]. New thermodynamics concepts that describe the steady-state have to be developed beyond standard equipartition and fluctuation-dissipation theorems [11,13].
Developing these new concepts is a very difficult task, and the present manuscript proposes an original theoretical approach to these questions, drastically different from existing literature. From Fig. 1 and what has been stated above about single-mode thermodynamics, it appears that there should be a method to derive the macroscopic thermal properties from a pure mechanistic approach. In order to do so, we shall consider a rather simplified textbook problem which is essentially analytically solvable: a cylindrical rod made of a single isotropic, homogeneous and simple material (e.g. non-magnetic) with no conduction electrons, connected to two reservoirs in vacuum (see Fig. 2 below). For the sake of concreteness, we shall present numerical estimates using properties of polycrystalline Silicon. The dimensions will be assumed large enough to be fairly macroscopic, with respect to the atomic size and lengthscales related to phonons (see discussion below). We will limit ourselves to the treatment of linear response for both thermal and mechanical properties. Finally, the average temperature T will be assumed high enough so that everything can be treated in the framework of classical physics.
Structure of the paper -In Section II, we will first remind the reader about fundamental theoretical aspects related to phonons, thermal equilibrium and thermal transport. Basic parameters will be introduced, and the temperature profile T (z) along the rod presented. Section III gives the mathematical formalism that defines energies (kinetic, potential, friction) in the rod in terms of its normal modes. These are obtained from the Pochhammer-Chree waves [14], which are described in the Appendix A for the sake of completeness. The simple case of in-equilibrium conditions is treated first in Section IV, introducing the standard equipartition and fluctuation-dissipation theorems. Section V then treats the case where a small thermal gradient is imposed, with a weak heat flow. The new implications for equipartition and fluctuation-dissipation theorems are then described in Section VI. In Conclusion, we discuss these findings and how they could be confronted to experiments.
Fundamental outcome of the theory -The modeling predicts nonzero amplitude cross-correlations between the modes, proportional to the heat flow, with a specific dependence to mode number. This actually means that energy transport is not carried by the modes themselves, but by their 2-wave mixing. The theory leads to a link between a microscopic lengthscale, namely the phononic mean-freepath, and a macroscopic one which characterizes the size of the clamping zone. This clamping zone is a key element of the modeling: it is where correlations between modes of nearby wavevectors are created, leading to nonzero mixing amplitudes. The full model is a linear theory, even though nonlinear processes within the clamp, which are responsible for correlations, are key.

A. Collective modes: phonons
Let us consider a bulk piece of a homogeneous material, at temperature T . We shall restrict this introductory discussion to the simplest types of solids: those with crystalline order. The situation of amorphous matter is much more involved and shall not be described here [15]; we will just comment on it at the end of this Subsection when addressing crystalline defects.
The elementary excitations related to motion, namely the phonons, are waves that propagate in the crystal [2,3,16]. They have dispersion relations ω( q) that depend on the direction of propagationq = q/ | q|. For small enough wavevectors | q| compared to the size of the first Brillouin zone (i.e. long wavelengths), these are linear: ω( q) = c n,q | q| with c n,q the wave (phase) velocity in the directionq, for mode n. In each direction, the motion sustains 3 types of waves: one longitudinal (n = l, the vibration is alongq), and two trans-  Figure 2: a) Schematics of the problem at hand (not to scale). A cylindrical rod is clamped between two thermal reservoirs, in vacuum; in inset, the relevant lengthscalesλ mf p ,λT are depicted qualitatively, together with microscopic defects (magenta dots, internal degrees of freedom). Geometrical parameters are discussed in the text; body forces in the rod due to in-built stress (appearing e.g. if one pulls on the clamps) are considered only in Appendix A. b) Theoretical temperature profile along the rod axis (under a heat flowQ); T is uniform within a section of the length L, but has a more complex pattern within the clamps (in dashed, see text).
verse ones (n = t 1 , t 2 with vibration perpendicular toq). The velocities of t 1 and t 2 can be degenerate, but are in principle distinct. When | q| gets close to the border of the first Brillouin zone (i.e. small wavelengths), the dispersion relations start to flatten out. Besides, if the crystallographic structure contains more than one atom per unit cell, phonons also support traveling modes at high frequencies, which frequencies do not extrapolate at zero for | q| → 0: these are the optical branches.
These concepts can be illustrated with Silicon, see e.g. Refs. [17][18][19] providing comparison between theory and experimental data. Monocrystalline Silicon has a fcc structure with 2 Si atoms in the cell; the lattice parameter a is about 0.54 nm, which is the fundamental lengthscale that fixes strict bounds to the modeling (see discussion at the end of this Subsection). Below about 4 THz (i.e. corresponding to 200 K which is "almost" room temperature), the dispersion relations are reasonably linear. This is the regime we shall concentrate on for the work presented in this manuscript.
In thermal equilibrium, phonons distribute in energy ac-cording to Planck's black body law [2,3]: expressed in (Joules/m 3 )/(Rad/s), with Plank's (reduced) constant.c is the average speed of sound in the solid, averaged over all c n,q values. This function is plotted on Fig. 3, withc ≈ 5 000 m/s for Silicon and T = 200 K. It is strongly peaked around ω T ≈ 2.821 (k B T )/ , which corresponds to an average wavelength ofλ T ≈ 2.227 c/(k B T ) if we consider a strictly linear dispersion law. This is our first relevant lengthscale known as the dominant phonon wavelength [2,3]. At T = 200 K, it essentially means that motional energy stored in phonons corresponds to distortions occurring over lengthcales of the order ofλ T ∼ 0.4 nm, i.e. about the interatomic distance. Summing up the energy stored in the whole phonon bath, and deriving it by temperature, one defines the specific heat C v (at fixed volume). It can be written [2,3,16]: where D(ω) is the density of states, i.e. the number of phononic states available between ω and ω + dω per unit volume [in m −3 (Rad/s) −1 ]. The function D(ω) is directly obtained from the dispersion relations ω( q); in practice this displays a rather complex structure with strong peaks. A particularly useful simplified treatment has been proposed by Debye [16], considering a strictly linear dispersion law with average velocityc and replacing the first Brillouin zone by a sphere of radius q D . This radius is chosen such that the total number of degrees of freedom is preserved: having N atoms moving in 3D space, this makes a total of 3N modes. One obtains: with Θ ω D (ω) the step function (1 for ω < ω D and zero above; ω D =c q D ). Eq. (2) is plotted in Fig. 4 in units of nk B (n being the number volumic density), with T D = ω D /k B = 620 K (1/q D ≈ 0.62 Å), the value corresponding to Silicon [16]. The calculated curve is particularly close to experimental data [20], despite the simplifications used. At high temperatures, it saturates at 3: this is known as the Dulong-Petit law, a limit reached when all modes are equally populated [2,3,16]. This is the "classical limit" which is relevant to the present paper. Around 200 K, in Si C v is already of the order of 65 % of this value.
Phonons are collective excitations that behave in the harmonic approximation as a gas of free particles. But even in an absolutely ideal crystal (free of any defects), they experience some interactions among themselves. These interactions are due to the anharmonicity of the actual potential V that describes the solid. Third order terms of the type ∂ 3 V /∂x i ∂x j ∂x k (with all permutations x i ←→ y i , z i ; x j ←→ y j , z j ; x k ←→ y k , z k ) in the Taylor expansion of V lead to 3-phonon processes. Similarly, the fourth order terms lead to 4-phonon processes, and so on. Because of the conservation of energy and momentum, selection rules also apply to these processes. The scattering rates can be calculated using perturbation theory and the Fermi Golden rule, see e.g. Refs. [2,3] for detailed discussions. These  Debye calculation of the Silicon specific heat Cv, with TD = ωD/kB = 620 K. At high temperature, it saturates at 3nkB with n the number volumic density (Dulong-Petit law, horizontal line).
mechanisms are known as inelastic scattering, since they redistribute energy among the phonon gas. They are strongly temperature-dependent [2,3].
Considering a more realistic situation, phonons also scatter from structural defects of the crystal. These are crystallographic mis-stackings, like vacancies (point-like defects), dislocations (line defects) or grain boundaries (surface defects). Even with a chemically pure solid, there is a natural isotopic variation within the material. This is the isotopic replacement of 28 Si (most abundant) with 29 Si and 30 Si in our example, at random places of the lattice. It creates mass point-like defects off which phonons can scatter [2]. These defects do not possess any internal degrees of freedom to which energy can be transferred: thus these processes are elastic, and only redistribute momentum. They are in principle temperature-independent, because they depend only on geometric parameters [2,3].
All these (independent) mechanisms sum up to define a total phonon scattering rate Γq ,n , for a given dispersion branchq, n. Averaging over all n,q one can defineΓ, and within the linear dispersion limit, an average mean-free-path λ mf p =c/Γ using the average velocityc. This is our second relevant lengthscale: it represents the typical distance over which motional energy is homogenized (in both ω and q) over the phonon bath.
Let us consider the high-temperature limit whereλ T < λ mf p and both are much smaller than the size of the solid. In this case, at any point r within the object we can define a phonon number distribution n T ( r) which depends on a local temperature T ( r): this is known as the diffusive limit [2,3]. There is local equilibrium, while the global solid might be out-of-equilibrium if a heat flow is imposed. This heat flow can be described by the Boltzmann transport equation [16]. After some manipulations, the heat current j E (energy flow density carried by phonons, in W/m 2 ) can be written in the form [2,3,16]: with κ = 1 3 C vcλmf p the thermal conductivity of the material. For Silicon, κ ≈ 200 W.m −1 .K −1 around 200 K [21], which leads toλ mf p ≈ 60 nm (the typical size at which local equilibrium is realized). Eq. (4) is the well-known Fourier's law, describing thermal transport at the macroscopic scale, which will be used in the next Subsection. Since equilibrium is reached at a local scale, the thermal properties appear isotropic on the scale of the whole solid.
The whole point of this Subsection is to identify the conditions which are targeted by our modeling, and what underlying limitations have to be considered. For the sake of simplicity, we shall consider in the following a polycrystalline material with a grain size large enough such that thermal and mechanical properties are well defined over it. Assuming the grains small enough compared to the size of the solid, and randomly oriented, the macroscopic mechanical properties reduce to averaged, and thus isotropic, quantities: the stiffness tensor [c ijkl ] is defined by introducing only 2 independent parameters, see Section III. The size of grains is thus our last relevant lengthscale.
In polycrystalline Aluminum thin films used for micromechanics or superconducting quantum bits, defects are known to impact strongly device properties [22,23]. These are Two-Level-Systems (TLS) that reside in the amorphous Aluminum oxide: atoms, or groups of atoms that can occupy different positions within the solid which are nearly equivalent energetically. The switching between these positions can be induced by tunneling or thermal activation, and the overall contribution of such defects generates a rather flat in energy distribution of low-energy states to which strain can couple [24]. This means that in a glass, phonons will behave essentially like in a crystal at long wavelengths (where the actual atomic order does not really matter), but they will be much less well-defined: their scattering rateΓ is strongly impacted by their coupling to TLSs [15,24]. Thus for a nonperfect material, there are internal degrees of freedom which can absorb or emit energy from/to the phonon gas ( Fig.  2 inset). This dissipative channel is clearly distinct from the phonon-phonon interaction already introduced. However, since defects interact (and thus thermalise) through phonons, and since these phonons obey local equilibrium, they shall be described by the same T ( r) profile.
For intrinsic Silicon, even at room temperature the thermal population of charge carriers remains low enough to be completely ignored [21]. This is obviously not the case for a metal, like e.g. Aluminum, were electron-phonon coupling impactsΓ. Besides, thermal conductivity is then dominated by electronic transport; this is a completely different regime than the one discussed here [2,16].
But there is a principle limitation in the description we do of collective modes: we assumed the dispersion relations to be strictly linear (i.e. the Debye approach). Obviously, the simple picture we have developed fails to represent reality at small (i.e. atomic) wavelengths and high energies. At that scale, even the hypothesis of local equilibrium fails; phonon frequencies are higher thanΓ, meaning that collisions are not frequent enough to thermalise over a typical period of the propagating mode. However, numerical estimates (i.e. the Debye specific heat) are in rather good agreement with experiments. For this reason, we shall use in the following the same pragmatic, and simplified approach for mechanical properties: elasticity theory will be assumed exact down to the smallest scales, and we shall bound the number of normal modes by the same argument as for the specific heat (imposing the total number of degrees of freedom to be 3N ).
The relevance of the modal description also depends on the damping each mode experiences, a point discussed explicitly in Subection III D and Section IV. As a matter of fact, local thermal equilibrium will be assumed within this extrapolation. This is a sort of "mechanical adaptation" of Debye's hypothesis.

B. Macroscopic description: Fourier's law
The simple geometry we consider is depicted in Fig. 2. A cylindrical rod of a homogeneous and isotropic material is connecting two thermal baths in vacuum. Its length is L and radius R. On each side, a small region of length ε L represents the clamping zone where the radius monotonically increases in order to adapt to the bulk heat reservoirs. These small portions of the rod will receive particular attention in Subsection III D below. For the time being, we shall treat them on the same footing as the bulk of the rod. Indeed, everywhere in the bar the hypotheses described in the previous Subsection do apply: there is thus a well defined temperature at every point T (r, θ, z), and a well defined thermal conductivity κ(T ). For symmetry reasons, T does actually not depend on θ.
Equilibrium writes, in steady state with no internal heat load div( j E ) = 0. Using Fourier's law Eq. (4) we obtain: In the above, we assumed the thermal gradient along the rod to be very small such that κ can be treated as a constant. For the section of the rod of length L and constant radius R, the problem is fairly easy to solve: the temperature is homogeneous across the rod (no radial heat flow), and depends linearly on the z abscissa. However, within the clamps it is much more complex (and outside of the scope of the paper): as the beam flares-out on the right of Fig. 2 a), the energy flow spreads over a larger cross-section while the total flux is preserved. This means that the temperature gradient along z should be reduced, and a local r-dependent gradient will appear; the reverse is true for the left side clamp. The profile T (0, z) on the axis of the rod is qualitatively drawn in Fig. 2 b).
We need now to specify concrete boundary conditions. The thermal bath on the left is maintained at T c , homogeneously over its entire (and very large) volume. The right bath is at T h , and we arbitrarily chose T h > T c . On the inner-sides of the clamps, the temperatures at each extremity of the homogeneous section of the rod are defined as T h,− and T c,+ (see Fig. 2). One of the key assumptions of the modeling is that the clamps are almost ideal, which implies T h,− ≈ T h and T c,+ ≈ T c . The temperature profile thus reads (for 0 < z < L): with ∆T = T h − T c T avg , and T avg = (T h + T c )/2, see Fig. 2 b). The heat flowQ = j E d S across a section of the rod is then defined as: the − sign reminding that energy flows from the hot side to the cold one. We recognize the standard expression for thermal conductance G = κ πR 2 /L. In the high-temperature limit where the specific heat is constant C v = 3nk B , the total motional energy density stored in the beam (defined per unit length) can be written as: having made use of n = N/(πR 2 L). This last expression is a straightforward implication of the local equilibrium assumption. Eqs. (6,7,8) completely describe thermal properties at the macroscopic level, and are the basic inputs for the mechanicsbased theory we develop below.

A. Propagating acoustic waves
We shall start the Section by recalling basics of continuum mechanics. The distortion of a solid body is described by the displacement field u( r, t), which represents by how much any point r moves under a given solicitation. This solicitation is due to forces acting upon the object, which propagate inside it as stresses described by the tensor [σ ij ( r, t)]. The strain tensor [ kl ( r, t)] is constructed from the first derivatives of the displacement vector u with respect to space coordinates, keeping only the symmetric components. Elasticity theory simply relates the two tensors through the stiffness [c ijkl ]: where the indexes ijkl run each through the 3 directions of space x, y, z [2]. [σ ij ( r, t)] and [ kl ( r, t)] are thus 3 × 3, while [c ijkl ] contains 81 coefficients. But this complexity can be simplified drastically by taking into account the symmetries of the problem at hand. For an isotropic material, only 2 components will be distinct. In the following, we chose to describe mechanical properties with the Young's modulus E Y and the Poisson ratio ν. The former represents the elasticity of the rod when subject to a force along its axis, while the latter quantifies by how much it expands in the radial direction when it is pressed upon. With a homogeneous solid, [c ijkl ] is also r independent. Within the "mechanical adaptation" of Debye's hypothesis mentioned in the previous Section, we shall also consider that [c ijkl ] is strictly constant up to the highest frequencies ω: this means, equivalently, that it is also t-independent. Note that by construction this modeling fails to reproduce the high frequency optical branches. These modes are phenomenologically accounted for in the prolongation of the elasticity law beyond the atomic size, similarly to Debye's phononic treatment.
Newtonian dynamics brings the fundamental equation of motion for the field u( r, t) in the absence of body forces [2]: which is the basis for acoustics, with ρ the mass density of the material, and ∆ the vector Laplacian. We assume R L so that the rod can be treated as infinite: then the displacement profile should display the same shape for all z, t. This means that we can separate the variables: written in cylindrical coordinates. U {η} represents a propagating solution along z, while the φ i{η} (adimensional, with i = r, θ, z) describe the shape of the distortion within a section. The set {η} corresponds to the parameters that will distinguish the different existing solutions (see below). Note the z-derivative in Eq. (13) [see Appendix A]. The mathematical formalism of the next Subsections is built on the assumption that all the existing normal modes within the rod are represented, at least within our "Debye" simplification. We shall therefore spend some time on describing how these modes are constructed. The solutions are found by injecting Eqs. (11 -13) into Eq. (10). The fundamental equation of motion is satisfied with: with K and W two constants. Propagating waves require W = −ω 2 (which we chose ω > 0 with no loss of generality), and can be classified in 3 categories (solving generically for K ∈ C) [14]: • Evanescent, K = +k 2 e , k e ∈ R, The real-valued solutions for the two first ones are: Evan. U {η},0 (z, t) = U 0 exp(±k e z) cos(ωt), (18) or U {η},π/2 (z, t) = U 0 exp(±k e z) sin(ωt), (19) having defined explicitly the two equivalent timedependencies with orthogonal phases 0, π/2. U 0 is a (yet arbitrary) amplitude, and the sign ± defines in which direction the wave travels or decays. For the mixed solutions we have: (20) or U {η},π/2 (z, t) = U 0 exp(±k e z) sin(±k p z − ωt), (21) with the complex solutions built as U {η},0 ± i U {η},π/2 . In the above expressions, we shall keep only solutions with k e k p > 0 for which the mixed wave travels and decays in the same direction (exponential growth being unphysical).
The functions φ r{η} (r, θ), φ θ{η} (r, θ) and φ z{η} (r, θ) that satisfy Eq. (10) are described in Appendix A; the impact of in-built stress is discussed therein. They are regrouped in 3 families: torsional (T), longitudinal (L) and flexural (F)   which are known as Pochhammer-Chree waves [14]. Similar results can be obtained for a rectangular beam (with a rather involved formalism describing the φ i{η} functions, i = x, y, z), see Ref. [25]. For torsional and longitudinal modes, a single index m ≥ 0 describes all the waves of the family; flexural waves require two indexes n > 0, m ≥ 0, and are represented in two degenerate sets (defined by a rotation angle θ 0 = 0, +π/2). Having chosen a wave within a family from its index(es), and having chosen a type of wave (traveling, evanescent, or mixed), solving the dynamics equation generates a relationship between angular frequency ω and wavevectors(s) k p , k e : this is a dispersion relation. All dispersion relations can be represented in 3D space with vertical axis ω and abscissae the complex plane k p + i k e (only first quadrant k e , k p > 0 since all properties are propagationdirection independent). This generates rather complex plots that require to be commented [14,[26][27][28].
We write ω {m,T} (k p ), ω {m,T} (k e ) for the traveling and evanescent torsional dispersion relations, and similarly ω {m,L} (k p ), ω {m,L} (k e ) for longitudinal, with index m. Flexural dispersion relations are designated by ω {n,m,F} (k p ), ω {n,m,F} (k e ) following the same conventions. Only the lowest torsional and longitudinal waves (i.e. m = 0) do not support evanescent propagation. Traveling wave dispersion relations are drawn in the (k p , k e = 0) plane, while evanescent ones are drawn in the (k p = 0, k e ) one. For a given wave, these two join at the vertical axis, when k p = k e = 0. They form what we shall call a branch, that can be represented by flattening the two planes onto the same page, see Appendix A and schematics in Fig.   5. Only the 3 lowest branches (m = 0 torsional and longitudinal, n = 1, m = 0 flexural) go all the way to ω = 0 for k p = 0. All other branches have a finite frequency ω = 0 at k p = ke = 0 (as is the case in Fig. 5).
The zoology of these modes needs to be discussed, since our modeling requires an exhaustive description of the rod's mechanical solutions. For longitudinal and flexural waves, branches with different index(es) connect on the evanescent side in a complex manner which depends strongly on the Poisson ratio ν [26,28]. They are also interconnected by mixed wave solutions (pieces of branches with both k p = 0, k e = 0 that are not contained in one of the above mentioned planes), which originate and end up at local minimas of the branches (dots in Fig. 5) [14,26]. Actually, most of the evanescent solutions possess also pieces of branches which are not directly connected to the traveling solution, but are linked to the solution network only by other branches or mixed wave solutions, see Refs. [14,26,28]. We shall call them secondary branches, as opposed to the main branch which connects to the k p plane. Mixed waves and secondary branches are not drawn in Appendix A, because they are not relevant to our discussion (see Subsection III B below): the key argument being that they always link a given main branch to lower ones, or directly to the ω = 0 plane [14]. Standard relevant configurations are thus shown in Fig. 5 (for branches strictly higher than the 3 lowest ones). At large k p 1/R, the dispersion relations are always linear. On the evanescent side, the branch is always limited in k e (with a rather weak span of available ω values). At connecting points, the group velocity dω/dk e is either zero or infinite.
A solution on a given branch is thus described by a set of parameters {η}, as introduced above. It regroups the index(es) defining the branch plus the wavevector (e.g. {m, L, k p } for propagating longitudinal wave m of wavevector k p , or {n, m, F, k e } for evanescent flexural wave n, m of wavevector k e ).

B. Normal modes
The normal modes of the rod are standing waves obtained from the previous propagating solutions. These are due to the boundary conditions chosen at each end of the rod. We should thus define these (idealized) boundary conditions for the amplitude parameter that describes each of the wave families (see Appendix A for details on this parametrization).
For torsional waves, the amplitude of motion of a mode is defined from the tangential displacement at the periphery u θ (r = R, θ, z, t). By construction, this means that φ θ{η} (r = R, θ) = 1 for any set {η}. The boundary conditions then naturally write: for any θ and t. Similarly, the longitudinal waves are characterized by the motion amplitude in the z-direction on the axis, u z (r = 0, θ, z, t), meaning φ z{η} (r = 0, θ) = 1 for any {η}. The boundary conditions become: Eqs. (11)(12)(13) are formally equivalent under the exchange: so that we can still seek solutions for longitudinal motion using the same Eqs. (16)(17)(18)(19)(20)(21); the boundary condition Eq. (23) is thus formally equivalent to Eq. (22). The flexural waves amplitude is more complex to define: indeed, the motion is characterized through a planar dis- The boundary conditions should then stipulate that u F vanishes on the two extremities: which involves now both u r and u θ . Eqs. (22,23,26) are referred to as idealized "weak clamping", which limits as little as possible the displacement field u at extremities. This modeling will be further discussed in Subsection III D.
Let us consider now that energy is fed into the rod at a given frequency ω, toward a specific branch. The different situations that can be encountered are depicted in Fig. 5. The simplest situation is the one obtained with ω 1 in both a) and b) sub-figures: only one (traveling) wavevector can afford this frequency. Since the φ i{η} functions (i = r, θ, z) do not depend on sign ±k p , using Eqs. (11,12,13) all the "weak clamp" conditions listed above reduce to U {η} (z = 0; L, t) = 0. It is then matched from Eq. (16) summing up two waves traveling in opposite directions ±k p , with opposite amplitudes U 0 and −U 0 , leading to: for which 2U 0 sin(ωt) factorizes out in Eqs. (11)(12)(13). Similarly form Eq. (17) one obtains: with the complementary phase for the time-oscillation. The wavevector k p should then verify: The same situation applies to the case ω 2 in Fig. 5 a) considering the full line, or to the lowest torsional and longitudinal branches (m = 0) which do not support evanescent solutions (see Appendix A). This is indeed trivially equivalent to what is usually done in textbooks when considering low frequency/long wavelength torsional/longitudinal modes in beam theory, see for instance Ref. [2]. Explicitly, the displacement field u {η} (r, θ, z, t) reads: for torsional (T) and flexural (F) solutions, and: in the longitudinal (L) case. Now on, we write U ω (t) = U 0 cos(ωt) + U π/2 sin(ωt) the generic time dependence; with the φ i{η} functions properly normalized, the amplitude of motion is encoded in the two quadratures of the motion U 0 , U π/2 (in meters). A given normal mode motion is completely defined by these amplitudes, together with the set {η} to which we add the standing wave parameter q. These expressions that we call "string-like" (in analogy with the high-stress solutions of Euler-Bernoulli theory) are the ones we will use below when defining mode energies.
Consider now the case ω 2 in Fig. 5 a) with the dashed line as evanescent contribution, or equivalently ω 2 in Fig.  5 b) with the full line branch (or the dashed one turned upwards). There are now two possibilities to support an excitation at ω 2 , one at finite k p and the other one at finite k e . Following the same philosophy as above, we should now combine Eqs. (16,17,19) with wavevectors ±k p , ±k e , factorizing out the same time-oscillation sin(ωt) [and similarly with Eqs. (16,17,18) for the π/2 shifted one cos(ωt)]; however now, the φ i{η} do not factorize out, since they are generally different for k p and k e waves. Besides, we need an extra boundary condition to fix the amplitudes of the evanescent components. This condition is chosen to be that the z-derivative of the amplitude parameter vanishes at the two extremities: We call this situation idealized "strong clamping" (as opposed to the previous one). Such a situation is never encountered by torsional (T) waves (see Appendix A). The approach matches the standard one performed in the framework of the Euler-Bernoulli modeling [2]. Besides, we shall discuss explicitly in Subsection III D how the boundary conditions are related to the energy flow toward the anchoring points.
After some manipulation, the solution writes for longitudinal (L) modes: with the generic time dependence U ω (t) and the difference between traveling and evanescent functions denoted by η p , η e respectively. For flexural (F) modes the situation at hand is again more complex. Eqs. (26,37) While it might seem to over-constrain the problem, the choice of the Θ F plane reference makes it guaranteed (see Appendix A). Postponing the discussion of this fact, we obtain: The constant writes: and the wavevectors k p , k e should satisfy: Considering k e (k p ) as a function implicitly parametrized through ω, Eq. (45) is the equivalent of Eq. (29) leading to a discrete set of k p possibilities that we index again with q > 0, and conveniently incorporate in {η}. This situation also applies to the first flexural branch n = 1, m = 0 for small enough frequencies, when both traveling and evanescent solutions exist (see Appendix A). It is in direct correspondence with the solving of the Euler-Bernoulli equation (at low in-built stress), which describes the same physics at small wavelengths/low frequencies (with the neat simplification k e = k p ) [2]. For this reason, We shall call "beam-like" the above mentioned expressions.
The situation encountered for ω 4 in Fig. 5 is rather trivial: only evanescent or mixed waves can sustain this frequency. One can show that in this case, it is impossible to match the boundary conditions; therefore for a given branch, no standing waves can exist as soon as there is no traveling solution available.
The situation depicted in Fig. 5 b) with ω 3 is the most interesting. In this case, two traveling waves match on the right panel. Consider first an evanescent part of the branch of the type depicted by the dashed curve turned upwards: there is then no evanescent solution. Strictly speaking, it is impossible to satisfy the boundary conditions at the same time for two k p having exactly the same ω: therefore, we shall have two distinct waves with very close ω, both verifying Eq. (29). If one considers the dashed evanescent branch turned downwards, there is then also one evanescent solution. In this case again we shall have two distinct solutions, with very close ω, but verifying this time Eq. (45). The most complex situation arises for the full line evanescent branch where two evanescent solutions match ω 3 . Four possible combinations will then coexist with very close frequencies, each related to a distinct equation of the type of Eq. (45). There is a sort of "degeneracy lifting" which we represented in Fig. 5 with a fuzzy and broader line for ω 3 than for the others.
Imposing the boundary conditions thus turned all the branches into infinite sets of standing wave modes, indexed through q > 0: these are the normal modes u {η} (r, θ, z, t) that describe all possible deformations of the rod in the limit of linear elasticity. Consider now an arbitrary standing wave u within the rod. It is obtained by a linear superposition of all the solutions we constructed, formally: This is our starting point for the mechanical description of the rod energetics.
The vast majority of the normal modes are of the "stringlike" nature, especially all asymptotic behaviors at large k p (which is of importance when dealing with convergence issues in the following Sections); only some branches support "beam-like" modes, at (reasonably) small k p . For this reason (and also pragmatically because all the analytics is much simpler), we will restrict in Sections IV and V the mathematical description to the one referring to Eqs. (30)(31)(32)(33)(34)(35). But all the presented concepts also apply to "beam-like" solutions, and we leave the derivation of the proper related expressions for a future work.

C. Stored energy: kinetic and potential
We can now define the mechanical energy stored in the rod by an arbitrary motion u. We thus start by recalling basics of continuum mechanics that can be found in textbooks, see e.g. Ref. [2]. The (volumic) density of kinetic energy writes: while the potential energy density E p is the sum of a bending term E E Y and a stretching term E σ defined by Eq. (A17), which is nonzero only when an axial stress σ 0 = 0 is inbuilt. The latter corresponds to the elastic energy stored in the rod by the length variations caused by motion, which we model by a body force. This approach is not exact but preserves the Pochhammer-Chree waves as being solutions. It neglects an addendum E 0 which is essentially irrelevant, see the discussion of Appendix A 4 for details. The bending energy density is defined in continuum mechanics from the strain [ kl ] and stress [σ ij ] tensors as: written in polar coordinates. The strain field is obtained from the displacement vector u: and the stress field from the elasticity relations, Eq. (9): with the bilinear forms E A {η},{η } (A = ρ, σ, E Y ) constructed from the expansions of the above expressions with respect to the u {η} of the sums. One can show that these expressions are symmetric under the exchange {η} ←→ {η } (Appendix C). The total energy density is simply: The bilinear forms introduced here are explicitly given in the following Sections IV and V for the simple case of "stringlike" modes.

D. Dissipation: internal and clamp friction
Energy can be fed into a mode (or can decay away) via two paths: through internal degrees of freedom or towards the two clamping ends. The internal contribution is chosen to reproduce a conventional viscous damping; in mechanics terms, this is usually realized by introducing a complex Young's modulus with an imaginary part representing the damping coefficient [31]. This simply means that locally, friction is modeled by a stress field [χ ij ] proportional to the rate of change of the strain [32]: for which we chose the same parametrization as for Eqs. (55-60) for simplicity. Two constants are thus introduced, the damping parameter Λ Y in Pa.s and ν λ that quantifies the directionality of the friction field with respect to the strain field (similarly to ν). The power lost through friction (per unit volume) is then obtained as: similarly to the stored bending energy Eq. (48). Note however the factor of 2 difference in definitions. The internal mechanisms leading to the definition of Λ Y are outside of the scope of this article. An introduction to high-temperature mechanisms can be found in Ref. [2]; conversely, internal two-level systems (defects for crystals, or structural for amorphous matter) are usually invoked as the dominant mechanism of low-temperature friction in NEMS and MEMS. For an illustration of this, see e.g. Ref. [33].
The modeling of the anchoring points (see Fig. 2) requires specific attention. First, one should point out that for −ε < z < 0 and L < z < L+ε, the writing of u based on Eqs. (11)(12)(13) cannot be valid anymore, since precisely in these regions of space the motion should gradually switch from rod-like waves to bulk-like waves. As such, we recast formally these equations into: where the "modal shape" now depends also on z. We assume that all the functions introduced are regular, so that we can write: for z around 0 (and a similar expansion for z ≈ L), and i = r, θ, z. Obviously ϕ i{η} (r, θ, 0) = ϕ (0) i{η} (r, θ) = φ i{η} (r, θ). Higher order functions n = 0 depend on the actual geometry of the clamp, and their description is outside of the scope of this manuscript. From these expressions, we can define the time-dependent stress field [σ c ij ] within the clamp, injecting the above in Eqs. (49-60). The power density flowing through the two regions of length ε is thus expressed similarly to Eq. (70): valid only for −ε < z < 0 and L < z < L + ε (and zero elsewhere by definition). The rate of change of the strain [∂ c kl /∂t] within the clamp is again defined from a linear response perspective. We therefore introduce two friction constants Λ c , ν c that fully characterize the ability of energy to radiate into the bulk supports. Formally, this is given by: keeping again the same type of parametrization as in the above. This approach generalizes the introduction of an admittance matrix when modeling only the lowest mechanical branches [35][36][37][38]; however, the calculation of this coefficient is outside of the scope of our manuscript. It suffice to point out that leakage of energy towards the anchor can be mininmised through a proper design of the clamp; this leads to the creation of "phononic shields" [39,40]. Note that Eq.
(75) contains all contributions that arise from the distortion, including the shear generated by the angle appearing between the anchoring point and the rod experiencing flexure [41,42]. Again, through a clever design distortions at the anchor can be minimised, which leads to a decrease in friction named "soft clamping" [32,43]. The final ingredient of the clamp modeling consists in the application of the boundary conditions already introduced in Subsection III B. Let us consider first the case of T and F solutions, for which we only impose U {η} (z = 0; L, t) = 0 yet. This leads to the expansions: in the vicinity of z ≈ 0 (and similar expressions for z ≈ L).
Using the replacement Eqs. (24)(25), the equivalent expansions for L modes write: theŨ in Eq. (85) referring to the primitive function with respect to the z variable (and evaluated at z = 0 here). In each set of equations, two functions which are constants-in-z (but functions of t) appear, which are defined from the mode shape function U {η} . One can argue that they represent two possible "loss channels" (for mode {η}), linked to the distortion pattern of the wave. Depending on the clamp condition (so-called "weak" or "strong"), which defines the type of mode confined in the rod ("string-like" or "beam-like" respectively), different "loss channels" are selected. This is intimately linked to the nature of the forces and torques applied onto the anchoring point by the rod's distortion. The concept of "channels" for energy flow is discussed further in the following Sections. Similarly to Eqs. (61,62), the friction terms can be expanded onto the normal modes as: with coefficients resembling those found for the bending expansion E E Y within Eq. (62). The total power density lost through friction is then: reminding the specificity that the power densityḊ clamp is nonzero only for −ε < z < 0 and L < z < L + ε (whilė D int is defined over 0 < z < L strictly). Note that for the modes to be well-defined, the total damping rate Γ {η} derived fromḊ tot should remain much smaller than the mode resonance frequency ω {η} (the high-Q limit, see following Subsection IV). Also, the friction coefficients Λ Y and Λ c which model specific damping mechanisms may appear to be frequency-dependent in the most generic case (as well as their respective ν λ and ν c ): we write them Λ Y (ω), Λ c (ω).
For simplicity, we shall drop the frequency dependence of the ν λ , ν c parameters (see Appendix B). Practically, within a high-Q approximation they can be treated as constant for each specific mode, and we simply define Λ Y ≈ Λ Y {η} (and similarly Λ c ≈ Λ c {η} ) when ω ≈ ω {η} . This is to be contrasted with the hypothesis made on the elastic constants E Y , ν which are assumed to be frequency-independent. As a consequence (and contrarily to the energetics coefficients discussed in Subsection III C), the difference between damping coefficients Λ Y {η} , Λ Y {η } and Λ c {η} , Λ c {η } leads in principle to an asymmetry in the bilinear forms introduced here under the exchange {η} ↔ {η }; these refined aspects are discussed in Appendix C. As for the energy expansions, the friction coefficients are detailed in Sections IV and V below for the simple case of "string-like" modes.

IV. IN-EQUILIBRIUM SITUATION
Let us first describe the simple case of thermal equilibrium. Each mode {η} experiences Brownian motion, such that the time-dependent motion U ω (t) in Eqs. (30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43) is not a simple sine-wave (at frequency ω) imposed by an external drive of some sort. It is from now on a stochastic time-dependent variable U ω {η} (t) which spectral weight lies mainly around ω {η} (Section VI). Consider the energy density expansions Eqs. (61,62). As will be shown in the following Section V, the sums {η} = {η } are the consequence of the presence of a thermal gradient; we shall thus take them to be zero for the discussion on which we focus here. For "string-like" solutions, injecting Eqs. (30-35) into Eqs. (47,A17,48), and integrating over the cross-section we obtain respectively: for torsional (T) and flexural (F) modes. For longitudinal (L) modes, one should swap cos ↔ sin. The writing of the introduced coefficients (which are k p -dependent) is given explicitly in Appendix B. We shall not discuss the much more involed situation encountered for "beam-like" solutions. We define: where the factor 1/2 appears because of the sinusoidal mode shape. These will enable to build the mode's effective dynamics equation in the following Section VI. Per construction these quantities verify: We write · · · the ensemble average on the stochastic variables. Summing up Eqs. (91-93), and using the property U 2 ω {η} = ω 2 {η} U 2 ω {η} justified by the spectral concentration of U ω {η} around ω {η} (see Section VI), we obtain the total average energy per unit length of the mode: where we made use of Eqs. (94,95), and ∆K {η} = ∆S {η} + ∆B {η} . We kept the two 1/2 in this writing to underline that one term refers to kinetic energy, and the other potential.
Postponing the discussion of the mode dynamics (Section VI), the equipartition theorem applied to mode {η} states that: with T {η} the temperature associated to its Brownian motion (and L the length of the rod). On a given branch, the "string-like" solutions verify k p = q π/L with q > 0 (and {η} containing both this index and the branch label, {m, T, q}, {m, L, q} or {m, n = 0, F, q} by definition; remember also that they are two degenerate flexural families θ 0 = 0 and +π/2, see Appendix A). Eq. (96) finally writes: Noticing that ∆M {η} /M {η} < 1 and ∆K {η} /K {η} < 1, the second line of this expression verifies: Summing up to an upper cutoff N cut i 1, we produce the total energy per unit length corresponding to this specific branch i (with i representing the branch label; now on we do not use the writing in {· · · } for it, in order to keep it clearly distinct from the mode labels): where T {η} is defined as the average temperature of the branch. For T branch i to be finite, all temperatures T {η} have to be bound by an upper limit. Therefore, from Eq. (99) we realize that the sum 1/N cut i Ncut i q=1 performed over the second line of Eq. (98) tends to zero for all 0 < z < L since 1/N cut i Ncut i q=1 cos 2π q z L → 0 for N cut i → +∞, excluding the two extremes z = 0; L. This is what is reminded by the 0 in the above expression.
We can now add together the energy per unit length of all the branches, leading to the total energy integrated over the cross-section as defined by Eq. (63). For simplicity, we will again neglect "beam-like" solutions and assume that all branches are described by "string-like" modes; this shall be enough to demonstrate our methodology. The calculated energy should match the macroscopic description's result 3N k B T avg /L of Eq. (8), here within the assumption T c = T h = T avg . It imposes that: 1 3N all branches i and wavevectors q i.e. all modes {η} with N the number of atoms in the solid rod. The first equation simply states that the total number of degrees of freedom is preserved, which is also what is performed in Debye's approach (see Section II). Physically, the existence of a finite number of branches and of (very large but) finite cutoffs N cut i is precisely due to the fact that a mode's shape cannot present distortions on a scale smaller than the inter-atomic distance (see Appendices A and B). The second equation means that the mode's mean temperature equals the macroscopic temperature of the two baths T avg .
We now compute the friction experienced by each mode, restricting again the discussion to "string-like" solutions. As previously, the {η} = {η } terms in our expressions shall be set to zero. After integration over the cross-section, the internal contribution Eq. (70) leads to: which is the power dissipated per unit length due to internal degrees of freedom. The ± sign stands for T and F modes (+), and L modes (−) respectively. As for the energetic coefficients, L {η} and ∆L {η} are given explicitly in Appendix B.
After integration over the length L of the rod, the cos(2k p z) expressions of Eq. (103) disappear. We can then identify an effective viscous force acting on mode {η}:  (82-87), we see that at lowest order this term is z-independent. Integrating over the clamp volumes 0 < z < ε and L < z < L + ε thus leads to: for ε R; higher order terms are outside of the scope of this article, and do not impact our methodology. The coefficient C {η} (clamp friction per unit length) is given in Appendix B for T, L and F modes. Depending on the type of solution considered (T, L or F, "string-like" or beam-like" with the corresponding boundary conditions), the expression of C {η} is different, selecting thus for each case a specific "channel" for the energy leakage towards the outside. This can be interpreted in terms of forces and torques applied onto the anchoring point, which depend strongly on the nature of the standing wave localized in the rod; an explicit discussion is presented in Appendix B. Since the dynamics of U ω {η} (t) is spectrally localized around ω {η} (Section VI), we can substitute ω {η} U ω {η} (t) ≈U ω {η} (t) within a term that does not carry energy on average. Eq. (105) therefore leads to the identification of an effective clamping viscous force acting on mode {η}: which physically decomposes symmetrically in two halves applied on each end of the rod; note that for a perfect clamping ε → 0, this term vanishes (see Appendix B and Section VI).
The two damping forces Eqs. (104,106) add up and fully characterize the total friction, Eq. (90), which will appear in the effective dynamics equation of mode {η} (Section VI). The fluctuation-dissipation theorem states that they are related to two stochastic Langevin forces, characteristic of each of the baths they refer to. Their correlation functions verify (δ 0 being Dirac's function): which simply means that the intrinsic response times of these baths are infinitely faster than any of the modes, and that they are uncorrelated. These forces are responsible for the Brownian motion of mode {η}; by definition the clamp is thermalised at a temperature T avg , while we define T int for the internal degrees of freedom. These two contributions can be conveniently recast in a single one, with: in which appears the total friction coefficient L {η} L + C {η} ε and the mode temperature written as: Therefore, when the internal degrees of freedom of the rod are thermalized to the same temperature as its two ends (which is the case if there is no extra source of heat or cooling acting upon them), we recover the simple case of T int = T avg . This leads to the in-equilibrium result: which obviously implies T branch i = T avg for all branches i, and verifies our former condition Eq. (102). The quality factor Q {η} associated to mode {η} is defined as: with Γ {η} the relaxation rate, and should verify Q {η} 1 in order to guarantee a well-defined resonance; which poses a condition to be fulfilled for the friction coefficients.

V. THERMAL GRADIENT SITUATION
Consider now that a weak thermal gradient has been established between the two ends of the rod. In order for the energy density to reproduce Eq. (8), the cross-terms in Eqs. (61,62) must be nonzero. Integrating over r − θ, we obtain: for the kinetic part, and: for the potential part; we define k p = q π/L and k p = q π/L with both integers q, q > 0 for each mode. The ± sign stands for T and F modes (+) and L modes (-). The modal coefficients introduced in this Section are listed in Appendix C.
Reminding that energy is localized around specific frequencies within the motion spectrum (see Section VI), we extend the in-equilibrium property mentioned in the previous Section: which we insert in Eq. (117) while taking the ensemble average of both contributions of the energy density. Besides, the scalar product < · · · • · · · > defined in Appendix A (not to be confused with the ensemble average just mentioned) tells us that modes from different branches are orthogonal, which means that they verify M {η},{η } = 0 (and ∆M {η},{η } as well). By construction, this implies that Eq. (117) is zero: such cross-terms cannot carry energy to create a thermal gradient, and shall disappear from the calculation. Only modes belonging to the same branch should be kept in our expansion, leading to the expressions in Eqs. (121,122) written below. The first term in Eq. (122) corresponds to our previous contribution Eq. (100), which is due to the mean temperatures T {η} ; the ± sign has the same meaning as above. We define: An educated guess already tells us that these correlations should quickly decay as |∆q| increases. The neighboring modes {η} , {η } in a pair with non-vanishing contribution to Eq. (122) should therefore be very similar; this actually means that the result Eq. (95) can be extended to the present case, with ω {η} ω {η } M {η},{η } /K {η},{η } ≈ 1 (see Appendix C for details). Furthermore, the sum on ∆q = q − q runs on negative and positive integers, and for most modes the actual boundary is irrelevant since correlations die away rapidly: only the very first modes on a branch (small q ≈ 1) and the very last (q ≈ N cut i ) will have a truncated, asymmetric sum on ∆q. This has no weight in the total expression Eq. (122) and can be safely neglected, which allows us to write: which is the generalization of Eq. (100), with T branch i,∆q = 1/N cut i Ncut i q=1 T {η},∆q (the average correlation function for branch i).
Summing-up all branches, and having imposed that the total number of degrees of freedom is preserved [as stated by Eq. (101)], we finally reproduce the thermal gradient shape Eq. (6) appearing in the total energy density expression Eq. (8) from the equality: where T avg,∆q = 1/3N all modes {η} T {η},∆q extends the previous in-equilibrium expression Eq. (102). The links between the different correlation parameters T {η},∆q , T branch i,∆q , and T avg,∆q shall be discussed in the following Section VI. Since the cos[ π L ∆q z] functions form a free family, this is easily solved producing: with ∆T the macroscopic temperature gradient between the two ends. All even terms are zero, and only odd ∆q allow correlations, which fall as ∝ 1/∆q 2 (such that only nearest neighboring modes are really relevant Similarly to the energy densities, the power friction densities Eqs. (88,89) also have {η} = {η } cross-terms. The internal friction contribution leads to: when integrated over the cross-section. The ± corresponds to T and F modes with +, and L modes with -(as above). The L {η},{η } , ∆L {η},{η } coefficients are given in Appendix C explicitly. Inspecting Eq. (126), we see that the internal friction density has a z-dependence similar to the potential energy Eq. (118). Again, integrating over the length of the rod these terms vanish from the effective dynamics equation (Section VI).
The situation encountered with the clamping terms {η} = {η } is in this respect slightly different. One obtains: with the C {η},{η } coefficients discussed in Appendix C. Eq. (127) means that for odd ∆q (the ones that carry correlations), these clamping cross-terms vanish too: again, the mode dynamics equation is free from such contributions (see Section VI). However, even though there is no net related viscous clamping force acting on the mode, there is a finite force acting on each end of the rod (with two opposite signs). This has to be contrasted with the in-equilibrium treatment which results in a finite clamping friction force (the two anchoring contributions have the same sign). Replacing Eq. (120) in Eq. (127), the total energy flowQ can thus be written as: considering only one end of the rod (in this writing, the sign refers to flow in the direction of z), and summing up all contributions. The findings of this Section shall finally be analyzed in the following one, in the framework of an extension of the fluctuation-dissipation theorem.

VI. EXTENDED FLUCTUATION-DISSIPATION THEOREM
The aim of this final Section is to interpret the previous findings in terms of effective stochastic forces δf {η} (t) acting upon each of the modes {η}. Putting together the ingredients derived in the two previous Sections and applying Newton's law, we derive the mode's effective dynamics equation: The standard procedure when considering a single mode is then to perform a Rotating Wave Approximation (RWA); we shall not do that here, and instead look for the exact solution. The RWA is discussed explicitly in Appendix E. Applying the same approach as the one of the standard Fluctuation-Dissipation Theorem, let us assume that the δf {η} forces are all centered ( δf {η} (t) = 0) and stationary. We define: the correlation function and corresponding spectrum of the force fluctuations. FT [· · · ] is the Fourier Transform (defined as +∞ −∞ · · · exp[−i ωτ ]dτ ). Similarly, we write: for the motion amplitudes of two arbitrary modes {η} , {η } (which might be equal or not). Eq. (129) leads to: We then have: in the high-Q limit, with Q {η} defined in Eq. (115). I Re and I Im are two functions of Q {η} /Q {η } described in Appendix D.
• When the frequencies ω {η} = ω {η } , the two resonance peaks are clearly separated: ω {η} − ω {η } is much larger than their respective damping rates Γ {η} , Γ {η } . This situation is encountered only on the lowest branches, at small wavevectors k p R, k p R 1. In this case one obtains:  The marginal situation ω {η} = ω {η } should however be commented, for practical reasons: this is the situation encountered with the modes that are easily at reach for the experimentalist [11][12][13]22]. In this case, the model leads to: with the normalization I Im (2.61) ≈ −1 chosen arbitrarily (Appendix D). This means that the stochastic force correlators verify C δf {η} ,δf {η } (τ ) = −C δf {η} ,δf {η } (−τ ): they are anti-symmetric by time-reversal. Such kinds of correlations have been introduced in the context of exotic superconductivity or magnetic orders [44], with a specific broken symmetry. If this is a plausible physical result for mechanical states that could lead to measurable correlations remains to be demonstrated.
To conclude the modeling, let us consider the energy floẇ Q as defined from Eq. (128). Inserting in this expression the mode parameters given in Appendix C, it can be recast into:Q with g {η},{η } a new set of adimentional parameters. Using again the argument that modes are very close, we can approximate g {η},{η } ≈ g {η} and rewrite the sum over q and ∆q following the same approach as in the previous Section V:Q Each term of the sum essentially acts as a conduction channel for the energy flow: the modes themselves are not truly transporting energy, but their correlations do. These correlations are nothing but the amplitudes of the mode's 2wave mixing, appearing because of the quadratic form of the energetic functionals (see Appendix E). Besides, surprisingly Eq. (146) does not depend specifically on the couplings C {η},{η } to the anchoring ends (which are our two baths), but only on the overall relaxation rates. There is therefore no reason to assume any specific difference between modes in their ability to carry the heat, and we will simply assume T {η},∆q = T branch i,∆q = T avg,∆q . The calculation simplifies then into:Q having introduced: 3N all branches i and wavevectors q i.e. all modes {η} which defines the average clamping loss parameterΛ c (in Pa.s), discussed in more details in Appendix C. The only requirement here is that it should be finite; injecting Eq. (125) into Eq. (152) we finally write: with n the atom volumic density. When modeling the clamps, we assumed that ε L; actually for the bracket to be finite, we realize that ε should even be second-order o(1/L 2 ). We therefore introduce a new lengthscale : which characterizes the anchoring zone. BothΛ c and are nonzero for a physical realistic system, and the ideal clamp limit ε → 0 appears as a trivial consequence of considering the rod infinitely long L → +∞. Comparing expression Eq.
(154) to Eq. (7), we identify the thermal conductivity κ: having introduced the Dulong-Petit specific heat C v = 3nk B , and the mean sound velocityc. The materialdependent intensive parameter in brackets corresponds to the microscopic parameterλ mf p , the phonon mean-free-path reminded in Section II. Eq. (156) therefore eventually links a microscopic lengthscale, the mean-free-path, to a macroscopic one (but nonetheless very small) which is the size that defines the clamping region . This is somehow a consequence of the diffusive limit of phonon transport: within this region, the nonlinear processes responsible forλ mf p generate the cross-correlations between modes; while the rod itself is treated within a completely linear theory.

VII. CONCLUSION
We present a theoretical treatment that links mechanics to thermal properties for a dielectric (thin-and-long) rod connected to two heat reservoirs. It is based initially on the basic observation that thermodynamics applied to mechanical modes must reproduce the outcomes of the microscopic phononic theory at the macroscopic level, in thermal equilibrium but also when a small heat flow is present.
The key hypotheses on which we built are: for mechanics, linear response theory up to the atomic inter-spacing, and a high enough temperature to ensure classicality to all the mechanical modes. For phononics, local equilibrium should be guaranteed within the system. From the Pochhammer-Chree waves that exhaustively describe the rod's mechanics [14], we reconstruct both the specific heat (mean temperature) and the thermal conductivity (from the T -gradient profile). The mathematics presented is limited to "string-like" modes, namely the ones constructed from pure traveling waves (and no evanescent contributions). As a result, the consequences of the model are captured within a modified version of the Fluctuation-Dissipation Theorem. It predicts cross-correlations between nearby mechanical modes, which should be measurable. These cross-correlations arise because what transports energy is not the modes themselves, but actually their mixing. This counter-intuitive outcome redefines what one calls a "propagation channel" for heat. Besides, the theory produces a fundamental link between a microscopic lengthscale, the phononic mean-free-path, and a macroscopic one which defines physically the clamping region (which is assumed to be small). The nonzero mixing amplitudes are generated there, through phonon-phonon interactions which are outside of the scope of the manuscript.
The theory should apply equally well to mechanical standing waves constructed from both traveling and evanescent ones ("beam-like" modes). Calculating the coefficients appearing in the model within this situation is of practical interest (but is particularly more difficult), since such modes are the ones which are the easiest to measure, like the first flexures. If such modes, which do not overlap in frequencyspace, can really present correlations is an open question that should be answered experimentally. Finally, the theory can also be extended to cantilevers and 2D systems (plates and membranes), which are extensively studied experimentally [13,22]. No matter what are the actual geometries and materials considered, the prediction of cross-correlations has to be robust; only the quantitative magnitude should depend on the details of each specific system considered. Eventually, applying these concepts to the quantum case would be extremely enlightening [45,46], making the bridge with the conventional mesoscopic physics treatment of conduction channels.
We look for propagative solutions to Eq. (10), and we will restrict our discussion to traveling and evanescent waves. The φ i{η} (r, θ) functions (i = r, θ, z) introduced in Eqs. (11)(12)(13) can be written [14]: with θ0 = 0 or +π/2, Jn Bessel's function of the first kind (of order n ∈ N) and J n its derivative. ke,p stands for either ke or kp, and the ± sign reads + for traveling and − for evanescent. We shall already point out that the wavevectors ke,p extend up to about ∝ 1/a where a is the inter-atomic distance; beyond this limit Eq. (10) is not relevant anymore. The parameters α {η} and β {η} verify: 2ν)] and ct = c0/ 2(1 + ν) are the longitudinal and transverse velocities respectively, with c0 = EY /ρ. A {η} , B {η} and C {η} are (dimensionless) constants which will be defined below from boundary conditions and normalization. Here we discuss only the situation where no body forces are present in the rod; the impact of in-built stress is explicitly addressed in Subsection A 4 below. For a given n and θ0, imposing the free-stress boundary condition on the periphery of the rod leads to families of solutions described in the following Subsections: T, L and F branches [14,[26][27][28]. From Eqs. (A4,A5), one finds the dispersion relations ω {η} (ke) and ω {η} (kp) which are depicted in this Appendix, which in turn define the α {η} , β {η} constants. Each can be indexed by m ∈ N.
On the vectorial space that contains the displacement field solutions (namely 3D-space vector u(r, θ, z, t) functions with r ∈ [0, R], θ ∈ [0, 2π] and {z, t} ∈ R 2 ) we define the scalar product (not to be confused with the ensemble average): written for solutions of the type described above. The explicit dependence to the z−derivative of U {η} functions comes from our writing using only real-valued solutions. When combining those in complex form (i.e. ∂/∂z → ±ke or ±ikp for evanescent and traveling respectively), U {η} U {η } can be factorized out. The key fact that should be pointed out is that all the branches are orthogonal with respect to this scalar product, namely: if η and η correspond to two distinct sets of {θ0, n, m}, even within the same family. Only within a branch do we have < u {η} • u {η } > = 0, obviously when η = η but also when considering two distinct wavevectors ke,p and k e,p of the same branch (T, L or F). This actually justifies why standing wave modes where created in Subsection III B within a single branch, and not mixing up the branches (and even the families); as a direct consequence, normal modes of different branches are also orthogonal. Implications of orthogonality are further discussed in Subsection V.
These solutions were first studied by Pochhammer and Chree [29,30]. An obvious question that comes to mind is: are there any other possibilities? In a more mathematical language, the solutions described being orthogonal they form a free family in the { u}-Hilbert space. The question is then: is this a complete set that forms a basis of this Hilbert space? This is indeed mandatory in order to write any solution u in the form of Eq. (46). When we constructed the solutions, imposing the z-derivative in Eq. (13) was a necessity to obtain propagating wave solutions, with both Eq. (14) and Eq. (15). The question then shifts to the construction of the φ i{η} functions. Counting these functions, we see that we have 2 × {n, m} solutions listed, with {n, m} ∈ N 2 . This means that the free family is (infinitely) countable, through a bijection with N; the cardinality of our set is ℵ0. Comparing with other problems involving second order differential equations with boundary conditions (e.g. Laplace's equation), which have the same cardinality (their solution space has dimension ℵ0), we can safely consider that our set is indeed complete and forms a basis; this is a simple (yet not fully rigorous) implication of the dimension theorem.
We now turn to explicit Pochhammer-Chree solutions. These functions are used to compute the energetics coefficients introduced in the core of the paper (and listed in their respective Appendices). They are produced using Mathematica; the presented numerical results are obtained with ν = 0.20.

Torsional (T) waves
The simplest solutions are torsional waves (T) obtained for θ0 = 0 and n = 0. Then, only φ θ{η} is nonzero, and reads: with the proper normalization φ θ{η} (r = R, θ) = 1 for any m: the parametrization of the motion is done through the tangential displacement at the periphery of the rod. Note that wavevectors (and ν) do not appear yet explicitly; the angular dependence to θ also vanished. The boundary conditions lead in this case to the simple equation: which defines β {m,T} for m ∈ N, and thus: with + for traveling and − for evanescent solutions. The first solutions are β {0,T} = 0 [then φ θ{0,T} (r, θ) = r/R], β {1,T} = 5.1356 · · · /R, β {2,T} = 8.4172 · · · /R. The lowest dispersion relations Eq. (A10) are plotted in Fig. 6. The first branch is strictly linear (velocity ct) and verifies ω {0,T} (kp = 0) = 0; it does not support evanescent solutions. This is the one discussed in beam theory, see e.g. Ref. [2]. The others have finite frequency at kp = 0; all evanescent solutions link their branches to ω = 0.
The asymptotic behavior of all branches is dω/dkp = ct at large kp: we then recover a linear dispersion law. Finally, note that none of the branches reconnect to one another, and that no local minima exist (apart at ke,p = 0); the branches never cross each other.

Longitudinal (L) waves
The next family is found by imposing θ0 = +π/2 and n = 0: these are longitudinal (L) waves. Then, only φ θ{η} is zero and motion exists both on the axis and radially. As for T waves, there is no θ dependence, but the solution is more complex.
The boundary conditions involve now two constants, A {η} and B {η} . They can be written in matrix form, and a solution exists only if the determinant is zero:     The last constant A {η} is defined through a proper normalization of the wave shape. Since the displacement along z on the axis essentially never vanishes (see final comment of this Subsection), and in order to match the common definition used for the first mode at small wavevectors, we chose φ z{m,L,ke,p} (r = 0, θ) = 1 for all modes m (and all ke,p). This means that applying the replacement procedure Eqs. (24,25), the amplitude of motion is actually encoded in U {m,L,ke,p} (z, t) → RdU {m,L,ke,p} (z, t)/dz. Example of displacement fields uz, ur are shown in Fig. 8.
The first branch m = 0 does not support evanescent solutions. It is linear at small wavevectors with velocity c l (and goes to ω = 0 at kp = 0), while for large wavevectors it is linear again with velocity cR, the Rayleigh velocity (the one of surface waves) [14]. A fairly good approximation (within typically 10 −4 for −0.45 < ν < +0.45) can be produced using: The kpR 1 situation is the one presented in the framework  of beam theory in textbooks [2]. All higher branches have finite frequency at kp, ke = 0, with a linear asymptotic dependence at large kp reaching the transverse velocity dω/dkp = ct. Fig. 7 shows a complex pattern of reconnections on the evanescent side, which strongly depends on ν [26]. On both kp and ke panels, one can see local minima dω/dkp = 0 at finite wavevectors in some of the wave dispersions (e.g. {1, L} right panel). These points connect to mixed wave solutions, which are outside of the scope of this paper [14]. Note that none of the branches cross each other.
A final comment on the numbering of the branches is in order. Actually, the value of ω {m,L} (kp = 0) for m > 0 does depend on the Poisson ratio ν, so there is no specific ordering per se: depending on ν, branches can exchange positions with respect to ω. On the other hand, there is a difference in nature between nearby branches which we will rely on for our labeling. We thus count branches while increasing frequency, ordering them with the extra condition: • even m (0, 2, 4, · · · ) are strongly longitudinal, • odd m (1, 3, 5, · · · ) are strongly radial.
The resulting labeling shown in Fig. 7 is not the one usually found in the literature, e.g. in Refs. [14,26]. Also as a consequence, our fixed parametrization of motion amplitude based on φ z{m,L,ke,p} (r = 0, θ) = 1 is perfectly adapted to even m branches, but not really to odd ones. It can also happen that for very specific ke,p values (on specific branches), the z component vanishes. But this happens only in a single point, and in all the vicinity around this point all calculations are well defined. Since in the end, the modeling presented in the paper is concerned only with intensive properties, the actual parametrization does not matter (and final parameters quoted in the core of the manuscript are always finite by construction), and we kept a fixed parametrization of motion amplitude (independent of m) for simplicity.

Flexural (F) waves
The situation encountered with flexural (F) waves is the most complex, with two indexes n > 0, m ≥ 0 and all three displacement components nonzero. The solutions produced by the π/2 rotation (i.e. θ0 = 0 or +π/2) are degenerate: they correspond to the same displacement field patterns (but rotated from each other), having the same frequencies.
The boundary conditions affect now 3 constants, which can be put in matrix form in a similar fashion to Eq. (A11). We write the resulting determinant as Eq. (A13) below, similarly to Ref. [14]. The same reasoning defines the dispersion laws ω {n,m,F} (ke) and ω {n,m,F} (kp), which then fix the α {n,m,F,ke,p} and β {n,m,F,ke,p} constants. From the matrix equation, we can define the A {η} and B {η} constants as a function of the last one C {η} . The lowest F branches are shown in Fig. 9 Again, the last constant (here, C {η} ) is defined from the normalization procedure. In order to match the intuitive result that suits the lowest F branch, the amplitude of motion has to be defined from the (planar) displacement uF = ur(r = R, θ, z, t) er + u θ (r = R, θ, z, t) e θ occurring on the periphery of the rod. This can be projected on the two x, y axes defining φF,x and φF,y (from φr and φ θ ) respectively; we show as an example these components for the {1, 0, F} wave in Fig. 10, calculated for the θ0 = 0 family. And indeed for small wavevectors, we see that the planar component of the motion is almost uniform in θ with φF,x ≈ 0 and φF,y constant: this is just the well-known flexure (along y) of the rod (with nonetheless a small φF,z component, see Fig. 10). Similarly with θ0 = +π/2 one obtains the same type of motion, but along x. This is what is studied using the simpler Euler-Bernoulli formalism [2]. a11 a12 a13 a21 a22 a23 a31 a32 a33 However for higher branches (or larger kp), the motion is much more complicated; in particular, there is a clear θ-dependent pattern (Fig. 10). Therefore, the normalization procedure must chose a given θ = ΘF angle. Again to keep it simple, we take ΘF = 0 for all n, m (and both evanescent or propagating solutions). The amplitude definition of U {η} is then fixed by constraining φ r{η} (R, ΘF ) 2 + φ θ{η} (R, ΘF ) 2 = 1 (or equivalently using the x, y notations). This is in particular what has been done for Fig. 10. Note that similarly to the L waves normalization, uF essentially never vanishes, which makes the procedure well defined for any branch (and almost any ke,p, see discussion of previous Subsection). As a result, φ θ{η} (r = R, θ) = +1 cos (n θ) and φ r{η} (r = R, θ) ∝ sin (n θ) for any {η}, which actually guarantees the boundary conditions introduced for the construction of standing waves, Subsection III B.
For each n subset, the F branches plotted in Fig. 9 display properties similar to the ones of the L branches shown in Fig. 8. We also clearly see on the left panel branches that directly reconnect to ω = 0, as in Fig. 7. However, noticeably the lowest branch {1, 0, F} does support an evanescent part. For small kp, ke 1/R, these dispersion relations are quadratic, which is what one expects from Euler-Bernoulli (with ω {1,0,F} [ke,p = 0] = 0) [2]. All {n, 0, F} (n > 0) tend asymptotically at large kp toward a linear dispersion law with dω/dkp = cR (the Rayleigh velocity). The other ones tend to the same dependence with ct velocity [14]. For a given n, the branches never cross.
Finally, we should mention that as for L waves the numbering chosen for the branches is again related to the specific nature of the motions, namely: • even m (0, 2, 4, · · · ) are strongly planar, • odd m (1, 3, 5, · · · ) are strongly longitudinal.
Such labeling issues are usually not discussed in the literature [14,27,28]. As in the previous Subsection, this means that the chosen normalization is not specifically adapted to odd m branches; but we use it for simplicity.

Effect of in-built stress
Consider that the rod is subject to an in-built axial stress σ0 (along z). We assume that the length L of the rod corresponds to the static rest shape under this load. An arbitrary displacement u(r, θ, z, t) causes then a stretching ∆L = δL, which stores an elastic energy (per unit volume) σ0 δL/δz. The stretched length L + ∆L = dL is defined from: which expansion (at lowest order) brings: The first term in Eq. (A15) leads to a linear term in the energetics, that causes a renormalization of the rest position for a given wave (a static term, see below); only the second one contributes to the total time-dependent potential energy with a quadratic term. We write the elastic energy as Eσ + E0: The truncation in the expansion Eq. (A15) therefore neglects nonlinear effects, which are outside of the scope of this article.
Eσ is used when defining the stored mechanical energy in the rod, see Subsection III C. The addendum E0 is discussed at the end of this Subsection; we shall show that this term can be reasonably neglected.
Eq. (A17) leads to a body force fσ that adds up to the righthand-side of Eq. (10). It reads: and in our definitions σ0 > 0 for tensile stress. The presence of this force modifies the dynamic equilibrium, but it turns out that the analytic expressions used to describe the Pochhammer-Chree waves without in-built stress are still solutions of this new Eq. (10). This is actually the reason behind the splitting of the stretching energy into a term which is analytically solvable, and the addendum E0 which should be analysed separately (see below). However, this comes at the cost of a modification of the dispersion laws Eqs. (A4,A5): introducing two correcting terms within the k 2 dependencies which are proportional to σ0/EY (a small parameter for realistic physical situations). ± signs again stand for traveling (kp, with +) and evanescent (ke, with −).
To our knowledge, the stress correction is not discussed in the literature devoted to waves in rods [14,27,28]. For simplicity, All the illustrating graphs in this Appendix have been realized with σ0 = 0. However, adding stress does not change significantly any of our conclusions, and it is therefore possible (but tedious) to take a finite σ0 = 0 in practical numerical evaluations. For all traveling waves except the first flexural one, the stored axial stress just slightly modifies the dispersion relation ω {η} . This is illustrated in Fig. 11 for the first traveling waves of each family, with a stored tensile stress σ0 = +0.1EY (beyond reach of conventional materials). With a compressive load, the graph is essentially mirror-imaged around the horizontal y = 1 line (frequencies decrease, instead of increase). At small wavevector kpR 1, the f i (norm.) 1.
2p correction is negligible for m = 0. At large wavevector kpR 1, the axial stress modifies the sound velocity (but does not affect the linear-in-kp dependence). On the evanescent panel, the connecting points of the branches are slightly displaced by the added stress; the dispersion curve is thus slightly above (compressive) or below (tensile) the unstressed solution.
The main impact of the stress is to turn the quadratic dispersion law for the lowest traveling flexural wave {1, 0, F} at small wavevector kp into a linear one, for σ0 > 0 (tensile). Reversely for a compressive load σ0 < 0, a cutoff kcutR ≈ 2 |σ0| /EY in kp appears and the smallest wavevectors cannot propagate anymore; this is shown in Fig. 12. For a rod of finite length L, identifying kcut = 2π/L produces the first buckling instability [34]. At large kpR 1, one obtains again a small renormalization of the sound velocity, as for the other waves (see Fig. 11). On the evanescent side, a similar graph can be produced while inverting the role of tensile and compressive loads. Therefore, there is no {1, 0, F} evanescent wave available at small wavevectors ke, for large σ0 > 0.
Let us now come back to the addendum term in the energetics E0 = σ0∂uz/∂z − 1 2 σ0 (∂uz/∂z) 2 . We first consider the quadratic term, which has been introduced in order to preserve the validity   of Pochhammer-Chree analytic solutions. We define the quantity: which quantifies the importance of the uz term with respect to the "normal" stretching contribution involving only ur and u θ . For torsional (T) waves, Rσ is exactly zero: stating E0 = 0 is exact. For all L and F waves, this expression tends to 0 for large kpR 1. However, for kpR 1 (and keR 1 as well) two situations have to be distinguished: • for L waves, Rσ tends to 0 when m is odd; it however grows as kpR 1 for m even, • for F waves, Rσ tends to 0 when m is even; it grows for m odd.
This is simply linked to the nature of the branch, as discussed in the previous Subsections: the correction is getting worse for waves which are strongly longitudinal, as could be expected. This is illustrated in Fig. 13 with the computed Rσ for traveling solutions {0, L} and {1, L}. The point is that for all m = 0 waves, the correction at kpR 1 is anyway negligible; the fact that neglecting E0 is not always rigorous becomes thus irrelevant. For the first flexure {1, 0, F} (which is m even), when kp,eR 1 the modeling discarding E0 is correct; and for all waves when kpR 1, neglecting the quadratic term in E0 is a good assumption. The main limitation of this approach is that the model essentially over-estimates the effect of stress around kpR ∼ 1 for all modes, and also for the first longitudinal mode {0, L} at small wavevectors [the correction in Fig. 11 b) should physically tend to zero when Rkp → 0, while it does not]. We shall not comment these facts any further, and will simply drop the E0 quadratic term.
We now come to the first term in the addendum E0, which is linear in displacement amplitude uz. Consider a modification of our solution: where we introduced a z-and-t independent contribution Φ {η} . As such, it does not impact the kinetic energy density Eq. (47) E offset depends only on the Φ r{η} , Φ θ{η} , Φ z{η} functions, while the other coefficients mix them with the φ r{η} , φ θ{η} , φ z{η} ones. We spare the reader the explicit expressions of these terms. Similarly, the expansion of the linear stretching term leads to: Comparing Eq. (A25) with Eq. (A26), we see that they can compensate each other if: for any {η}, provided the sums {η} ={η } vanish. This set of 3 equations defines the 3 introduced Φ {η} functions. The transformation Eqs. (A22-A24) therefore removes the linear term of E0 from the energy expansion, and replaces it with the energy density E offset .
For torsional (T) waves (or when we impose the situation σ0 = 0 with L and F waves), one simply has Φ {η} = 0. Physically, in all other situations the Φ {η} (r, θ) vector corresponds to a new definition of the rest position for wave {η}, which is not homogeneously zero anymore for all components; it has a specific shape within the r, θ coordinates in the rod. This new rest position should derive from a (static) body force entered in Eq. (10). Its energetic consequence is the static energy density E offset (r, θ). Note that the length L (and similarly the radius R) has been chosen to be the rest dimension under the axial load σ0. This means that the sum of all the static distortions should vanish, per definition: {η} Φ {η} = 0. The closure relations mentioned above for Eq. (A25) with {η} ={η } · · · = 0 should thus be a consequence of this property. One can define a static stress [σ 0 {η} ij ] associated to each mode by injecting the Φ {η} solutions into Eqs. (55-60). By definition one has {η} [σ 0 {η} ij ] = [σ0 ij ], which has a single nonzero component σ0 zz = σ0. The explicit calculation of these terms and the demonstration of these rules are outside of the scope of the manuscript. It suffice to say for our purpose that the linear term cannot couple different waves, and it does not affect the t-dependent energetics. In the core of the paper, the linear contribution of E0 (and the related Φ {η} functions) is thus neglected.

Appendix B: Modal coefficients, in-equilibrium
In this Appendix we present the coefficients entering the energetics expansions obtained for "string-like" solutions. The derivation of the corresponding "beam-like" ones is much more complex, and shall not be addressed here. We write: for the mass and spring constant per unit length coefficients, for torsional (T) and flexural (F) solutions. The corresponding "deviation terms" introduced in Section IV are: again written in the case of T and F modes. For Longitudinal (L) solutions, all the right-hand-sides should be divided by (kpR) 2 because of our choice of amplitude definition, see Appendix A.
The dimensionless coefficients appearing in the above expressions are explicitly listed in Tab. I. They are defined from the φ i{η} functions (i = r, θ, z), which means that in general they also depend on kp (which is included in the mode index {η}) and the Poisson ratio ν. Despite the similarities (see e.g. between some a i {η} and s i {η} ), we kept distinct notations for all terms in order to make the approach as clear as possible. From the (numerical) knowledge of the mode functions, all coefficients can be calculated (see thereafter).
Integrating the kinetic and potential energies over the length L of the rod, we define respectively the effective mass M {η} L/2 and the effective spring constant K {η} L/2 that describe the dynamics of Uω {η} (see Section VI). By construction, they verify ω 2 {η} = ( 1 2 K {η} L)/( 1 2 M {η} L). For "string-like" solutions, kp = (π/L) q with q > 0 an integer. On a given branch i, we see that this index runs up to about L/a, with a the lattice spacing. This defines roughly the cutoffs Ncut i, which have to fulfill our "Debye" assumption. For two nearby modes separated by ∆kp = (π/L) ∆q (with ∆q = 0 integer), we have ∆kpR 1: at large kp on all branches, or for small kp on high-order branches (which are essentially flat near kp = 0), the mode parameters are almost equal for kp and kp + ∆kp wavelengths. This is especially true for the frequencies ω {η} .
coefs.  The internal friction per unit length coefficient is given by: and similarly: for torsional (T) and flexural (F) solutions. As above, the longitudinal (L) case is obtained by dividing the right-hand-sides by (kpR) 2 . From the choice we made in the parametrization of the internal dissipation mechanism, the l i {η} are identical to the p i {η} given in Tab. I (i = 0, 1, 2, c), provided one applies the modification ν → ν λ . The last coefficients introduced in Section IV concern the clamping losses. Because of the specificity of the stress profile at the anchoring point for each mode, each situation is somehow different. We shall discuss here only the "string-like" solutions, which verify: for torsional (T) and flexural (F) modes, and: for longitudinal (L) ones [following the replacement procedure of Eqs. (85-87)]. q is the mode number related to the wavevector kp = q π/L that was introduced for "string-like" modes. Replacing these within the power density expressions, and integrating over the clamp volumes leads at first order in ε to: with: for T and F modes, and: for L ones, where we deliberately kept as prefactor a term (equal to 1) reminding that the two ends are included in the calculation. The c i {η} coefficients (i = 0, 1, 2, c) depend on both ν and νc; they are defined from the φ i{η} functions (i = r, θ, z), but contain also the first order expansion correction describing the clamping strain profile [functions ϕ (1) i{η} (r, θ) introduced in Eq. (74), Subsection III D]. They are listed in Tab. II. As for the other terms of Tab. I, they depend on kp through the index q (which is part of the mode label {η}). Note that the idealized clamping conditions chosen in Subsection III B for our description remain rather generic: the actual strain profile within the anchors is by construction taken into account by the ϕ (n) i{η} functions, n > 0. How to calculate them is another matter, outside of the scope of this manuscript. Another simplification made was to consider the damping parameters ν λ , νc independent of frequency. While there is no fundamental reason to impose this, it simplifies the writing of our coefficients in Tabs. I,II. Pragmatically, all frequency dependencies can be incorporated into the Λ Y {η} , Λ c {η} parameters without losing much generality. Technically, the perfect clamping is obtained from a geometrical argument ε → 0 (vanishing clamp zone), while the material-and-geometry-dependent Λ c {η} shall remain finite (characteristic of the way rod-waves radiate into the bulk support; a complex problem outside of our scope). The actual meaning of the clamp modeling through ε is further discussed in Section VI. coefs.
∂φ r{η} ∂r (r, θ) i{η} , the first order correction within the clamping zone). They depend on both ν and νc, by construction; for details, see text.
To conclude this appendix, we shall give some straightforward examples illustrating Tabs. I and II. The case of torsional (T) modes {m, T, kp} is rather enlightening, and fully analytic. The mass and spring constant coefficients reduce to: and: z{η} , the clamp dissipation constant writes: which is the only one required for T modes [see Eq. (B18)]. The internal damping terms l i {η} (i = 0, 1, 2, c) are obtained from Eqs. (B22-B25) with the replacement ν → ν λ .
The motion amplitude parameter Uω {η} , which is defined as the tangential displacement of the rod's periphery, is directly linked to the usual angular motion parameter Θω {η} by: We then easily derive the moment of inertia ρ Ω(m) π 2 R 4 and torsional rigidity E Y 2(1+ν) + σ0 Ω(m) π 2 R 4 , which match well-known definitions for the lowest m = 0 branch without stored stress σ0 = 0 (as they should) [2]. Note the presence of the extra term ∝ p 0 {η} in the potential energy for T modes m > 0. The clamp friction Eq. (105) can then be recast into C {η} εR 2Θ2 ω {η} , which leads to an effective friction torque τ z {η} acting onto each anchor: with the brackets containing all material-dependent parameters. Looking into Eq. (75), this torque is actually due solely to a shear σ c θz . This is the only clamp "loss channel" available to torsional waves within the first-order in ε modeling.
For flexural (F) modes, the problem at hand is much more complex and non-analytic. In Fig. 14 we show the coefficients a i {η} and p i {η} computed for modes {n = 1, m = 0, F, kp}, for small kp wavevectors and no axial stress applied. When Rkp 1, we derive a mode mass ρ πR 2 and a mode spring constant EY (1 + 1.5 ν 2 + 0.3 ν 3 + 15.0 ν 4 + 35.0 ν 5 ) π 4 R 4 k 4 p (both per unit length), which match usual Euler-Bernoulli definitions (in the limit ν → 0, as they should) [2]. The ν-expansion of the spring constant coefficient is fit on numerical data within 1 % in the range −0.35 < ν < +0. 35. The clamp damping term c 1 {η} (neglecting the ϕ (1) z{η} function) is plotted in Fig. 16 a). As for T modes, this is the only one relevant for a first-order in ε modeling. When Rkp 1, it behaves as ∝ (kpR) 4 , which is actually due to both σc rz and σ c θz components of the clamp's stress tensor. In Cartesian terms, these reduce to a stress tensor σc yz , and the clamp friction can be interpreted as an effective friction torque τ x {η} acting upon each anchor (for a flexure in the y direction, Fig. 10). This is the only clamp- ing "loss channel" available to these modes, in the limit of small ε.
In Fig. 15 we present the coefficients a i {η} and p i {η} computed for modes {m = 0, L, kp}, at small kp wavevectors (and no axial load). The solution is again found numerically; asymptotically at Rkp 1, the mass per unit length tends to ρ πR 2 (kpR) 2 and the spring constant per unit length EY (1 + 1.5 ν 2 + 1.0 ν 3 + 15.0 ν 4 + 30.0 ν 5 )(kpR) 2 πR 2 k 2 p . The (kpR) 2 term appearing in both expressions is due to our choice of motion amplitude, see Appendix A. As for F modes, the ν expansion has been fit on the calculated results, within 1 % in the range −0.35 < ν < +0. 35. Dividing mass and spring constants (per unit length) by (kpR) 2 , one recovers the usual beam theory expressions (in the limit ν → 0, as it should be) [2]. The clamp friction is more involved than in the previous cases, and requires 3 coefficients plotted in Fig.  16 b) [while assuming ϕ (1) i{η} = 0, i = r, θ]. It is dominated by c 2 {η} , which tends to ≈ π for kpR 1. Its origin lies actually in the components σc ii of the stress tensor (i = r, θ, z), with a clear predominance of the σc zz contribution for small kp. One can interpret it as an effective friction force Fz acting on each of the anchors, which represents the clamping "loss channel" available (at lowest order in ε) to longitudinal (L) modes. normalisation comes from our definition of motion amplitudes, see Appendix A. The stretching and internal damping terms can be obtained from these, following the rules described in the text. Plot realized with ν = 0.2 and no axial stress (σ0 = 0). In the limit Rkp 1, the terms a 1 {η} and p 2 {η} dominate, leading to the usual mass and spring constant expressions (see text).

Appendix C: Modal coefficients, thermal gradient
In this Appendix we extend the formalism of the previous one, in order to be able to deal with the situation where a thermal gradient is present. For "string-like" solutions (the much more complex "beam-like" case shall not be discussed), we define the cross-terms in the energetics {η} = {η } to be: for L ones. With modes of relevance being nearby, we simplify again the problem into g {η},{η } ≈ g {η} by replacing in the above that we repeat here for convenience. For our model to be meaningful, this sum should be finite. Experimentally, Λ c {η} is found to grow with increasing mode number (see e.g. Ref. [47] for a study on NEMS beams). We can therefore verify the convergence of the sum by checking the asymptotic behaviour of the coefficients appearing in Eqs. (C14,C15). For torsional (T) modes, one easily shows that g {η} → Cste as kp → +∞ (for any mode index m; for m = 0 it is strictly a constant). For F and L modes, the problem has again to be solved numerically (see Figs. 14, 15, and 16 for an example of parameter sets).

Appendix D: Integrated generalized spectra
The properties of the generalised stochastic Langevin forces acting onto the modes are derived from the integrals: which we will study here numerically. Two situations have to be considered: either the mode frequencies ω {η} , ω {η } are much closer than their relaxation rates Γ {η} , Γ {η } , or they are far apart.
Consider first the case ω {η} ≈ ω {η } , which is the most common one. Typically, this is a good approximation (within a few %) as By construction, we chose I Re (1) = 1 which reproduces the conventional case where the mode labels {η} , {η } refer to the same one; besides, I Im is bound and remains I Im < 0.5.
Because of the Q-dependence of the imaginary component of Eq. (D3), we see that this contribution will disappear in the high-Q limit. This is what justifies keeping only the real part of the force noise spectrum in Section VI, ensuring thus time-reversal invariance for the correlation functions. We fit I Re (x) = 2/(x 0.5 + 1/x 0.5 ) on the numerical data of Fig. 17 (dashed line).