Anomalous scaling of flexural phonon damping in nanoresonators with confined fluid

Various one and two-dimensional (1D and 2D) nanomaterials and their combinations are emerging as next-generation sensors because of their unique opto-electro-mechanical properties accompanied by large surface-to-volume ratio and high quality factor. Though numerous studies have demonstrated an unparalleled sensitivity of these materials as resonant nanomechanical sensors under vacuum isolation, an assessment of their performance in the presence of an interacting medium like fluid environment is scarce. Here, we report the mechanical damping behavior of a 1D single-walled carbon nanotube (SWCNT) resonator operating in the fundamental flexural mode and interacting with a fluid environment, where the fluid is placed either inside or outside of the SWCNT. A scaling study of dissipation shows an anomalous behavior in case of interior fluid where the dissipation is found to be extremely low and scaling inversely with the fluid density. Analyzing the sources of dissipation reveals that (i) the phonon dissipation remains unaltered with fluid density and (ii) the anomalous dissipation scaling in the fluid interior case is solely a characteristic of the fluid response under confinement. Using linear response theory, we construct a fluid damping kernel which characterizes the hydrodynamic force response due to the resonant motion. The damping kernel-based analysis shows that the unexpected behavior stems from time dependence of the hydrodynamic response under nanoconfinement. Our systematic dissipation analysis helps us to infer the origin of the intrinsic dissipation. We also emphasize on the difference in dissipative response of the fluid under nanoconfinement when compared to a fluid exterior case. Our finding highlights a unique feature of confined fluid–structure interaction and evaluates its effect on the performance of high-frequency nanoresonators.

S.1. RESERVOIR DENSITY (ρ * b ) AND CONFINEMENT DENSITY (ρ * i ) Figure 1a shows the change in bulk density, ρ * b of argon as the pressure in the reservoir is increased between ∼ 10 − 10 4 bars. For the pressure range considered in this study, argon exists either in the vapor or the supercritical state 1 . For any higher pressure, argon turns solid and will prevent the filling of the SWCNT. A cubic spline fit is performed to the data points obtained from equilibrium MD simulations to relate the pressure and the density, ρ * b at any intermediate point in the range. Corresponding to the reservoir pressure (or density) range, the confinement density, ρ * i of argon in the filled SWCNT varies between ∼ 0.01−0.57 as shown in Figure 1b. It can be seen that ρ * i increases with ρ * b (or pressure) linearly with the linear slope changing at ρ * i ∼ 0.20. A bilinear fit is used to capture the relation. The fit can be used to obtain the bulk density corresponding to any confinement density or vice versa. This is later used to plot and compare the interior and exterior cases of dissipation in terms of the density. i of argon in the SWCNT versus bulk density, ρ * b (lower X-axis) and bulk pressure (upper X-axis) of argon in the reservoir. For argon, ρ = ρ * × 1687.5 kg/m 3 . The green squares with error bars (extremely small and not visible) in (a) and the red circles with error bars in (b) are the data points from equilibrium molecular dynamics (MD) simulations. In (a) and (b), the black line is a cubic spline fit and bilinear fit, respectively, to the equilibrium MD data. The light shades of blue, green and red in the background indicate the states of argon, which are vapor, supercritical, and solid, respectively.

S.2. COMPUTATION OF PHONON MODE SHAPES
The mode shapes of the doubly-clamped SWCNT resonator can be obtained by carrying out real space quasiharmonic (QHMR) analysis 2,3 of the equilibrated SWCNT structure. Using QHMR analysis, for a structure with N atoms, a mass-weighted 3N × 3N force constant matrix, Φ can be defined with its elements given by where V is the potential energy, j and k are the atom indices, α and β are the directions, where B is the number of basis atoms in the lattice structure of SWCNT. The large size of the force constant matrix makes the QHMR analysis computationally expensive.

S.3. REAL SPACE TO PHONON SPACE
The resonant motion of any phonon mode, k, of the structure involves the collective movement of all of its atoms. To understand the motion, a coordinate transformation using the phonon mode shapes, φ k is required, from the real space to the phonon space, which can be written as Here, s i,α can be displacement, velocity or force corresponding to the i th resonator atom along the α th direction in the real space coordinate system, where α takes values 1, 2, or 3 corresponding to the x, y, or z directions, and S k corresponds to the same quantity (displacement, velocity or force) for the k th mode in the modal space coordinate. Following this transformation, at any instant for some mode k, the displacement of the phonon mode (q k ), the velocity of the phonon mode (q k ) and the hydrodynamic force by the fluid on the phonon mode ((P k ) flu ) can be expressed as Here, u and v and (F ) flu is displacement, velocity, and hydrodynamic force, respectively, on the resonator atom i along the α th direction. k = 1, . . . , 3N corresponding to the 3N phonon modes and k = 1 is the fundamental flexural mode.

A. Amplitude
The phonon modes obtained from the real space quasiharmonic (QHMR) analysis 2,3 are usually orthonormalized, i.e., the inner product,φ kφl = δ kl , whereφ k andφ l are the k th and l th modal vectors (or modeshapes) and δ kl is the Kronecker delta function. The amplitude of any phonon mode obtained using the orthonormalized modal vector may not be its true measure of amplitude. Conventionally, the amplitude of oscillation of a resonator is measured at the point of maximum displacement in its volume. The oscillation amplitude of that point is regarded as the amplitude of the resonator. If the modeshape φ k is orthonormalized and its magnitude is maximum at the i th m atom along the α th m direction, then corresponding to any amplitude A o k of the displacement q k (t) of the orthonormalized mode φ k , its true measure of amplitude, A k is given by is the component of the modal vector, φ k evaluated at the i th m atom along the α th m direction.

B. Effective mass
A resonator being an extended object, one can choose to record the displacement at any point along its volume, and consequently, the amplitude of resonance is not defined uniquely. However, the system's potential energy is a unique quantity and should be independent of the measurement criteria. Thus, the mass parameter in the damped harmonic oscillator model is defined to be the effective mass of the resonator where the effective mass is an adjustment over the true mass of the resonator to ensure that the potential energy of the motion is uniquely defined for the system 6 .
Here, a derivation of the effective mass of a resonator is shown corresponding to the assumption that the amplitude of oscillation of the resonator is measured at the point of maximum displacement in its volume. The effective mass µ k corresponding to any phonon mode k can be used to calculate the maximum potential energy stored by the mode as (U p ) max = µ k ω 2 k A 2 k /2. µ k is a fraction (ξ k ) of total mass m of the resonator with N atoms as µ k = ξ k m. The maximum potential energy can also be calculated as sum of the maximum potential energies (U p i ) max stored by each atom (with mass m/N ) during vibration under the mode shape φ k . (U p i ) max can be expressed as ( where α = 1, 2 and 3 corresponds to the x, y and z directions, and A o k is substituted as During flexural vibration of long slender structures like a carbon nanotube, nonlinearity is prevalent. The nonlinear mode coupling effect on the fundamental flexural mode of an isolated SWCNT can be effectively represented by introducing a damping term and an amplitude dependent frequency term in the harmonic oscillator model of the flexural motion.
The presence of fluid additionally contributes to the damping. In this context, since we have restricted our study to smaller amplitudes, the damping and the frequency parameter can be treated as average quantities constant over the motion. Thus, a damped harmonic oscillator model is sufficient to describe the ringdown of the flexural mode as where q k represents the displacement of the flexural phonon. µ k , k k , and ν k are the effective mass, effective stiffness, and damping coefficient, respectively of the phonon mode. For an empty SWCNT in vacuum, µ k = ξ k m (refer Sec. S.4 B), where m is the mass of the resonator and ξ k is a fraction (< 1) which depends on the mode shape. Also, k k = µ k (ω * k ) 2 where ω * k is the resonant frequency of the phonon mode in a case where there is no damping, and is related to its resonant frequency ω k in the present case with damping as ω k = (ω * k ) 2 − ν 2 k 4 . The damping coefficient ν k and the inverse quality factor Q −1 k are related as ν k = ω k Q −1 k (refer Sec. S.4 E). The fundamental flexural phonon is represented by k = 1. In the case of empty SWCNT in vacuum, it is found that Q −1 1, which implies ω 1 ν 1 /2. Consequently, ω 1 ∼ ω * 1 and k 1 can be written as k 1 = ξ 1 m(ω 1 ) 2 . We represent m and ω 1 in the case of empty SWCNT in vacuum as m 0 and ω 0 1 , respectively. When the SWCNT is coupled with the interior or exterior argon, the effective mass of the flexural phonon can change due to the added mass by the fluid. However, fluid coupling doesn't alter the effective stiffness of the phonon mode. Considering this and accounting for the damping, the resonant frequency of the flexural phonon can be expressed as 7 where ∆m is the added mass due to the fluid during the fundamental flexural motion.

D. Damping time
The solution to the equation of motion (Eq. 1) is q k (t) = A k e −ν k t/2 cos(ω k t − θ). Here, A k is the initial amplitude, and A k e −ν k t/2 denotes the time decay of the initial amplitude.
The exponential decay can be characterized by a damping time τ d k as τ d k = 2/ν k .

E. Quality factor
The total energy of the oscillator at any instant is E k = µ kq 2 k /2 + µ k (ω * k ) 2 q 2 k /2, where µ k is the mass of the oscillator. Using the equation of motion, the rate of change of energy can be expressed as dE k /dt = −µ k ν kq 2 k . The energy lost during a period of oscillation will be Under weak damping (ν k 2ω * k ), the exponential factor in the integral can be approximated as unity (e −ν k t/2 ∼ 1) and consequently carrying out the integration would lead to ∆E k = . The maximum energy stored by the oscillator is at the initial instant as E = µ k (ω * k ) 2 A 2 k /2. The inverse quality factor, defined as the ratio of energy dissipated over an oscillation period to the maximum energy stored, can thus be

F. Relaxation time
The phonon relaxation time, (τ k ) ph denotes the decay time of the energy of the phonon modes. From the expression of rate of change of energy, dE k /dt = −µ k ν kq 2 k = e −ν k t (. . .), (τ k ) ph can be approximated as (τ k ) ph = 1/ν k = τ d k /2.

S.5. FLUID DAMPING CONSTANT FROM HYDRODYNAMIC FORCE
The hydrodynamic force, (P 1 ) th flu on the fundamental flexural mode (k = 1) can be calculated by projecting the hydrodynamic force, (F i,α ) flu on each resonator atom along the mode shape φ 1 (refer Sec. S.3). (P 1 ) th flu can be expressed, in the most general way, as the sum of a conservative force proportional to the displacement q 1 (t) of the flexural phonon and a damping force with memory, proportional to the velocityq 1 (t), such that 8 Here, . neq denotes average over different time instants during the nonequilibrium ringdown simulation. We note that q 1 (t) at any instant can be expressed as q 1 (t) =ȃ 1 (t)e −iω 1 t , wherȇ a 1 (t) = a 1 (t)e −iθ is the complex time-dependent amplitude of vibration and i = √ −1. For a small value of dissipation, a 1 (t) varies slowly over the periods of oscillation such that it can be approximated as a constant, i.e., a 1 (t) ∼ A 1 . Inserting q 1 (t) and its time derivative, q 1 (t) in Eq. 3, we can write Here,c f andc I f are the real and imaginary Fourier transform of c f (t), respectively, from time to frequency domain and evaluated at ω 1 .c f can be calculated from the slope of (P 1 ) th flu (t) neq with q 1 (t) neq . The fluid dissipation, D th flu can be computed using the hydrodynamic force from Eq. 4 and expressed in terms of c f as Using Eq. 5 and E sto = 1 2 ξ 1 m 0 (ω 0 1 ) 2 A 2 1 , the inverse quality factor can be calculated as

S.6. FLUID DAMPING KERNEL USING LINEAR RESPONSE THEORY
Linear response theory (LRT) 9,10 states that the response of the thermally equilibrated system to a small external perturbation can be predicted from the thermal fluctuations of the system at equilibrium. In case of the SWCNT resonator with frozen atoms, flexural motion is set by artificially moving the atoms along the mode shape φ 1 with a time dependence as A 1 sin(ω 1 t). During the nonequilibrium process, the Hamiltonian of the fluid in the coupled system can be written as H = H 0 + Aλ(t). Here, H 0 is the Hamiltonian before imparting any perturbation. λ(t) is the external perturbation and A is any observable conjugate of λ with respect to the Hamiltonian. In the present case, λ(t) ≡ A 1 sin(ω 1 t) and A is the hydrodynamic force by the fluid on the resonator atoms along mode φ 1 in response to the perturbation, i.e, A ≡ P 1 . Under this setup, LRT can be used to express the response of another physical observable B in the nonequilibrium state as 9 ∆B(t) neq = ∞ −∞ dt χ AB (t − t )λ(t ). Here ∆B(t) neq = B(t) neq − B(t) eq , and . neq and . eq are nonequilibrium and equilibrium ensemble averages, respectively. χ AB is the after-effect function given by χ AB (t) = −β d dt δA(0)δB(t) eq , where β = (k B T ) −1 , δA = A − A eq and δB = B − B eq are the thermal fluctuations in variables A and B respectively. Presently, we are interested in the response of P 1 (t) itself such that ∆B(t) neq ≡ ∆P 1 (t) neq = P 1 (t) neq − P 1 (t) eq and χ AB (t) ≡ −β d dt δP 1 (0)δP 1 (t) eq . Noting that ∆P 1 (t) neq can also be parameterized in terms of any c f (t) as ∆P 1 (t) neq = t −∞ dt c f (t − t )λ(t ), it can be shown that c f (t) = β δP 1 (0)δP 1 (t) eq .
Similar to Eq.