In the last two decades, nano-mechanical resonators have risen as highly promising devices for mass sensing due to their ultrahigh sensitivity. They can be used to measure the mass and stiffness of single particles like small pollution particles, viruses, bacteria, or even proteins. These capabilities of the nano-mechanical resonators have allowed the birth of a new type of mass spectrometry with no need of fragmentation or ionization of the sample and therefore ideal to measure big masses, where conventional mass spectrometers have important problems. The shape and modes of vibration of a nano-mechanical resonator can be very different and the advantages and drawbacks of one respect to another is still a hot topic. A unified theoretical framework to describe the effect of particle adsorption on a mechanical resonator is fundamental but still lacks in the literature. In this work, we present such a framework and examine the particular case of a rod-like particle adsorbed on a cantilever beam for flexural and torsional vibrations as well as on a disk resonator for radial breathing vibrations and compare the results with finite element simulations.

## I. INTRODUCTION

Nano-mechanical resonators have attracted the attention of researchers for more than two decades due to their extraordinary sensitivity and their wide range of applications such as temperature sensors,^{1,2} fundamental quantum studies,^{3,4} biological sensors,^{5–7} or mass spectrometry^{8–12} to name a few. Their ability to change the mechanical properties as a response to an external perturbation makes them ideal devices as extremely accurate sensors. The improvement in fabrication techniques and the readout systems made possible to reduce the dimensions of these devices to nanoscale increasing drastically the sensitivity.^{13} One of the most popular applications of nano-mechanical systems is mass sensing. It is possible to measure the mass of small particles with extraordinary accuracy by detecting changes in the natural frequencies of the nano-mechanical device. This possibility opened the door to a new type of mass spectrometry called nano-mechanical mass spectrometry that does not need fragmentation and ionization of the sample and that covers a range of mass from hundreds of MDa to pg where conventional mass spectrometers cannot be used. In this range, we can find important microbiological entities like viruses, bacteria and fungi, as well as non-organic pollution particles like PM10, PM2.5, or PM0.1 that cause serious respiratory problems for people, especially, in big cities.^{14–17} The popularity of nano-mechanical spectrometry has increased rapidly in recent years and important achievements have been made from detecting individual proteins in real time^{9} to neutral nanoparticles^{12} or viruses.^{8} Apart from the mass, it has been proved that these systems can also measure the stiffness of individual particles.^{18} This is due to the local resistance that the adsorbed particle produces to the vibration of the resonator device. In a recent study, the mass and stiffness of *Escherichia coli* bacteria cells were measured.^{10} This is a truly remarkable achievement that will improve drastically the capability of these systems to distinguish two particles with similar masses. A lot of different shapes have been proposed for the resonator as a mass sensor, from cantilever beams (which are the most common) to nanowires,^{19} membranes,^{20} toroids,^{21} or disks.^{22} The mass sensitivity is inversely proportional to the mass of the resonator, and therefore, the smaller the resonator the more sensitive it will be.^{23} However, for nano-mechanical spectrometry, there is another important factor, which is the efficiency to capture the particles and therefore, the wider the area, the better the capture efficiency.^{24} In order to maximize both mass sensitivity and capture efficiency, the most logical resonator is a thin plate of any shape. In order to extract the mass or the stiffness of the adsorbed particle from the frequency shifts of the resonator, it is necessary to have a theoretical framework that relates these quantities. This relation is widely studied for flexural modes of cantilever beams;^{18,25} however, there is not much information about other devices and other modes of vibration and none considering the stiffness effect of the particle. In this work, we present a theoretical framework that describes the mass and the stiffness effect of an adsorbed particle on the eigenfrequencies of a resonating plate of any shape that is valid for any type of mode of vibration. We then examine the particular case of the effect of a rod-like particle adsorbed on the plate with different orientations on the flexural and torsional modes of cantilever beams as well as on the radial breathing modes of a disk resonator and compare the theory with finite element simulations (FEMs).

## II. GENERAL THEORY OF PARTICLE ADSORPTION ON A MECHANICAL RESONATOR

We study a uniform isotropic plate with volume $Vp=Sphp$, where $Sp$ is the area of the plate and $hp$ is thickness. We define a coordinate system where the *x* and *y* directions are in the plate plane (Fig. 1).

The origin of the z axis is set at the top surface of the plate. We will develop all the theoretical frameworks in the frequency domain so that all the functions that will appear later in the text are the Fourier transform of the original functions in the time domain. For the sake of simplicity, we will omit frequency dependency. Assuming the thickness of the plate much smaller than any of its planar dimensions, the stresses normal to the plate can be neglected and the stress–strain relations of plane stress approximation can be used.^{26} The displacement of every point of the plate is described by the vector $rp=[up,vp,wp]=A\psi p(x,y,z)$, where *A* is an arbitrary amplitude and $\psi p=[\psi x,\psi y,\psi z]$ is the eigenmode vector that satisfies the following normalization condition:

The strain energy of the plate per oscillation cycle at the $nth$eigenfrequency is then given by^{27}

where $Ep$ is Young's modulus,$\nu p$ is Poisson's ratio, and$\beta n$ is given by

where $\epsilon ijpn$are the strains inside the plate at the $nth$eigenfrequency. For convenience, in the above integral, the original coordinates $[x,y,z]$ have been transformed into the normalized coordinates$[X,Y,Z]$ such that $x=lpxX$, $y=lpyY$, and $z=hpZ$, where $lpx$ and $lpy$ are the characteristic lengths of the plate in the *x* and *y* directions such that $lpxlpy=Sp$. The integration domain $\Gamma p$ in Eq. (3) is the three-dimensional domain of the plate in the normalized coordinate system $[X,Y,Z]$. The kinetic energy at the $nth$eigenfrequency is given by

where $Mp$ is the total mass of the plate and $\omega n$ is the $nth$angular eigenfrequency. Exact solutions for the eigenfrequencies of the plate can be obtained by using the Rayleigh–Ritz method^{28} if the vibration mode shape is known. This method relies on the conservation energy principle and states that the mean kinetic energy per oscillation cycle is equal to the mean strain energy per oscillation cycle at the eigenfrequencies.

When a particle adsorbs on the plate surface, two different effects arise: (i) the total mass of the system increases; therefore, the frequency must decrease in order keep the kinetic energy constant and (ii) the particle in contact with the plate surface deforms along with the plate; thus, the corresponding increase in strain energy must be compensated by an increase in the frequency to keep the balance between the kinetic and strain energies. In a nutshell, the mass of the particle induces a frequency decrease and the stiffness of the particle induces a frequency increase. Assuming the size of the particle is much smaller than the size of the plate, the relative frequency change can be expressed as

where $Tan$ and $Uan$ are the kinetic and strain energies of the adsorbed particle at the $nth$ eigenfrequency, respectively. The displacement throughout the particle can be decomposed into rigid displacement due to the motion of the plate and the displacement due to the deformation of the adsorbate. However, assuming that the adsorbate is much smaller than the plate, the displacement due to the deformation of the adsorbate is much smaller than the rigid displacement coming from the movement of the plate, and therefore, it can be neglected in the calculation of the kinetic energy of the adsorbate,

where $Ma$ is the mass of the particle and $(x0,y0)$ are the coordinates of the adsorption position on the plate. If the adsorbate is much thicker than the plate, we must include the rotatory inertia of the adsorbate in Eq. (6) that is proportional to the square of the slope of the mode shape. However, in this work, we consider small adsorbates that are not much thicker than the plate, and therefore, we will neglect rotatory effects of the adsorbate. In general, the adsorbate can have any shape, and therefore, the plane stress approximation cannot be used to describe the deformation inside the adsorbate. Instead, the general constitutive relations of an isotropic and homogeneous elastic material must be used,^{27}

where $\epsilon klan$ and $\sigma ijan$ are the strains and stresses inside the adsorbate at the $nth$eigenfrequency; indices *i*, *j*, *k*, and *l* can be *x*, *y*, and *z*; and $Cijkl$ is a fourth order tensor given by

where $\delta ij$ is the Kronecker delta and Einstein's notation of repeated indices is being used.^{29} The strains and stresses inside the adsorbate are produced by strain transmission from the plate, and therefore, every component of the strain in the adsorbate can be expressed as a linear combination of the in-plane components of the strain of the plate at the point of adsorption,

where the indices *m* and *q* can be *x* and *y*, and the functions $aijmq$ satisfy the equilibrium equations,

The boundary conditions that must be satisfied are as follows: (i) the displacements, the in-plane strains, and the shear stresses $\sigma xz$ and $\sigma yz$ at the interface plate-adsorbate must be continuous, and (ii) the normal components of the stress at the free surface of the adsorbate must be zero.

Using Eq. (10), the strain energy in the adsorbate at the $nth$eigenfrequency can be calculated and expressed as

where the indices *r* and *s* can be as well *x* and *y*, and $\gamma mqrs$ is a fourth order tensor given by

where, for convenience, the coordinates inside the integral have been transformed such that $xyz=VaX\u2032Y\u2032Z\u2032$ where $Va$ is the volume of the adsorbate and $\Gamma a$ is the three-dimensional domain of the adsorbate in the normalized coordinate system $[X\u2032,Y\u2032,Z\u2032]$. The tensor $\gamma mqrs$ has the so-called major and minor symmetries of a fourth order tensor,^{30} i.e., $\gamma mqrs=\gamma qmrs=\gamma mqsr=\gamma rsmq$. Therefore, the entire tensor $\gamma mqrs$ can be described by six different coefficients $\gamma xxxx$, $\gamma yyyy$, $\gamma xyxy$, $\gamma xxyy$, $\gamma xxxy$, and $\gamma yyxy$. The coefficients $\gamma mqrs$ can be understood as indicators of how well a given component of the plate strain can be transmitted to the adsorbate and contribute to its strain energy. The strain transmission occurs through the area of contact between the adsorbate and the plate, and therefore, the coefficients $\gamma $ are extremely sensitive to the area of contact.^{18} If the contact area is large, more strain can be transmitted and the coefficients $\gamma $ will be larger. On the other hand, the coefficients $\gamma $ will be smaller if the contact area is smaller until the case when the contact between the adsorbate and the plate is just a point where no strain can be transmitted whatsoever and the coefficients $\gamma $ will be zero [Fig. 2(a)]. From the contact area, the strain propagates inside the adsorbate making a distribution dominated by the free boundaries and, therefore, depending on the shape of the adsorbate.

Substituting all the expressions calculated for the strain and kinetic energies in Eq. (5), we obtain a general expression for the relative frequency shift due to the adsorption of the particle,

One key feature of the coefficients $\gamma $ is that, for a given type of vibration, out-of-plane or in-plane, they are fully independent of the adsorption position and mode of vibration, and therefore, they can be obtained by multimode measuring techniques. However, it is important to note that, apart from the shape of the adsorbate, due to the boundary conditions at the interface between plate and adsorbate, the coefficients $\gamma $ might depend as well on the ratio between the thickness of the adsorbate and the thickness of the plate and on the ratio between Young's modulus of the adsorbate and that of the plate.^{25} This dependence is weaker as these ratios become smaller, i.e., $\gamma =\gamma 0+O(hahc,EaEc)$. It is also different for out-of-plane modes than for in-plane modes because the deformation in the first type comes from bending and in the latter type from stretching, and therefore, $\gamma $ will change depending on the type of modes that is considered. Another important aspect is that, for a non-axisymmetric particle, it is obvious that the stiffness parameters will depend on the angle of orientation of the particle because the coefficients $aijmq(x,y,z)$ will be different due to the different boundary conditions of the rotated particle. However, we can easily incorporate this effect to the theory because as a tensor, the parameters $\gamma mqrs$ transform into new parameters $\gamma mqrs\u2032$ in a coordinate system that is rotated an angle $\theta $ around the z axis [Fig. 2(b)] like a rotation of a fourth order tensor,

where the indices $\alpha ,\delta ,\lambda $, and $\mu $ can be *x* and *y*, and $\Theta \alpha \delta $ are the elements of the rotation matrix given by

It is interesting to note that the invariants of $\gamma mqrs$ will not depend on the angle of orientation of the particle. Therefore, these invariants can be very good candidates for a stiffness fingerprint of particles of any shape. There are two linear invariants that come from the two contractions $\gamma mmqq$ and $\gamma mqmq$ that can be expressed as

Apart from these two linear invariants, second order invariants can be also found by doing the four following contractions $B1=\gamma mqrs\gamma mqrs$, $B2=\gamma mmrs\gamma qqrs$, $B3=\gamma mmrs\gamma qrqs$, and $B4=\gamma mqqr\gamma mssr$. From Eq. (14), depending on the symmetry of the adsorbed particle, tensor $\gamma $ can be simplified. At this point, we would like to distinguish between three different particle geometries. The first one is a particle that has less than two orthogonal planes of symmetry perpendicular to the plate plane. For this type of irregular particles, tensor $\gamma $ cannot be simplified and, in general, the six different parameters must be taken into account independently. The second type of particles corresponds to the particles that have two orthogonal planes of symmetry perpendicular to the plate plane without being axisymmetric, like, for instance, rod particles. Assuming that for $\theta =0$, the planes of symmetry are aligned to the *x* and *y* axes, and then for $\theta =45$, the coefficients must satisfy that $\gamma xxxx\u2032=\gamma yyyy\u2032$ and $\gamma xxxy\u2032=\gamma yyxy\u2032$. From Eq. (14), it can be deduced that $\gamma xxxy=\gamma yyxy=0$, and therefore, for these kind of particles, we have only four stiffness parameters. The last type of particles corresponds to axisymmetric particles, where $\gamma xxxx=\gamma yyyy$ for any value of $\theta $, like spherical nanoparticles or nanodisks. Evidently, the stiffness effect for these particles is not affected by the angle of orientation, and then, from Eq. (14), it can be deduced that $\gamma xxyy=\gamma xxxx\u22122\gamma xyxy$, and therefore, we have only two independent stiffness parameters.

### A. Out-of-plane modes of vibration

For out-of-plane modes of vibration, the in-plane displacement components *u* and *v* can be neglected. The strains inside the plate take the following form:^{31}

Using Eq. (18) for the strains, the parameter$\beta n$ can be expressed as

where $\alpha n$ is given by

where $\Omega p$ is the surface of the plate in the normalized coordinate system $[X,Y]$. Using Eqs. (18) and (19), the relative frequency shift for the out-of-plane modes of a plate is given by

### B. B. in-plane modes of vibration

For in-plane modes of vibration, the displacement *w* can be neglected and the displacements *u* and *v* at first order do not depend on *z*. The strains inside the plate take the following form:

In this case, the parameter$\beta n$ can be expressed as

where $\alpha n$ is given by

And the relative frequency shift is given by

In Eq (27), it is important to notice that the parameters $\gamma mqrs$ are, in general, not equal to parameters $\gamma mqrs$ for out-of-plane vibrations because the functions $aijmq(x,y,z)$ are different for in-plane modes than for out-of-plane modes.

## III. RESULTS AND DISCUSSION

### A. Flexural modes of rectangular cantilever beams

Two very simple cases of out-of-plane modes of a plate are the flexural and torsional vibration modes of a cantilever beam with the length *L* much greater than width *b*. In this case, we choose $lpx=L$ and $lpy=b$. For a flexural mode of such a structure, $\u22022\psi zn(X,Y)\u2202Y2=\u2212\nu p(bL)2\u22022\psi zn(X,Y)\u2202X2$ and $\u22022\psi zn(X,Y)\u2202X\u2202Y=0$. The dependence on *Y* can be neglected, and the final expression for the relative frequency shift is given by

where $\gamma flex=\gamma xxxx+\nu p2\gamma yyyy\u22122\nu p\gamma xxyy$ and $\lambda n$ is the $nth$ root of the eigenvalue equation $cosh\u2061(\lambda n)cos(\lambda n)+1$ that corresponds to the eigenvalue associated to the $nth$flexural mode given by the following equation:

where $X=0$ at the clamped edge and $X=1$ at the free edge. Finite element simulations have been performed in order to validate Eq. (28). The aim of the simulations is to calculate the relative frequency shift of the first three flexural modes due to the adsorption of a small particle. The cantilever used in the simulation was a 60 × 8 × 0.1 *μ*m^{3} silicon cantilever, and the particle was a small parallelepiped of 400 × 70 × 70 nm^{3} with Young's modulus, density, and Poisson's ratio of 20 GPa, 400 kg/m^{3}, and 0.25, respectively. Two sets of simulations were done. In the first one, the relative frequency shift was calculated as a function of adsorption position $X0$ for three different angles of orientation, $\theta =0$ in which the long axis of the adsorbate is oriented along the longitudinal axis of the cantilever, $\theta =45$, and finally $\theta =90$. In the second set of simulations, the adsorbate was placed at $X0=0.405$, and the relative frequency shifts were calculated as a function of the angle of orientation. As expected, the stiffness effect is higher for $\theta =0$ because the main component of the strain for flexural modes of cantilevers is $\epsilon xxpn$. Figure 3 shows the comparison between the theory and simulations. In the first set of simulations, we used Eq. (28) with the analytical modes (29) (purple dotted line) and the agreement between the simulations (empty circles) and the theory is excellent except for the zone close to clamping. In this region, the cantilever starts to behave like a 2D plate, and in order to fit the simulations, the more general plate [Eq. (21)] must be used (black dashed line). In this case, we used the modes calculated by finite element simulations of the cantilever without any attached particle. It is remarkable how well the theory fits the relative frequency shifts calculated by finite element simulations for all positions and for all angles of orientation. For places where the curvature is higher, the stiffness effect is higher and the difference between different angles of orientation becomes larger. However, where the curvature is zero, only the mass effect is present and the relative frequency shift is negative and independent of the angle of orientation.

### B. Torsional modes of rectangular cantilever beams

For torsional modes, $\u22022\psi zn(X0,Y0)\u2202Y2\u226abL\u22022\psi zn(X0,Y0)\u2202X2\u226a\u22022\psi zn(X0,Y0)\u2202X\u2202Y$; therefore, the relative frequency shift can be expressed as

where, in this case, the eigenmode can be expressed as^{32}

where $X=0$ at the clamped edge, $X=1$ at the free edge, and $Y=0$ at the middle of the cantilever. The same sets of finite element simulations performed for flexural modes were also performed for torsional modes for $Y0=0$ with the same cantilever and the same adsorbate. Figure 4 shows the comparison between finite element simulations and theory, where now, the analytical modes (31) (purple dotted lines) clearly fail to accurately predict the relative frequency shift calculated by FEM. However, using Eq. (21), and the modes calculated by FEM (black dashed lines) again fit perfectly the results obtained by FEM. Notice that for torsional modes, the relative frequency shift at $\theta =0$ is the same as that for $\theta =90$, that is why the blue circles cannot be seen in Fig. 4. Also notice that the maximum stiffness effect occurs for $\theta =45$. Because we chose to place the adsorbate at $Y0=0$, there are no negative values of the relative frequency shift since the mass effect vanishes at $Y0=0$. Notice how the relative frequency shift for torsional modes is in general smaller than for flexural modes. This is due to the fact that $\gamma xyxy$ is smaller than $\gamma xxxx$, which means that, for this adsorbate, the strain $\epsilon xxpn$ can be transmitted to the adsorbate more efficiently than the shear strain $\epsilon xypn$. However, it is interesting that using torsional modes, we have direct access to the coefficient $\gamma xyxy$ that can be very valuable information for spectrometry applications.

### C. Radial breathing modes of disk resonators

A good example of in-plane modes is the case of radial breathing modes (RBMs) of disk resonators. In these modes, the disk expands and contracts radially, and for a homogeneous and isotropic material, the problem is axisymmetric. This means that all points of the disk that are equally distanced from the center of the disk are in the same mechanical state. This example is very illustrative to show how the above equations can be transformed to other coordinate systems depending on the symmetries of the problem. We use cylindrical coordinates $[\rho ,\varphi ,z]$ with the origin at the center of the disk to describe the movement of the disk. At first order, for RBMs, the displacements along $\theta $ and *z* directions can be neglected. The displacement along the radial coordinate at the $nth$ eigenfrequency is given by $rn=A\psi \rho n=un2+vn2$, where $\psi \rho n$ can be expressed as^{33}

where $P=\rho Rd$ being $Rd$ the radius of the disk, $Ji$ are $ith$order Bessel functions of the first kind, and $\lambda n$ is the eigenvalue associated to the $nth$ RBM of the disk that satisfies the eigenvalue equation $\lambda nJ0(\lambda n)+2\nu pJ1(\lambda n)\u2212\lambda nJ2(\lambda n)=0$. The corresponding in-plane strains are given by

The relative frequency shift is then given by

where $\alpha n$ is given by

Finite element simulations have been performed in order to validate Eq. (36). The disk used in the simulation was 10 *μ*m radius and 100 nm thickness GaAs isotropic disk with Young's modulus, density, and Poisson's ratio of 85.9 GPa, 5320 kg/m^{3}, and 0.31, respectively, and the particle was the same particle used for the simulations of the cantilevers. Again, two sets of simulations were done. In the first one [Figs. 5(a)–5(c)], the relative frequency shift was calculated for the first three RBMs as a function of radial adsorption position $P0$ for three different angles of orientation, $\theta =0$ in which the long axis of the adsorbate is oriented along the radial direction of the disk, $\theta =45$, and finally $\theta =90$. In the second set of simulations, the adsorbate was placed at $P0=0.3$, and the relative frequency shifts were calculated as a function of the angle of orientation. Figure 5 shows the comparison between the theory and the simulations. In the first set of simulations, we used Eq. (36) with the analytical modes (32) (purple dotted line) and the more general plate Eq. (27) for in-plane modes (black dashed line). In this case, the analytical solution fits perfectly the simulations for all the cases. A very interesting property of the radial breathing modes of disks is that we have access to $\gamma \rho \rho \rho \rho $, $\gamma \rho \rho \varphi \varphi $, and $\gamma \varphi \varphi \varphi \varphi $ and, therefore, to the first invariant $I1$. This is extremely useful for spectrometry applications because it will allow the identification of non-axisymmetric particles with much more accuracy because now the angle of orientation does not have to be taken into account. Interestingly, the strains $\epsilon \rho \rho pn$ and $\epsilon \varphi \varphi pn$ become equal as we approach the center of the disk, and therefore, the dependence of the relative frequency shift with the angle of orientation vanishes as we approach to the center. On the other hand, despite the ratio $EaVaEpVp$ is higher for these disks than for the cantilever used in other simulations, the relative frequency shift is smaller than for cantilevers. The reason for this is twofold, the strains produced in the disk are smaller for radial breathing modes, and the transmission of these strains to the adsorbate is also smaller for in-plane modes than for out-of-plane modes making the coefficients $\gamma $ also smaller. Therefore, with respect to the cantilevers, we can say that the radial breathing modes of disks resonators have some advantages in terms of spectrometry applications for non-axisymmetric particles, but the stiffness effect in the relative frequency shift is smaller, and therefore, we will need more precision to measure these changes. It would definitely be interesting to design a structure with out-of-plane modes of vibration to maximize the stiffness effect and, at the same time, that has access to either the invariant $I1$or $I2$ for spectrometry applications.

## IV. CONCLUSIONS

In conclusion, we have presented a theoretical framework that describes the effect of particle adsorption on the resonance frequencies of a vibrating plate of any shape. The frequency variation for out-of-plane or in-plane modes departs from the same fundamental equation involving on one hand the displacements of the plate at the adsorption point (mass effect) and, on the other hand, the in-plane strains at the surface of contact between the plate and the adsorbed particle (stiffness effect). The theory shows how the stiffness effect depends naturally on the orientation of the adsorbed particle and how the stiffness coefficients form a fourth order tensor with major and minor symmetries that allows the determination of some stiffness invariants that will not depend on particle orientation and therefore will be ideal candidates for particle identification applications. We have tested the theory with finite element simulations of the adsorption of a rod-like particle on a cantilever beam, for flexural and torsional modes, and also on a disk resonator for radial breathing modes with excellent agreement between simulations and theory. This theoretical framework, which so far was lacking in the literature, sets the basis for the particle adsorption theory in nano-mechanical spectrometry.

## ACKNOWLEDGMENTS

This work was supported by the European Union's Horizon 2020 Research and Innovation Program under Grant Agreement No. 731868—VIRUSCAN.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.