A Review on Theory and Modelling of Nanomechanical Sensors for Biological Applications

Over the last decades, nanomechanical sensors have received signiﬁcant attention from the scientiﬁc community, as they ﬁnd plenty of applications in many different research ﬁelds, ranging from fundamental physics to clinical diagnosis. Regarding biological applications, nanomechanical sensors have been used for characterizing biological entities, for detecting their presence, and for characterizing the forces and motion associated with fundamental biological processes, among many others. Thanks to the continuous advancement of micro- and nano-fabrication techniques, nanomechanical sensors have rapidly evolved towards more sensitive devices. At the same time, researchers have extensively worked on the development of theoretical models that enable one to access more, and more precise, information about the biological entities and/or biological processes of interest. This paper reviews the main theoretical models applied in this ﬁeld. We ﬁrst focus on the static mode, and then continue on to the dynamic one. Then, we center the attention on the theoretical models used when nanomechanical sensors are applied in liquids, the natural environment of biology. Theory is essential to properly unravel the nanomechanical sensors signals, as well as to optimize their designs. It provides access to the basic principles that govern nanomechanical sensors applications, along with their intrinsic capabilities, sensitivities, and fundamental limits of detection. produced by the adsorption of a typical E. coli bacterium on a silicon nitride cantilever and a doubly clamped beam of 200 × 20 × 0.6 µ m for the ﬁrst three ﬂexural modes. It can be seen that, because of the edge effects of the bacterium, the stiffness effect is considerably smaller when the bacterium is transversally oriented.


Introduction
Advances in nanotechnology have led to the development of plenty of novel functional devices that are revolutionizing science and their applications at all levels. Nanomechanical sensors belong to this group of novel and highly promising devices and are attracting great attention from the scientific community owing the variety of applications they may find [1,2]. Just to mention some examples, nanomechanical sensors have provided access to physical quantities such as single spin [3,4] and fundamental quantum state [5,6], detected the presence of chemical and biological species at extremely low concentration [7,8], characterized the mass and mechanical properties of biological entities with extraordinary precision [9,10], and characterized the forces [11] and motion associated with biological processes [12][13][14][15].
Regarding the working principle, nanomechanical sensors, which are based on the study of the deformation of flexible structures at the nanoscale, can be split into two main groups: static mode and dynamic mode sensors [16,17] (Figure 1a). While the static mode is based on detecting changes in the deformation state of the sensor, the dynamic one is based on detecting changes of the mechanical resonances of such flexible structures. Although the use of both operation modes dates from the mid 1990s, the dynamic mode spread more rapidly, and nowadays, it encompasses around 80% of the research works in this field (Figure 1b). The physical origin of the signals may come from a wide variety of sources depending on the operation mode and the particular application. Nanomechanical sensors were first applied for the atomic force microscopy (AFM) technique [18][19][20]. AFM serves to reconstruct the topography and physicochemical composition of surfaces at the nanoscale. In this case, the physical origins of the signals that enable characterizing the surfaces deal with the existing interaction forces between the tip of a micro-cantilever and the surface. Nanomechanical sensors comprise biological applications at every scale [21], going from monitoring vital signs to proteomics. At the macroscale, they have been applied for monitoring the corporal temperature [22], the heart rate [23], the respiratory rhythm, the blood pressure, and the level of sucrose in blood. Going down to the micro-and nanoscale, nanomechanical sensors can characterize the motion and forces associated with biomolecular processes and detect and identify the presence of biological entities at extremely low concentrations based on their mass and mechanical properties. In these cases, each operation mode is based on very different detection principles [24,25]. The static mode relies on inducing differential surface stress on the sensor [26][27][28], which is commonly a thin micro-cantilever. One surface of the sensor is functionalized with molecular receptors that exhibit high affinity to targeted molecules [14,29]. As molecular recognition takes place, differential surface stress leads to the deflection of the sensor. On the other hand, the dynamic mode relies on inducing changes on the sensor mechanical properties [30,31]. In most cases, these changes are associated with the addition of mass. In this operation mode, every surface of the sensor can be functionalized. When molecular recognition occurs, mass addition modifies the sensor resonance frequencies. Although cantilever beams are also a common choice for the sensor in this case, many other forms can be frequently found in the literature, from nanowires to membranes or more complex resonant structures [32][33][34][35]. One important aspect of nanomechanical sensors when applied to biological detection is that energy dissipation hampers their capabilities, such as their sensitivity and limit of detection. When immersed in fluids, dynamic nanomechanical sensors suffer a great amount of energy dissipation caused by the density and viscosity of the surrounding The physical origin of the signals may come from a wide variety of sources depending on the operation mode and the particular application. Nanomechanical sensors were first applied for the atomic force microscopy (AFM) technique [18][19][20]. AFM serves to reconstruct the topography and physicochemical composition of surfaces at the nanoscale. In this case, the physical origins of the signals that enable characterizing the surfaces deal with the existing interaction forces between the tip of a micro-cantilever and the surface. Nanomechanical sensors comprise biological applications at every scale [21], going from monitoring vital signs to proteomics. At the macroscale, they have been applied for monitoring the corporal temperature [22], the heart rate [23], the respiratory rhythm, the blood pressure, and the level of sucrose in blood. Going down to the micro-and nanoscale, nanomechanical sensors can characterize the motion and forces associated with biomolecular processes and detect and identify the presence of biological entities at extremely low concentrations based on their mass and mechanical properties. In these cases, each operation mode is based on very different detection principles [24,25]. The static mode relies on inducing differential surface stress on the sensor [26][27][28], which is commonly a thin micro-cantilever. One surface of the sensor is functionalized with molecular receptors that exhibit high affinity to targeted molecules [14,29]. As molecular recognition takes place, differential surface stress leads to the deflection of the sensor. On the other hand, the dynamic mode relies on inducing changes on the sensor mechanical properties [30,31]. In most cases, these changes are associated with the addition of mass. In this operation mode, every surface of the sensor can be functionalized. When molecular recognition occurs, mass addition modifies the sensor resonance frequencies. Although cantilever beams are also a common choice for the sensor in this case, many other forms can be frequently found in the literature, from nanowires to membranes or more complex resonant structures [32][33][34][35]. One important aspect of nanomechanical sensors when applied to biological detection is that energy dissipation hampers their capabilities, such as their sensitivity and limit of detection. When immersed in fluids, dynamic nanomechanical sensors suffer a great amount of energy dissipation caused by the density and viscosity of the surrounding fluid. Thus, for liquids, the natural environment of biology, reducing energy dissipation constitutes a very important task. In order to circumvent this problem, researchers have suggested the application of many novel devices and techniques, such as the use of channeled [36][37][38][39] or partially immersed resonators [40], monitoring higher order modes [41], modes with less dissipation [42], or implementing active feedback loops [43], among many others. To date, the ones providing better results are the channeled resonators [44]. These devices comprise internal channels through which the liquid flows. Such configuration significantly mitigates the energy dissipation caused by the presence of the fluid, keeping the sensors' capabilities unaltered. Unfortunately, it is very difficult to miniaturize channeled devices because their associated microfluidic resistance rapidly increases, and because of limitations in the fabrication process. More recently, researchers have suggested the application of optomechanical platforms [45][46][47][48][49][50][51] as an encouraging alternative. Optomechanics provide access to mechanical modes that possess extraordinary properties, in terms of sensitivity and speed, even when immersed in the most dissipative fluids.
Since their origin, experimentalists and theorists have worked closely, aiming at the implementation of more sensitive, rapid, effective, robust, and integrable devices. In this spirit, advances in micro and nanofabrication techniques have enabled the production of continuously smaller devices, made of plenty of different materials and geometries. Together with the rapid evolution of the devices, researchers have made a huge effort in developing more accurate theoretical models, which enable accessing the basic principles that govern their applications and allow the correct interpretation of the sensors' signals, providing access to more, and more precise, information about the particular entities or processes of interest. Moreover, theory is essential to assist the design of the devices, as it serves to predict their capabilities regarding a given application, as well as their intrinsic sensitivities and their fundamental limits of detection [52]. This article focuses on the evolution of theoretical modeling in the field of nanomechanical sensing for biological applications, particularly gathering the most popular theoretical models currently applied.
Theories and models for the static operation mode date as far back as 1909, when G.G. Stoney published his formula to relate the surface stress with the bending curvature of a free plate [53]. Although this formula was not intended for biological applications, it later became the basis for the surface stress calculation due to a biological layer adsorption on a plate [14,[54][55][56]. However, this formula did not account for the clamping effect of cantilevers. Much later, Sader et al. proposed a solution for a cantilever plate in the asymptotic limits of high and low aspect ratio of the plate [57]. Some years later, Tamayo et al. obtained very simple solutions for the transversally averaged curvatures of the plate that can be used to calculate the surface stress from the experimental measurements of the cantilever deflection in a very easy way [58]. Regarding the dynamic mode, the first time that the use of the change in frequency of a resonator for bio-sensing applications was published was in 1997 in a work presented by Thundat et al. [59]. Three years later, in 2000, Ilic et al. demonstrated the detection of Escherichia coli bacteria using an array of low-stress silicon nitride resonant cantilevers [60], where they used simple formulas to relate the change in frequency to the mass of the analyte. This change was inversely proportional to the mass of the sensor and directly proportional to the mass of the analyte, and it was calculated for a homogeneous layer of material adsorbed on the cantilever surface. For small particles, the adsorption position plays a crucial role and must be known in order to calculate the mass of the analyte correctly. In 2005, Dohn et al. proposed a method to calculate the mass and position of the adsorbed particle using different modes of vibration [41] and, years later, in 2012, Hanay et al. proposed a method for the same task based on probability density functions [9]. In 2006, Ramos et al. proposed a model where, not only the mass, but also the stiffness effect of the analyte on the resonant frequency of the cantilever was accounted for [30]. This model did not account for the edge effects of the analyte, which were finally included by Ruz et al. in 2014 [61] and in a more general framework in 2020 [62]. These and other theoretical models that will be revised and explained later in the text can be considered as the basic theory for the field of nanomechanical sensing for biological applications. As we have emphasized, nanomechanical sensors can provide very valuable information about a very broad variety of physical properties of the biological entities or processes under study. Theoretical models are essential to properly process this information, as well as to design more sensitive and effective sensors. In this review, we aim to summarize the main theoretical models applied in this field. We first address the static operation mode and then we move to show basic principles and models of the dynamic mode. We then focus on particular models applied when nanomechanical sensors operate in liquids. In each section, the theoretical models increase in complexity, in the same way as they evolve along the history.

Nanomechanical Sensors in Static Mode: Effect of Surface Stress
The static bending of a thin plate can be caused by an applied force perpendicular to the plate plane or by a bending moment. The latter can be produced by different means. For instance, a plate formed by two layers of different materials will suffer bending whenever it is exposed to a temperature change because of the different thermal expansion coefficients of the two layers [63,64]. Larger temperature changes will produce larger bending on the plate, thus the static bending of the plate can be directly related to the temperature change, making it a highly sensitive temperature sensor [13,65,66]. A similar phenomenon occurs when one face of the plate is covered by a layer of a biological material [25,67,68]. The interaction forces between molecules within the biological layer produce surface stress on the face of the plate [69][70][71]. If the other face of the plate is uncovered, the difference of surface stress between the upper and lower faces of the plate will produce a bending moment that is proportional to the difference of surface stress between the two faces, or in this case, to the surface stress generated by the biological layer. In order to extract the surface stress from the measurement of the bending of the plate, the relation between the bending and the surface stress must be known. The differential equation that governs the bending of a rectangular plate is the so-called biharmonic equation, which, in Cartesian coordinates, is given by [72] where w(x, y) is the out-of-plane displacement of the plate, which is set to be along the z coordinate direction. Equation (1) assumes that there are no net external forces applied to the plate. Let us consider a rectangular plate of length L, width b, and thickness h ( Figure 2a). The most common cases of boundary conditions are free and clamped. When there is a bending moment M applied to the plate, the free boundaries satisfy the following conditions: where n and s refer to the components normal and parallel to the boundary, respectively, and E c and ν are the Young's modulus and Poisson ratio of the plate, respectively. The first to obtain the solution of Equation (1) together with boundary conditions (2) and (3) and relate it to the difference of surface stress between the upper and lower faces of the plate was George Gerald Stoney in 1909 [53]. The solution is a quadratic function of x and y and can be expressed as Equation (4) states that the curvature produced by the differential surface stress is constant throughout the whole plate. The bending moment is related to the differential surface stress Δ as = Δ .
Therefore, the curvature of the plate κ can be expressed as Equation (5) is the well-known Stoney's equation and has been used to quantify the differential surface stress on plates since it was first published in 1909. However, for biological applications, it is common to use a plate in a cantilever configuration, i.e., one of the edges is clamped so that its movement is completely restricted. Assuming the clamped edge at = 0, the boundary conditions for a clamped edge are Solution (4) is not valid anymore for a cantilever plate and Stoney's equation loses accuracy to describe the curvature, especially near the clamping region ( Figure 2b). Because of the lesser influence of the clamped edge, it is obvious that Stoney's equation will describe better the curvature for high aspect ratio cantilevers. The analytical solution of Equation (1) with cantilever boundary conditions was attempted by Zeng et al. [73]. They proposed a solution in the form of an infinite series of hyperbolic and trigonometric functions. However, the calculation of the coefficients of the series requires solving an ill-conditioned infinite system of equations with slow convergence, which suggests the clamped boundary condition and the free boundary conditions are in conflict at the common corners. John Elie Sader [57] found approximated solutions to the problem in the asymptotic limits of high and very low aspect ratios. For ≫ , Sader's solution takes the following form: where the functions ( ) and ( ) are given by Equation (4) states that the curvature produced by the differential surface stress is constant throughout the whole plate.
The bending moment is related to the differential surface stress ∆σ as M = h 2 ∆σ. Therefore, the curvature of the plate κ can be expressed as Equation (5) is the well-known Stoney's equation and has been used to quantify the differential surface stress on plates since it was first published in 1909. However, for biological applications, it is common to use a plate in a cantilever configuration, i.e., one of the edges is clamped so that its movement is completely restricted. Assuming the clamped edge at x = 0, the boundary conditions for a clamped edge are Solution (4) is not valid anymore for a cantilever plate and Stoney's equation loses accuracy to describe the curvature, especially near the clamping region ( Figure 2b). Because of the lesser influence of the clamped edge, it is obvious that Stoney's equation will describe better the curvature for high aspect ratio cantilevers. The analytical solution of Equation (1) with cantilever boundary conditions was attempted by Zeng et al. [73]. They proposed a solution in the form of an infinite series of hyperbolic and trigonometric functions. However, the calculation of the coefficients of the series requires solving an ill-conditioned infinite system of equations with slow convergence, which suggests the clamped boundary condition and the free boundary conditions are in conflict at the common corners. John Elie Sader [57] found approximated solutions to the problem in the asymptotic limits of high and very low aspect ratios. For L b, Sader's solution takes the following form: where the functions g 1 (x) and g 2 (x) are given by where the coefficients τ i depend on the Poisson's ratio and are given by Equation (8) states that the curvature of the cantilever follows Stoney's equation far from the clamped edge and then decays exponentially with a characteristic length of the order of b as we approach the clamped edge.
A different approach was published by Tamayo et al. [58] for cantilevers with aspect ratio L b > 1. They obtained solutions for the averaged curvatures along the transversal coordinate y. In this manner, the clamped boundary conditions are relaxed and the conflict at the corners disappears. They obtained very simple formulas for the averaged transversal and longitudinal curvatures: where the coefficient α depends on Poisson's ratio and is given by Equations (12) and (13) are extremely simple and can be used to calculate κ with very high precision from the experimental measurement of the cantilever deflection.

Nanomechanical Sensors in Dynamic Mode
Nanomechanical sensors have different modes of vibration, each one associated to its corresponding resonance frequency [74]. Although the resonator is a mechanical structure made of a continuum material and its movement should be described by the theory of elasticity, for certain applications, it is possible to model it as a harmonic oscillator with a certain effective mass (m e f f n ) and spring constant (k n ). In this manner, the resonance frequencies can be defined as Let us define an external perturbation, which can be of very different nature, from temperature, pressure, or humidity changes to force gradients, adsorption of individual molecules, or surface stress produced by a complete layer. This perturbation can change the elastic constant or the effective mass of the resonator, causing a shift in the resonance frequencies. These changes can be measured and carry valuable information about the external perturbation itself. This basic principle has allowed the nanomechanical sensors in dynamic mode to become an extremely useful tool for detecting biological entities and measuring its mechanical properties. In this section, we will review the theoretical models for the effect of the adsorption of analytes such as individual particles, as well as complete layers, on the resonance frequencies of the nanomechanical resonator, the two main effects used for biological applications.

Effect of a Complete Layer
There are two different effects that arise from the adsorption of a layer of biological material on the resonator surface. The first one is the surface stress produced by the Processes 2021, 9, 164 7 of 25 interactions between the molecules that formed the adsorbed layer. The second is the result of the added mass and stiffness due to this new layer of material.

Effect of Surface Stress
As we have mentioned before, surface stress generated by a biological layer on one of the faces of a cantilever plate can be obtained by measuring the static deflection [75]. However, surface stress also alters the dynamics of the plate, even if both faces are covered with the biological material and there is no bending at all. When a string is under tension, it vibrates at higher frequencies because the net force increases the effective stiffness of the string. Similarly, for a doubly clamped beam with a net surface stress applied σ T , the frequency of the beam increases proportionally to the surface stress and the relative frequency shift for the fundamental mode can be expressed as [76,77] ∆ f where f 0 and f are the resonance frequencies before and after deposition, respectively. However, cantilever plates can release the stress through the free edge, thus it is reasonable to think that the frequency will not change with surface stress. Whether or not surface stress affects the frequency of cantilever plates has generated a long debate among scientists of the field for the last decades [78][79][80][81]. Recent experimental studies have clearly shown that surface stress influences the stiffness of cantilever plates, changing the resonant frequency [76,82], and several theoretical models have been proposed to explain this phenomenon. At this point, it is important to distinguish between net surface stress, which is the sum of the surface stress on the upper and lower faces of the cantilever plate, and differential surface stress, which is the difference between the surface stresses of both faces. A net surface stress produces in-plane deformation of the cantilever, while a differential surface stress produces bending. A theoretical model to describe the effect of net surface stress on cantilevers was proposed by Lachut et al. [77,83]. In these works, they showed that there are two different effects, one coming from the geometrical changes of the structure due to the stress release and another one from the fact that the clamped edge cannot release the stress along the transversal direction, thus a residual stress is accumulated near the clamping region, affecting the resonant frequency of the cantilever plate. They proposed a formula for the relative frequency shift of the fundamental mode of cantilevers: As an example, Figure 3a shows the relative frequency shift due to net surface stress for a silicon nitride cantilever with a thickness of 100 nm and length of 20 µm, for different values of the aspect ratio, and as a function of the net surface stress. Equation (17) is formed by two terms that correspond to the two different effects described above. The geometric term is inversely proportional to the thickness of the cantilever, while the residual stress term is proportional to b 3 Lh 3 . It is interesting to note that these two effects have opposite signs; therefore, a cantilever could be designed with specific characteristics in order to avoid net surface stress effects on the frequency or, on the other hand, to maximize the effect in order to have a better measurement of the surface stress.
More recently, a theoretical framework to describe the effect of differential surface stress on the stiffness of cantilever plates has been developed [84,85]. In these works, it was found that two-dimensional bending of thin plates leads to the development of in-plane stresses that highly increase the flexural stiffness of the plate. In addition, a nonlinear effect due to large static bending was also included in the theoretical model. The equation of the effect of differential surface stress on the relative frequency shift of cantilever plates for the first three flexural modes of vibration is given by where c = 0.7745 is a constant, and the coefficient α n accounts for the nonlinear effect due to large static bending and takes the values of ≈ 0.02029, −0.11553, and − 0.04670 for the first, second, and third flexural modes, respectively. The function Λ n (ν) can be approximated by a second-order polynomial of Poisson's ratio for the first three flexural modes: The two-dimensional bending effect is proportional to b 4 h 6 ; therefore, it will be the dominant effect for very thin plates. On the other hand, as the aspect ratio of the cantilever increases and the nonlinear term becomes more important until, finally, in the asymptotic limit L b 1, the two-dimensional bending term is negligible with respect to the nonlinear term. Figure  3b-d show the relative frequency shift due to differential surface stress for a silicon nitride cantilever with a thickness of 100 nm, length of 20 µm, and different values of the aspect ratio as a function of the differential surface stress for the first three flexural modes.  The theoretical prediction of the relative frequency shift due to the differential surface stress Δ for the same cantilever and for the first, second, and third flexural modes, respectively, as well as for different values of the width. Insets show a representation of the mode shape.
More recently, a theoretical framework to describe the effect of differential surface stress on the stiffness of cantilever plates has been developed [84,85]. In these works, it was found that two-dimensional bending of thin plates leads to the development of inplane stresses that highly increase the flexural stiffness of the plate. In addition, a nonlinear effect due to large static bending was also included in the theoretical model. The equation of the effect of differential surface stress on the relative frequency shift of cantilever plates for the first three flexural modes of vibration is given by where = 0.7745 is a constant, and the coefficient accounts for the nonlinear effect due to large static bending and takes the values of ≈ 0.02029, −0.11553, and − 0.04670 The theoretical prediction of the relative frequency shift due to the differential surface stress ∆σ for the same cantilever and for the first, second, and third flexural modes, respectively, as well as for different values of the width. Insets show a representation of the mode shape.

Mass and Stiffness Effects
When a homogeneous layer of a material is deposited on one of the surfaces of the sensor, the total mass of the system increases, thus the resonance frequency must decrease in order to maintain the energy balance [86]. However, the new layer deforms along with the resonator, increasing the strain energy of the system, thus the resonance frequency must increase in order to keep the mechanical energy constant. Depending on the properties of the new layer, it might even change the strain distribution inside the resonator, considerably moving the neutral axis towards the adsorbed layer. In a nutshell, the mass effect produces a decrease in the resonance frequency and the stiffness effect produces an increase in the resonance frequency. For a general layer, the relative change in frequency for any flexural mode of a cantilever sensor can be expressed as where ρ is the density; η is the ratio between the thickness of the layer and the thickness of the cantilever; and subscripts a and c refer to analyte and cantilever, respectively. It must be noticed that Equation (20) is valid for both cantilever plates and doubly clamped beams. One important feature of Equation (20) is that the relative change in frequency does not depend on the mode of vibration; therefore, it is the same for any flexural mode. This is a particular characteristic of the adsorption of homogeneous layers that cover the whole surface of the resonator. For a thin layer of a biological material where η 1 and E a E c , Equation (20) can be simplified and the relative frequency shift can be approximated by where we have considered that the change in frequency is small and, therefore, Figure 4 shows the relative frequency shift as a function of η due to the adsorption of an homogeneous layer of a biological material with density of 800 kg/m 3 on a silicon cantilever for different values of E a . sensor, the total mass of the system increases, thus the resonance frequency must decrease in order to maintain the energy balance [86]. However, the new layer deforms along with the resonator, increasing the strain energy of the system, thus the resonance frequency must increase in order to keep the mechanical energy constant. Depending on the properties of the new layer, it might even change the strain distribution inside the resonator, considerably moving the neutral axis towards the adsorbed layer. In a nutshell, the mass effect produces a decrease in the resonance frequency and the stiffness effect produces an increase in the resonance frequency. For a general layer, the relative change in frequency for any flexural mode of a cantilever sensor can be expressed as where is the density; is the ratio between the thickness of the layer and the thickness of the cantilever; and subscripts and refer to analyte and cantilever, respectively. It must be noticed that Equation (20) is valid for both cantilever plates and doubly clamped beams. One important feature of Equation (20) is that the relative change in frequency does not depend on the mode of vibration; therefore, it is the same for any flexural mode. This is a particular characteristic of the adsorption of homogeneous layers that cover the whole surface of the resonator. For a thin layer of a biological material where ≪ 1 and ≪ , Equation (20) can be simplified and the relative frequency shift can be approximated by where we have considered that the change in frequency is small and, therefore, ≅ . Figure 4 shows the relative frequency shift as a function of due to the adsorption of an homogeneous layer of a biological material with density of 800 kg/m on a silicon cantilever for different values of .

Individual Particles
One of the most important applications of nanomechanical resonators is particle identification. Because of its extremely high sensitivity, nanomechanical resonators have

Individual Particles
One of the most important applications of nanomechanical resonators is particle identification. Because of its extremely high sensitivity, nanomechanical resonators have allowed the birth of a new type of mass spectrometry known as nanomechanical mass spectrometry [9,10,32,[87][88][89][90]. This type of mass spectrometry does not require ionization and fragmentation [90] of the sample, and thus is ideal for biological entities like viruses [32] or bacteria [10]. Like complete layers, individual particles that adsorb on the resonator surface have two opposite effects on the resonance frequency; that is, the mass effect, which decreases the frequency, and the stiffness effect, which increases the frequency. However, unlike complete layers, the frequency shift is in general different for different modes of vibration. The mass effect is associated with the kinetic energy of the adsorbed particle, thus it will be proportional to the square of the mode shape at the adsorption position [24,91]. On the other hand, the stiffness effect is associated with the strain energy of the adsorbed particle as a consequence of the bending of the resonator, thus it will be proportional to the square of the curvature of the mode shape at the adsorption point. The first formula that included the stiffness effect in the relative frequency shift due to the adsorption of a small particle type of analyte was proposed by Ramos et al. [30] and is given by where V is the volume; X = x/L is the longitudinal position normalized to the length of the resonator; X 0 is the normalized adsorption position; and ψ n . is the nth mode of vibration, given by where λ n is the eigenvalue that satisfies the following equation: where the + sign corresponds to cantilevers and the − sign to doubly clamped beams. In Equation (22), the mode of vibration has been normalized so that However, Equation (22) does not account for the relaxation of the stress through the free edges of the adsorbed particle. This relaxation is very important for small particles with a small surface to volume ratio, thus the stiffness part of Equation (22) is not accurate for this kind of particle. In 2014 [61], the edge effects of the particle were carefully studied and introduced into the stiffness part of Equation (22). These effects are responsible for the dependence of the stiffness effect on the shape, the contact to the cantilever surface, or the angle of orientation of the adsorbed particle with respect to the cantilever main axis. This proposed model includes a stiffness coefficient γ f lex that measures the strain transmission from the cantilever to the adsorbed particle, leading to the following equation for the relative frequency shift: γ f lex can be expressed as a function of the angle of orientation of the particle with respect to the main axis of the cantilever θ, just by applying a rotation along the vertical axis: where γ 1 , γ 2 , and γ 3 are coefficients that depend on the shape of the particle, the radius of contact, and η. Figure 5 shows the relative frequency shift due to the adsorption of a typical E. coli bacterium on a cantilever and a doubly clamped beam made out of silicon nitride as a function of the adsorption position X 0 for the first three flexural modes. As the E. coli bacterium has a cylindrical form, the relative frequency shifts for both longitudinal and transversal orientation were represented. nitride as a function of the adsorption position for the first three flexural modes. As the E. coli bacterium has a cylindrical form, the relative frequency shifts for both longitudinal and transversal orientation were represented. Recently, a more general framework to describe the stiffness effect of an individual particle adsorbed on a resonant plate of any shape has been developed [62]. In this work, it was found that the coefficient is actually a fourth-order tensor with major and minor symmetries that can be described by six independent coefficients. The framework describes the relative frequency shift due to the adsorption of a particle on a plate of any shape vibrating at any type of mode, in-or out-of-plane: where is a number related to the eigenvalue of the mode; ( , ) is the strain component of the plate at the adsorption position ( , ); the indexes , , , and can take the values and ; and Einstein's summation of repeated indexes is used. Note that the mode shape is in bold, indicating that it is a vector that, in a general case, can have , , and components. For flexural modes of cantilevers or doubly clamped beams, Equation (28) reduces to Equation (26). It is interesting to note that the tensor has two linear invariants to rotations around the vertical axis: These two quantities are very useful for particle identification applications because, as invariants, they do not depend on the angle of orientation of the particle, a variable that cannot be controlled. Recently, a more general framework to describe the stiffness effect of an individual particle adsorbed on a resonant plate of any shape has been developed [62]. In this work, it was found that the coefficient γ is actually a fourth-order tensor with major and minor symmetries that can be described by six independent coefficients. The framework describes the relative frequency shift due to the adsorption of a particle on a plate of any shape vibrating at any type of mode, in-or out-of-plane: where β n is a number related to the eigenvalue of the mode; ε mq (X 0 , Y 0 ) is the strain component mq of the plate at the adsorption position (X 0 , Y 0 ); the indexes m, q, r, and s can take the values x and y; and Einstein's summation of repeated indexes is used. Note that the mode shape is in bold, indicating that it is a vector that, in a general case, can have x, y, and z components. For flexural modes of cantilevers or doubly clamped beams, Equation (28) reduces to Equation (26). It is interesting to note that the tensor γ mqrs . has two linear invariants to rotations around the vertical axis: These two quantities are very useful for particle identification applications because, as invariants, they do not depend on the angle of orientation of the particle, a variable that cannot be controlled.

Multimode Measuring
The fact that the relative frequency shift is different for different modes of vibration when an individual analyte is adsorbed on the resonator surface allows for the determination of the adsorption position, mass, and stiffness of the particle just by measuring the relative frequency shift of several modes simultaneously [9,10,32]. This is commonly known as solving the inverse problem. Mathematically, Equation (26) using N modes constitutes a system of N equations with three unknowns, which are ∆ m = m a 2m c , ∆ s = 3E a V a 2E c V c γ f lex , and X 0 . In principle, three modes could be enough to solve the problem because we have three unknowns. However, the dependency on the adsorption position X 0 is nonlinear, thus the system might have more than one solution. Very different scenarios can occur depending on the adsorption position. For instance, for cantilever beams, ∆ m can never be determined for X 0 = 0 and ∆ s can never be determined for X 0 = 1. Another critical for any value of n. In this case, we can determine ∆ m + ∆ s , but never both of them separately. In practice, for real-time applications, the relative frequency shifts can be measured with a certain level of precision that is ultimately determined by the frequency stability of the system [92]. Assuming that the frequency noise is Gaussian, we can express the probability density function (PDF) of the relative frequency shift for the nth mode as follows [9,10]: f is the relative frequency shift measured experimentally and σ Allan is the Allan deviation of the corresponding frequency for the integration time used in the experiment. For N modes measured simultaneously, we must use the joint probability density function (JPDF) of the PDFs of all the modes. This function is a multi-normal distribution given by where∆ f f 0 is a vector formed by all the relative frequency shift variables,∆ f is a vector formed by all the relative frequency shifts measured experimentally, and Σ is the covariance matrix of the N modes. Obviously, Equation (32) has a maximum for∆ f f 0 =∆ f , and this information is not very useful. However, if∆ f f 0 is substituted by Equation (26), the new function will depend on the three unknowns of the problem X 0 , ∆ m , ∆ s . Then, maximization of Equation (32) for the new variables will lead to the solution of the problem. If the number of modes is not enough to obtain a unique solution, the JPDF will display several maxima for different combinations of the variables X 0 , ∆ m , ∆ s . The precision of the solution will be directly related to the width of the peak of the JPDF. Notice that, as a PDF, Equation (32) can be integrated in any of the three variables in order to eliminate the dependence on that variable. For instance, if we are interested in the mass probability density function, we can integrate in X 0 and ∆ s as follows: Similarly, different probability density functions can be obtained by integrating in the variables in which we are not interested. Figure 6 shows a schematic of the multimode method for the adsorption of a particle on a cantilever beam. Figure 6. Schematic of the particle's mass and stiffness determination using multimode measuring. In a first step, the frequencies of the resonator are tracked and the relative frequency shifts produced by the particle adsorption are recorded. In a second step, the joint probability density function ( ) is formed using the frequency noise for each of the modes being tracked. Finally, the can be integrated in the normalized position and either in Δ or Δ to obtain the mass or the stiffness , respectively.

Inertial Imaging
Multimode measuring can also be used to calculate the mass distribution of an analyte that is not considered as a point mass particle. Hanay et al. proposed a formalism based on the calculation of different moments of the mass distribution of the adsorbed particle [87]. The method neglects the stiffness effect of the adsorbed particle and is based on the fact that any polynomial function ( ) can be approximated by a linear combination of the square of the mode shapes . The relative frequency shift of the nth mode due to the adsorption of a particle that is not considered as a point can be expressed as where ( ) is the longitudinal mass distribution of the particle along the beam axis.
Equation (34)  , and so on. It is clear that ( ) corresponds to the total mass of the adsorbed particle, ( ) / ( ) is the center of mass of the particle or , and ( ) / ( ) is related to the size of the adsorbed particle projected to the axis. Higher moments can be calculated in the same way and give information about the shape of the particle. In order to calculate these moments from the relative frequency shifts, the functions ( ) ( ) are expressed as linear combinations of the mode shapes as Figure 6. Schematic of the particle's mass and stiffness determination using multimode measuring. In a first step, the frequencies of the resonator are tracked and the relative frequency shifts produced by the particle adsorption are recorded. In a second step, the joint probability density function (JPDF) is formed using the frequency noise for each of the modes being tracked. Finally, the JPDF can be integrated in the normalized position X 0 and either in ∆ s or ∆ m to obtain the mass or the stiffness PDF, respectively.

Inertial Imaging
Multimode measuring can also be used to calculate the mass distribution of an analyte that is not considered as a point mass particle. Hanay et al. proposed a formalism based on the calculation of different moments of the mass distribution of the adsorbed particle [87]. The method neglects the stiffness effect of the adsorbed particle and is based on the fact that any polynomial function p(x) can be approximated by a linear combination of the square of the mode shapes ψ n . The relative frequency shift of the nth mode due to the adsorption of a particle that is not considered as a point can be expressed as where µ a (X) is the longitudinal mass distribution of the particle along the beam axis. Equation (34) is reduced to the case of the point-like particle choosing µ a (X) = m a δ(X − X 0 ), where δ(X − X 0 ) is the Dirac delta function. The different moments m (k) of the mass distribution of the adsorbed particle are defined as m (k) = 1 0 µ a (X)g (k) (X)dX (35) where g (0) (X) = 1, g (1) (X) = X, g (2) , and so on. It is clear that m (0) corresponds to the total mass of the adsorbed particle, m (1) /m (0) is the center of mass of the particle or X 0 , and m (2) /m (0) is related to the size of the adsorbed particle projected to the x axis. Higher moments can be calculated in the same way and give information about the shape of the particle. In order to calculate these moments from the relative frequency shifts, the functions g (k) (X) are expressed as linear combinations of the mode shapes as where α (k) n are real numbers. Using Equations (34) and (36), the moments m (k) can be expressed as a function of the relative frequency shift measured experimentally as In practice, it is not possible to measure infinite modes, however, with a finite number of modes, it is possible to form a good approximation of the functions g (k) (X). Of course, the accuracy of the calculation will be higher as more modes are used and the approximation of the functions g (k) (X) is better.

Hydrodynamic Loading
As stated before, fundamental science and novel applications in technology still need the development of sensors with higher sensitivity and selectivity than the current state-of-the-art accomplishments. In that context, the advances in nanofabrication tools have placed the nanomechanical resonators as the crucial elements of sensing techniques, reaching unprecedented sensitivities. The application idea behind these outstanding sensitivities is the measurement of the changes induced in the mechanical parameters, resonant frequencies, and/or quality factor, of the resonators. However, in order to reach these impressive milestones in the sensing field, researchers are forced to work in highly controlled environments, i.e., ultra-high vacuum conditions and cryogenic temperatures. Hence, the application of these developments is limited to fundamental studies, leaving aside the majority of the biological and chemical reactions that take place in physiological environments. This is because the performance of nanomechanical resonators is heavily deteriorated in liquids, mainly owing to the double side effect of viscous environments [93][94][95][96]; the effective mass of the resonator is dramatically increased because of the dragged mass of the displaced fluid during the oscillation period and the viscous damping degrades the mechanical quality factor. Therefore, in recent years, there has been a great effort by the scientific community to elucidate new ways to combine the excellent dynamical properties of resonators oscillating in vacuum and the new sensing scenarios offered by the sensors operated in liquids. Here, we are going to distinguish between two different operation modes: the operation of the nanomechanical resonator immersed in liquid and the developments of suspended nano-and micro-channels, where the liquid flows inside the resonator, which is placed into a vacuum chamber.

Nanomechanical Resonators Immersed in Fluid
According to previous studies [97][98][99], the frequency dependence of the damping force for an oscillator immersed in a viscous fluid is related to the well-known hydrodynamic force, i.e., the force that the fluid exerts on the oscillator. For a cantilever beam vibrating in a fluid, the dominant geometric parameter is the width and the hydrodynamic force in the frequency domain can be expressed aŝ whereŵ is the out-of-plane displacement; f is the frequency; ρ f is the density of the fluid; and Γ( f ) is the hydrodynamic function that, for a circular cross section, can be obtained analytically [100] and can be written as where i = √ −1, K 0 , and K 1 are the modified Bessel functions of the third kind and re = ρ f f πb 2 /2η f is the Reynolds number, where η f is the viscosity of the fluid. However, for a rectangular cross section, the solution is not analytical. The previous function can be corrected by a geometrical factor Ω(ω) so that the new hydrodynamic function is expressed as the product of the geometrical factor and the hydrodynamic function for a circular cross section: Γ rec ( f ) = Ω( f )Γ( f ). This geometrical factor Ω(ω) was calculated in [99] as a rational function of log 10 re for a range re ∈ 10 −6 , 10 4 . The real and imaginary parts of Ω( f ) can be expressed as where τ = log 10 re. Figure 7a shows the real and imaginary parts of the hydrodynamic function of a rectangular cantilever of 200 µm by 20 µm typically used in sensing applications. immersed in liquid is not the same as the one oscillating in vacuum with a lower frequency and quality factor. Figure 7b shows the resonant peak for a constant quality factor (solid line) and for the quality factor in Equation (47)   The hydrodynamic function Γ ( ) mentioned above was obtained for a cantilever with ≫ . The accuracy of this function decreases as / increases, especially for higher order modes. However, an improvement of the hydrodynamic function was proposed in order to increase the accuracy for these cases [100]. The hydrodynamic function for a finite aspect ratio and higher order modes can be calculated by solving a system of linear equations that depend on the mode order and on the aspect ratio of the cantilever: where the coefficient can be calculated solving the following system of linear equations: where = and the coefficients , and , are given by The movement of the cantilever in the fluid is impeded by two different mechanisms: the added mass (m ad ) of the fluid displaced by the resonator during its oscillation cycle and the damping factor (γ f ). By expanding the hydrodynamic force into its real and imaginary parts, it is possible to obtain both factors [98,101] in terms of the hydrodynamic function: where Re and Im account for the real and imaginary, parts respectively, and in order to calculate Equations (42) and (43), the normalization condition (25) was applied to the mode shape of the cantilever. Considering these results, we can model the cantilever movement immersed in the fluid as a damped harmonic oscillator as follows: where F th (t) is the non-correlated Langevin thermal force and k e f f n is the effective spring constant of the cantilever, given by k e f f n = 1 12 The resonance frequency and quality factor can be approximated by the following expressions [99] where f R,n is the nth resonance frequency of the cantilever immersed in fluid, f vac,n is the corresponding resonance frequency in vacuum, and Q n is the quality factor. Looking at the different parameters in Equation (44), we realize that the mass and the damping factor are not constants anymore; there is a functional dependency with the frequency by means of the hydrodynamic function. Thus, the behavior of a nanomechanical resonator immersed in liquid is not the same as the one oscillating in vacuum with a lower frequency and quality factor. Figure 7b shows the resonant peak for a constant quality factor (solid line) and for the quality factor in Equation (47) (dashed line). The hydrodynamic function Γ rec ( f ) mentioned above was obtained for a cantilever with L b. The accuracy of this function decreases as b/L increases, especially for higher order modes. However, an improvement of the hydrodynamic function was proposed in order to increase the accuracy for these cases [100]. The hydrodynamic function for a finite aspect ratio and higher order modes can be calculated by solving a system of linear equations that depend on the mode order and on the aspect ratio of the cantilever: Γ rec ( f , n) = 8a 1 (48) where the coefficient a 1 can be calculated solving the following system of linear equations: where κ n = λ n b L and the coefficients A κ q,m and A Re q,m are given by A re q,m = −κ n 2 2 4q−5 √ π G 21 13 κ n 2 −i4re 16 0 where G 21 13 are Meijer functions and the integer p is to be increased until convergence is achieved.

Suspended Microchannel Resonators
The current state-of-the-art nanofabrication capabilities allow the development of micro-channels, where the liquid flows inside the resonator, which is placed into a vacuum chamber [38,102]. This is the most promising technique for nanomechanical sensing in liquid environments. Using this configuration, the nanomechanical resonator retains the good mechanical properties of a resonator in vacuum (high frequency and mechanical quality factor), while allowing the sensing in physiological environments, which is crucial for most biological or chemical reactions [103]. These new accomplishments allow to decrease the minimum detectable mass to reach tens of attograms resolution in liquid [104]. However, these hollow resonators are not capable of measuring the analyte inertial mass, giving us the so-called buoyant mass. When an analyte is flowing inside a liquid passing through the hollow resonator, the measured frequency shift is the result of the difference between the mass of the displaced volume of the carrier liquid and the mass of the analyte, which could be problematic when trying to distinguish between two different analytes of different materials and/or volumes. This problem has been addressed in previous studies by measuring the analytes while changing the carrier liquid, which means that it is necessary to measure each analyte several times to unambiguously distinguish between two analyte populations. Nevertheless, this technique remains costly, slow, and not user friendly, because, far from being a universal method, for biological analysis, both carrier liquids need to be carefully selected depending on the specific sample.
Recently, transparent hollow capillary resonators have been introduced to solve this problem by adding a new source of information: the optical properties of the flowing analytes [105]. Using a laser beam, it is possible to optically probe the analytes while measuring the mechanical displacement of the resonator [56]. In this technique, an interferometric measurement is used: light bounces between the suspended capillary and the substrate underneath, which provides an extraordinary displacement sensitivity [106]. The phase difference between the light reflected by the vibrating suspended capillary and the substrate provides the needed intensity modulation in the recombined beam reaching a photodetector. This approach uses a doubly clamped suspended resonator (Figure 8a) instead of a single clamped cantilever-like structure. The laser is focused on the center of the suspended area, where the displacement amplitude is maximal. As the phase difference between the reflected beams is not exactly π/2, the light collected by the photodetector carries double-sided information; that is, a harmonic intensity modulation given by interference caused by the mechanical displacement of the resonator (named as AC signal) over a non-modulated signal, which is proportional to the device reflectivity (named as DC signal). When a particle flows through the capillary, the mechanical resonance frequency consequently decreases owing to the added mass as well as the reflectivity when passing through the illuminated area as a result of the light scattering. From the reflectivity change (DC signal), it is possible to calculate the particle optical absorption, which depends on the material properties. However, in order to calculate intensive parameters such as the particle mass density or the refractive index, it is necessary to measure the particle size. In order to calculate it, it is necessary to know the velocity of the analytes passing through the device. From the mechanical frequency shift (AC signal) as a function of time ( Figure  8b), the velocity of the particle is calculated by dividing the length of the suspended area (fixed by the device design) by the time it takes to recover the original frequency. The particle size is subsequently obtained by analyzing the optical scattering signal (DC signal) as a function of time. Whenever the flowing particle is illuminated by the laser spot, its edges scatter light, which is measured as a reflectivity signal drop. Thus, the diameter of the particle is finally calculated as the measured velocity multiplied by the elapsed time between two consecutive scattering signal drops. The concept of buoyant mass comes from the Archimedes' principle: an object immersed in liquid experiences a force equal to the weight of the displaced fluid. Thus, the buoyant mass of a spherical particle is defined as [107] = − Hence, attending to the particle material and dimensions, it is possible to have particles of different materials and dimensions with the same buoyant mass. The problem of the buoyant mass becomes even more dramatic in the case of disperse particle populations given that two populations of different particles may be mechanically undiscernible, despite having a different average buoyant mass in the case wherein both distributions overlap. As an analyte flows through the suspended area, it produces a frequency shift given by where is the mass of the resonator filled with fluid and ( ) is the mode shape of a doubly clamped beam given by Equation (23). For the data analysis, the dips in the frequency signal produced by flowing analytes can be fitted to the shape of the mode by assuming that particles maintain a constant velocity while passing through the suspended region. These fits not only allow obtaining the frequency shift, but also give the transit time for each particle, which allows for directly calculating the velocity of the particle. Suspended microchannel resonators have also been proposed to measure the compressibility of the particles that flow through them. Biological particles may have a high compressibility contrast with respect to the surrounding fluid; therefore, accessing the compressibility of the particle can add very valuable information for subsequent identification and analysis. Using an opto-mechano-fluidic resonator, Bahl et al. managed to measure the compressibility of small particles based on the change in the frequency asso- The concept of buoyant mass comes from the Archimedes' principle: an object immersed in liquid experiences a force equal to the weight of the displaced fluid. Thus, the buoyant mass of a spherical particle is defined as [107] Hence, attending to the particle material and dimensions, it is possible to have particles of different materials and dimensions with the same buoyant mass. The problem of the buoyant mass becomes even more dramatic in the case of disperse particle populations given that two populations of different particles may be mechanically undiscernible, despite having a different average buoyant mass in the case wherein both distributions overlap. As an analyte flows through the suspended area, it produces a frequency shift given by where m 0 is the mass of the resonator filled with fluid and ψ n (X) is the mode shape of a doubly clamped beam given by Equation (23). For the data analysis, the dips in the frequency signal produced by flowing analytes can be fitted to the shape of the mode by assuming that particles maintain a constant velocity while passing through the suspended region. These fits not only allow obtaining the frequency shift, but also give the transit time for each particle, which allows for directly calculating the velocity of the particle. Suspended microchannel resonators have also been proposed to measure the compressibility of the particles that flow through them. Biological particles may have a high compressibility contrast with respect to the surrounding fluid; therefore, accessing the compressibility of the particle can add very valuable information for subsequent identification and analysis. Using an opto-mechano-fluidic resonator, Bahl et al. managed to measure the compressibility of small particles based on the change in the frequency associated with radi-ally symmetric mechanical breathing modes of the resonator as the particle passes through it [108]. The formula that they proposed for the relative frequency shift contains the inertial term due to the mass and a new term including compression effects of the particle: where χ is the compressibility and subscripts a and f refer to analyte and fluid, respectively. The coefficients C and D provide the sensitivity to the compressibility and density contrasts between the liquid flowing inside the channel and the solid particle, and are given by where W ap and W ak are the acoustic potential and kinetic energies, respectively, in the fluid evaluated over the volume displaced by the particle. W f p and W f k are the acoustic potential and kinetic energies, respectively, in the fluid evaluated over the entire volume of the fluid, and W wp and W wk are the strain and kinetic energies, respectively, associated with the solid part of the resonator. Their capability to simultaneously access the analytes' buoyant mass and compressibility, together with their extremely high speed, which enables characterizing a larger amount of analytes in a shorter time, settles optomechanical capillaries as one of the most promising devices in the field. The combination of optical and mechanical sensing methods in channeled platforms will revolutionize not only biology, biomedicine, and clinical diagnosis fields, but also the food industry and environmental analysis. However, in order to achieve success in accessing more and more precise information about the analytes, experimentalists still require more accurate theoretical models.

Summary
Nanomechanical sensors are rapidly gaining more and more attention from researchers and have been proven to be extraordinarily useful tools to monitor biological processes and to detect and measure the fundamental properties of biological entities. The advancements and developments of these systems are possible thanks to theories that allow the interpretation and understanding of the experimental observations. In this review, we have explored the main theories and models that are the base of nanomechanical sensors and that are used on a daily basis by researchers all over the world regarding this field. We have approached the topic by separating nanomechanical sensors into two main groups: static and dynamic. For the nanomechanical sensors in static mode, we have presented the main theories regarding the effect of differential surface stress on the bending of a rectangular cantilever plate, starting from the well-known Stoney's equation, the improved solution proposed by Sader et al., and the simplified averaged curvatures expressions obtained by Tamayo et al. The effect of surface stress can also modify the dynamic properties of the sensor, allowing for quantifying the stress by measuring changes in the resonance frequencies of the nanomechanical resonator. This effect has been covered as well in this review, where we have mentioned both the model to explain the effect of net surface stress on cantilevers and doubly clamped beams and the model to explain the effect of differential surface stress on cantilever beams. Other features of the biological systems that can produce a shift in the natural frequencies of nanomechanical resonators are the mass and the stiffness. When the biological system is attached to the resonator surface, the resonant frequencies of the sensor change as a result of two opposite effects: the added mass decreases the frequencies, while the stiffness effect increases the frequencies. We have presented the main theoretical models for the effect of the mass and stiffness of analytes on the resonant frequencies of nanomechanical resonators, where we have emphasized the difference between the ad-sorption of a complete layer and the adsorption of individual particles. While the relative frequency shift is the same for all flexural modes in the case of a complete layer, for the individual particles, the relative frequency shift is different for different modes because it depends on the adsorption position. This particularity led to the multimode measuring technique that consists of measuring the relative frequency shifts of different vibration modes simultaneously in order to extract the mass, stiffness, and adsorption position of the attached particle. This process is known as the inverse problem and is the base for nanomechanical mass spectrometry. Multimode measuring can also be used to obtain the size and shape of a particle that cannot be considered as a point-like particle. The formalism to do so is called inertial imaging and has also been included in this review. Finally, in the last part of the review, we focused the attention on the effects of hydrodynamic loading when the sensor is immersed in a viscous fluid. This part is especially important when the nanomechanical sensor is used for biological applications because water is the natural environment of biology. We reviewed the main theories and models that describe the effects of hydrodynamic loading on the frequency and quality factor of cantilevers, as well as the new microchannel resonators that overcome the degradation of the quality factor as a result of viscous damping for sensing applications in liquid environments. To conclude, we would like to emphasize that theory and modeling are fundamental for the development and advancements of nanomechanical sensors, as they not only allow for the understanding of the phenomena, but also help to optimize the devices for a particular application. It is then a must that researchers can keep advancing the existing models and creating new ones in order to keep up the growth of this fascinating field.

Data Availability Statement:
The data presented in this review are available from the corresponding author upon reasonable request.