# Multi-mode Brownian Dynamics of a Nanomechanical Resonator in a Viscous
Fluid

H. Gress  Department of Mechanical Engineering, Division of Materials Science
and Engineering, and the Photonics Center, Boston University, Boston,
Massachusetts 02215, USA  J. Barbish  Department of Mechanical Engineering,
Virginia Tech, Blacksburg, Virginia 24061, USA  C. Yanik  SUNUM,
Nanotechnology Research and Application Center, Sabanci University, Istanbul,
34956, Turkey  Faculty of Engineering and Natural Sciences, Sabanci
University, Istanbul, 34956, Turkey  I. I. Kaya  SUNUM, Nanotechnology
Research and Application Center, Sabanci University, Istanbul, 34956, Turkey
Faculty of Engineering and Natural Sciences, Sabanci University, Istanbul,
34956, Turkey  R. T. Erdogan  Department of Mechanical Engineering, Bilkent
University, Ankara, 06800, Turkey  National Nanotechnology Research Center
(UNAM), Bilkent University, Ankara, 06800, Turkey  M. S. Hanay  Department of
Mechanical Engineering, Bilkent University, Ankara, 06800, Turkey  National
Nanotechnology Research Center (UNAM), Bilkent University, Ankara, 06800,
Turkey  M. González  Aramco Americas, Aramco Research Center–Houston, Houston,
Texas 77084, USA  O. Svitelskiy  Department of Physics, Gordon College,
Wenham, Massachusetts 01984, USA  M. R. Paul  Department of Mechanical
Engineering, Virginia Tech, Blacksburg, Virginia 24061, USA  K. L. Ekinci
<EMAIL_ADDRESS>Department of Mechanical Engineering, Division of Materials
Science and Engineering, and the Photonics Center, Boston University, Boston,
Massachusetts 02215, USA

###### Abstract

Brownian motion imposes a hard limit on the overall precision of a
nanomechanical measurement. Here, we present a combined experimental and
theoretical study of the Brownian dynamics of a quintessential nanomechanical
system, a doubly-clamped nanomechanical beam resonator, in a viscous fluid.
Our theoretical approach is based on the fluctuation-dissipation theorem of
statistical mechanics: We determine the dissipation from fluid dynamics; we
incorporate this dissipation into the proper elastic equation to obtain the
equation of motion; the fluctuation-dissipation theorem then directly provides
an analytical expression for the position-dependent power spectral density
(PSD) of the displacement fluctuations of the beam. We compare our theory to
experiments on nanomechanical beams immersed in air and water, and obtain
excellent agreement. Within our experimental parameter range, the Brownian
force noise driving the nanomechanical beam has a colored PSD due to the
“memory” of the fluid; the force noise remains mode-independent and
uncorrelated in space. These conclusions are not only important for
nanomechanical sensing but also provide insight into the fluctuations of
elastic systems at any length scale.

##  I Introduction

Brownian fluctuations of mechanical systems have been a topic of active
research in physics since the 1920s [1, 2]. Early electrometers [3] and
galvanometers [4] that featured proof masses attached to linear springs
displayed irregular movements around their equilibrium points despite all
“precautions and shields” [5]. These early experiments eventually led to the
realization that the observed fluctuations, namely, Brownian motion, were of
fundamental nature and limited the overall precision of mechanical
measurements [5]. A century later, Brownian motion still remains centrally
relevant to precision metrology and sensing based on mechanical systems—in
particular, nanoelectromechanical systems (NEMS) [6] and AFM microcantilevers
[7]. These state-of-the-art miniaturized mechanical systems are even more
susceptible to Brownian noise than their macroscopic counterparts since they
tend to be extremely compliant to forces.

The Brownian dynamics of a nanomechanical resonator can be formulated using
elasticity theory and statistical mechanics. Elasticity theory provides a
dissipationless equation of motion, such as the beam equation. Solving this
equation under a harmonic ansatz and subject to boundary conditions maps the
dynamics of the nanomechanical resonator onto that of a collection of
eigenmodes, i.e., spring-mass systems, with discrete eigen-frequencies and
mode shapes (eigenfunctions) [8, 9]. In the simplest approximation of Brownian
dynamics, each eigenmode is assumed to have a constant and spatially uniform
dissipation, resulting in a Brownian force noise that is delta-function
correlated in both time and space. These assumptions result in a theoretical
expression for the power spectral density (PSD) of the displacement
fluctuations of the beam as a sum of the PSDs of the individual uncorrelated
eigenmode fluctuations [8, 10]. In the limit of small dissipation, multi-mode
noise measurements on cantilevers [11], microdiscs [12], microtoroids [13],
nanowires [14] and macroscopic elastic systems [15] all agree with this first-
order approximation.

Most mechanical systems, however, do come with some “memory,” making the
above-mentioned assumption of temporally uncorrelated force noise inaccurate
[8]. If the dissipation is spatially non-uniform, i.e., position dependent,
the force noise between different eigenmodes becomes correlated [16], with the
eigenmode expansion of the force noise becoming non-trivial. For elastic
systems with spatially non-uniform dissipation, alternative approaches to
calculate the noise PSD have been developed [16, 17] and experimentally tested
[18, 19, 20].

For a nanomechanical resonator immersed in a viscous fluid, memory comes from
the flow-resonator interaction [21]. Here, the presence of the viscous fluid
allows for a viable path to formulate the Brownian dynamics of the
nanomechanical resonator consistently [22, 23, 24]: first, the dissipation of
the resonator is found from fluid dynamics; then, the fluctuation-dissipation
theorem is used for the calculation of the PSD of the resonator fluctuations.
Since the fluidic dissipation is frequency dependent, the force noise PSD is
also “colored,” and the resonator fluctuations deviate substantially from the
first-order approximation discussed above. In nearly all work so far, the
dissipation in viscous fluids has been assumed to be mode independent and
spatially uniform, resulting in a Brownian force noise that is spatially
uncorrelated. This assumption of spatial homogeneity again leads to formulas
expressible as a sum in terms of the individual eigenmodes. There are notable
papers, where the experimental noise data have been successfully fitted with
colored PSDs. However, these experiments typically do not extend beyond the
first mode of the elastic structure [25, 26, 27, 28] and are thus not very
insightful on the spatial nature of the force noise.

The topic of this manuscript is the Brownian dynamics of a nanomechanical
resonator in a viscous fluid. In particular, we investigate how the Brownian
force driving the nanomechanical resonator is correlated in time and space in
a viscous fluid. To this end, we derive an expression for the PSD of the
displacement fluctuations of an elastic nanomechanical beam under tension in a
viscous fluid, assuming a frequency-dependent but spatially homogeneous
viscous dissipation. This results in a PSD that is the sum of the PSDs of the
uncorrelated fluctuations of individual eigenmodes. We validate this theory by
experiments performed on nanomechanical beams under tension immersed in air
and water. Using solely experimental parameters, we obtain excellent agreement
between the experimental data and theory. This agreement, up to the twelfth
eigenmode in air and seventh eigenmode in water, validates our overarching
assumptions: the Brownian force noise has a colored PSD due to the memory of
the fluid but can be approximated to be uncorrelated in space.

Figure 1: (a) A doubly-clamped beam with length $L$, width $b$, thickness $h$,
and gap $g$ between the beam and substrate. (b) Illustration of the first
three eigenfunctions $\phi_{n}(x)$. (c) SEM image of a doubly-clamped silicon
nitride beam under tension with $L\times b\times h\\!=\\!50~{}\rm\mu m\times
900~{}nm\times 93~{}nm$ and $g\\!=\\!2~{}\rm\mu m$. The silicon nitride layer
is shown in green, with the suspended part in light green. The silicon
substrate and the gold layer are gray and yellow, respectively.

##  II Theory

We start with the equation of motion for a beam under tension driven by a
deterministic external force in a viscous fluid. The respective linear
dimensions of the beam are $L\times b\times h$ along the $xyz$ axes [Fig.
1(a)], and $\mu=\rho_{s}bh$ is the mass per unit length with $\rho_{s}$ being
the density. The flexural displacement of the beam, $W(x,t)$, along the $z$
axis at position $x$ and time $t$ is given by

$\frac{EI}{L^{4}}\frac{\partial^{4}W(x,t)}{\partial
x^{4}}-\frac{F_{T}}{L^{2}}\frac{\partial^{2}W(x,t)}{\partial
x^{2}}+\mu\frac{\partial^{2}W(x,t)}{\partial t^{2}}=\\\
F_{f}(x,t)+F_{d}(x,t).$ (1)

The $x$ coordinate has been normalized by $L$ such that $0\\!\leq x\\!\leq 1$.
In Eq. (1), $E$ is the Young’s modulus, $I$ is the area moment of inertia,
$F_{f}(x,t)$ is the force per unit length of the fluid acting on the beam,
$F_{d}(x,t)$ is the external drive force per unit length, and $F_{T}$ is the
tension [29]. The beam has fixed boundaries such that
$W(0,t)\\!=\\!W(1,t)\\!=\\!W^{\prime}(0,t)\\!=\\!W^{\prime}(1,t)=\\!0$, where
a prime indicates an $x$ derivative.

It will be useful to proceed in the frequency domain [30] using the Fourier
transform pair

$\displaystyle\hat{W}(x,\omega)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}W(x,t)e^{i\omega t}dt,$ (2)
$\displaystyle W(x,t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{W}(x,\omega)e^{-i\omega
t}d\omega,$ (3)

where $\omega$ is the angular frequency. This leads to the transformed
differential equation

$\frac{EI}{L^{4}}\frac{\partial^{4}\hat{W}(x,\omega)}{\partial
x^{4}}-\frac{F_{T}}{L^{2}}\frac{\partial^{2}\hat{W}(x,\omega)}{\partial
x^{2}}\\\
-\omega^{2}\left[\mu+\frac{\pi}{4}\rho_{f}b^{2}\Gamma(\omega)\right]\hat{W}(x,\omega)=\hat{F}_{d}(x,\omega)$
(4)

for the Fourier component $\hat{W}(x,\omega)$ at $\omega$. In Eq. (4), we have
described the force due to the fluid as

$\hat{F}_{f}(x,\omega)=\frac{\pi}{4}\rho_{f}\omega^{2}b^{2}\Gamma(\omega)\hat{W}(x,\omega),$
(5)

where $\rho_{f}$ is the fluid density and $\Gamma(\omega)$ is the complex
hydrodynamic function for a blade, derived from Stokes’ oscillating cylinder
theory [31, 32, 33]. In Eq. (5), $\Gamma(\omega)$ quantifies the mass loading
and viscous damping of the fluid acting on the beams. To obtain
$\Gamma(\omega)$, one starts with the complex hydrodynamic function for an
infinitely-long oscillating cylinder,

$\Gamma_{c}\left(\rm{Re}_{\omega}\right)=1+\frac{4iK_{1}\left(-i\sqrt{i\rm{Re}_{\omega}}\right)}{\sqrt{i\rm{Re}_{\omega}}K_{0}\left(-i\sqrt{i\rm{Re}_{\omega}}\right)},$
(6)

where $K_{0}$ and $K_{1}$ are respectively the zeroth and first order modified
Bessel functions of the second kind [34]. The argument of $\Gamma_{c}$ is the
frequency-dependent Reynolds number, ${\rm{Re}}_{\omega}=\frac{\rho_{f}\omega
b^{2}}{4\eta_{f}}$, where $\eta_{f}$ is the dynamic viscosity of the fluid. To
account for the rectangular cross-section of the beams, one then applies a
small frequency-dependent correction factor to $\Gamma_{c}$ [33]. It can be
deduced from Eq. (6) that the only parameters in $\Gamma(\omega)$ are $b$ and
$\eta_{f}$, which are both constants. Thus, $\Gamma(\omega)$ is assumed to be
independent of position $x$ as well as the mode-shape of the beam.

We solve Eq. (4) using the eigenfunction expansion [33],

$\hat{W}(x,\omega)=\sum_{n=1}^{\infty}f_{n}(\omega)\phi_{n}(x),$ (7)

where $n$ is the mode number, $f_{n}(\omega)$ describes the frequency
dependence, and $\phi_{n}(x)$ are the orthonormal eigenfunctions of the beam
with tension [Fig. 1(b)]. Expressions for $\phi_{n}(x)$ and the eigen-
frequencies $\omega_{n}/2\pi$ for a doubly clamped beam with tension are
available [30, 35, 36, 28]. We note that both $\phi_{n}(x)$ and
$\omega_{n}/2\pi$ are found from the dissipationless equation of motion. The
influence of the tension force on $\phi_{n}(x)$ and $\omega_{n}$ can be
quantified in terms of the nondimensional tension parameter $U$, where
$U=\frac{F_{T}}{2EI/L^{2}}$. The dynamics respectively becomes that of an
Euler-Bernoulli beam and a string for $U\to 0$ and $U\gg 1$.

Using the orthogonality of $\phi_{n}(x)$, the solution to Eq. (4) can be
expressed as

$\hat{W}(x,\omega)=\frac{L^{4}}{EI}\sum_{n=1}^{\infty}\frac{\int_{0}^{1}\hat{F}_{d}(x^{\prime},\omega)\phi_{n}(x^{\prime})dx^{\prime}}{\Omega_{n}^{2}-B^{4}(\omega)}\phi_{n}(x).$
(8)

$\Omega_{n}$ are the nondimensional eigen-frequencies defined as

$\Omega_{n}\\!=\\!\frac{\omega_{n}}{\alpha/L^{2}},$ (9)

where $\alpha\\!=\\!(EI/\mu)^{1/2}$. The complex function $B(\omega)$ contains
the dissipation and added mass, and is given by

$B^{4}(\omega)=\Omega_{1}^{2}\left(\frac{\omega}{\omega_{1}}\right)^{2}\left[1+T_{0}\Gamma(\omega)\right],$
(10)

where $\omega_{1}$ is the fundamental eigen-frequency of the beam in the
absence of fluid, i.e., dissipation and added mass. The mass loading
parameter, $T_{0}\\!=\\!\frac{\pi}{4}\frac{\rho_{f}b}{\rho_{s}h}$, is the
ratio of the mass of a cylinder of fluid with diameter $b$ to the mass of the
beam.

In order to connect with the fluctuation-dissipation theorem, we next
calculate the susceptibility, ${\chi}(x_{0},t)$, which we define as the time-
dependent displacement of the beam measured at position $x_{0}$ due to the
application of a unit impulse of force at the same position $x_{0}$. Thus, we
specify

$F_{d}(x,t)=\frac{1}{L}\delta(x-x_{0})\delta(t),$ (11)

which becomes

$\hat{F}_{d}(x,\omega)=\frac{1}{L}\delta(x-x_{0})$ (12)

in the frequency domain, with $\delta$ being the Dirac delta function. Hence,
$\hat{\chi}(x_{0},\omega)=\hat{W}(x_{0},\omega)$, and we obtain

$\hat{\chi}(x_{0},\omega)=\frac{L^{3}}{EI}\sum_{n=1}^{\infty}\frac{\int_{0}^{1}\delta(x^{\prime}-x_{0})\phi_{n}(x^{\prime})dx^{\prime}}{\Omega_{n}^{2}-B^{4}(\omega)}\phi_{n}(x_{0}),$
(13)

which can be expressed as

$\hat{\chi}(x_{0},\omega)=\frac{L^{3}}{EI}\sum_{n=1}^{\infty}\frac{\phi_{n}^{2}(x_{0})}{\Omega_{n}^{2}-B^{4}(\omega)}.$
(14)

Using $\omega_{n}/\omega_{1}=\Omega_{n}/\Omega_{1}$, defining
$\tilde{\omega}_{n}=\omega/\omega_{n}$, and simplifying further yields

$\hat{\chi}(x_{0},\omega)=\sum_{n=1}^{\infty}\frac{1}{k_{n}(x_{0})}\\\
\times\frac{1}{1-\tilde{\omega}_{n}^{2}(1+T_{0}\Gamma^{\prime}(\omega))-i\tilde{\omega}_{n}^{2}T_{0}\Gamma^{\prime\prime}(\omega)},$
(15)

where $\Gamma^{\prime}(\omega)$ and $\Gamma^{\prime\prime}(\omega)$ are the
real and imaginary parts of $\Gamma(\omega)$, respectively. The effective
spring constant of mode $n$, when measured at $x_{0}$, is represented as
$k_{n}(x_{0})$; $k_{n}(x_{0})$ can be consistently determined by ensuring that
the kinetic energy of the spatially extended oscillating beam with mode shape
$\phi_{n}(x)$ equals that of a lumped system measured at $x_{0}$. This yields

$k_{n}(x_{0})=\frac{m\omega_{n}^{2}}{\phi_{n}^{2}(x_{0})},$ (16)

where $m=\mu L$ is the nominal mass of the beam.

The PSD of the Brownian fluctuations of the beam at axial position $x_{0}$ can
be directly expressed using the fluctuation-dissipation theorem [37, 38] as

$G_{W}(x_{0},\omega)=\frac{4k_{B}T}{\omega}\hat{\chi}^{\prime\prime}(x_{0},\omega),$
(17)

where $\hat{\chi}^{\prime\prime}(x_{0},\omega)$ is the imaginary part of
$\hat{\chi}(x_{0},\omega)$, $k_{B}$ is Boltzmann’s constant, and $T$ is the
temperature. The subscript $W$ on $G_{W}$ indicates that this is the spectral
density of the fluctuations in the flexural displacement $W(x,t)$ of the beam
along the $z$ axis.

Taking the imaginary part of Eq. (15) to find
$\hat{\chi}^{\prime\prime}(x_{0},\omega)$ and inserting this into Eq. (17), we
obtain the desired result as

$G_{W}(x_{0},\omega)=4k_{B}T\sum_{n=1}^{\infty}\frac{1}{k_{n}(x_{0})\omega_{n}}\\\
\times\frac{\tilde{\omega}_{n}T_{0}\Gamma^{\prime\prime}(\omega)}{\left[1-\tilde{\omega}_{n}^{2}(1+T_{0}\Gamma^{\prime}(\omega))\right]^{2}+\left[\tilde{\omega}_{n}^{2}T_{0}\Gamma^{\prime\prime}(\omega)\right]^{2}}.$
(18)

This expression yields the total PSD for the displacement fluctuations of the
beam at frequency $\omega$ and position $x_{0}$. We emphasize that
$G_{W}(x_{0},\omega)$ in Eq. (18) is obtained as a sum over individual
eigenmodes. This is because $\Gamma(\omega)$ is assumed to be spatially
homogeneous and independent of mode number $n$ [39].

Figure 2: PSDs of the displacement fluctuations $G_{W}(x_{n},\omega)$ plotted
using semi-logarithmic axes for the first twelve modes of a beam with
$L=30~{}\rm\mu m$ in air. The measurement is performed at $x\\!=\\!x_{n}$,
i.e., an antinode of mode $n$. The continuous lines are fits based on Eq.
(18). The insets for $n\\!=\\!9$ to $12$ show $G_{W}(x_{n},\omega)$ using
linear axes. Because the peak of the twelfth mode is barely resolved, we do
not determine values for $\omega_{12}/2\pi$, $k_{12}(x_{12})$, and
$Q_{12}^{(a)}$.

##  III Experiments

Our experiments are performed on silicon nitride doubly-clamped beams of
$b\\!=\\!900~{}\rm nm$, $h\\!=\\!93~{}\rm nm$, and three different lengths of
$L\\!=\\!30,~{}40$ and $50~{}{\rm\mu m}$; there is a gap of
$g\\!=\\!2~{}\rm\mu m$ between each beam and the substrate. Fig. 1(c) shows a
scanning electron microscope (SEM) image of a beam with $L\\!=\\!50~{}\rm\mu
m$. All the beams are from the same fabrication batch. The beams are under
tension as inferred from their resonance frequencies in vacuum [29]. We
determine $\mu\\!=\\!\rho_{s}bh\\!=\\!2.66\\!\times\\!10^{-10}~{}\rm kg/m$,
with the density measured as $\rho_{s}\\!=\\!2960~{}\rm kg/m^{3}$ [29]. The
beams also have u-shaped gold nanoresistor patterns near their anchors for
other experiments [28, 29].

We measure the displacement fluctuations of the beams using a path-stabilized
homodyne Michelson interferometer. The diffraction limited HeNe laser spot
with a diameter of $\sim 650\pm 10$ nm (FWHM) is positioned on the beam using
an XYZ precision stage. The typical powers incident on the beam and the
photodetector are $\sim 600~{}\rm\mu W$ and $\sim 1~{}\rm mW$, respectively,
with a shot noise limited displacement sensitivity of $\sim 5~{}\rm
fm/\sqrt{Hz}$. We calibrate the system against the wavelength of the laser
[40] and operate at the point of optimal sensitivity. For each measurement
taken at a given position $x_{0}$ on a beam, a second measurement is taken
with the same parameters near the anchor of the beam to determine the
background noise level in our measurements, e.g., due to the low-frequency
laser noise or cable resonances. We assume that the beam fluctuations and the
background noise are uncorrelated, and subtract this background from the noise
PSD measured on the beam [41]. This allows us to resolve the beam fluctuations
down to $\sim 2~{}\rm fm/\sqrt{Hz}$.

Figure 3: PSDs of the displacement fluctuations $G_{W}(x_{n},\omega)$ plotted
using semi-logarithmic axes for the first twelve modes of a beam with
$L=40~{}\rm\mu m$ in air. The measurement is performed at $x\\!=\\!x_{n}$,
i.e., an antinode of mode $n$. The continuous lines are fits based on Eq.
(18). The insets for $n\\!=\\!9$ to $12$ show $G_{W}(x_{n},\omega)$ using
linear axes. Figure 4: PSDs of the displacement fluctuations
$G_{W}(x_{n},\omega)$ plotted using semi-logarithmic axes for the first twelve
modes of a beam with $L=50~{}\rm\mu m$ in air. The measurement is performed at
$x\\!=\\!x_{n}$, i.e., an antinode of mode $n$. The continuous lines are fits
based on Eq. (18). The insets for $n\\!=\\!9$ to $12$ show
$G_{W}(x_{n},\omega)$ using linear axes.

The finite size of the optical spot introduces errors into the measurements
[42]. A significant source of error is the curvature of the beam, especially
in higher modes. By computing the overlap of the Gaussian optical spot with
the beam mode, we estimate this error to be less than $10\%$ for the twelfth
mode of our $30$-$\rm\mu m$-long beam, which has the largest curvature in all
our experiments. In water measurements, it also becomes problematic to
position the optical spot precisely at the desired locations on the beam. To
find $x_{0}=0.50$ (center) and $x_{0}=0.25$ ($L/4$) positions on the beam, we
maximize the $n=1$ and $n=2$ peaks, respectively.

##  IV Results

###  IV.1 Vacuum

We first measure the resonance frequencies of the eigenmodes of the beams in
vacuum ($p<10^{-6}~{}\rm bar$), as listed in Table 1. These peak frequencies
in vacuum, $\omega_{n}^{(v)}/2\pi$, should be very close to the eigen-
frequencies of the dissipationless beam, $\omega_{n}/2\pi$. We also extract
quality factors in vacuum from Lorentzian fits and find that all
$Q_{n}^{(v)}\gtrsim 10^{3}$.

Figure 5: (a) Eigen-frequencies $\omega_{n}/2\pi$ and (b) effective spring
constants $k_{n}(x_{n})$ for beams with $L=30~{}\rm\mu m$, $L=40~{}\rm\mu m$,
and $L=50~{}\rm\mu m$. Experimental values and theoretical predictions are
shown by symbols and continuous lines, respectively. Typical error bars for
$k_{n}(x_{n})$ are smaller than the symbols unless shown explicitly. Figure 6:
PSDs of the displacement fluctuations $G_{W}(x_{0},\omega)$ at $x_{0}=0.50$
(a-c) and $x_{0}=0.25$ (d-f) for three beams with different lengths in water.
The continuous lines are predictions of Eq. (18) based on experimental
parameters. The dotted lines show the contributions of individual modes $n$ to
$G_{W}(x_{0},\omega)$. Modes not shown do not contribute significantly at the
particular $x_{0}$.

###  IV.2 Air

Next, we examine the displacement fluctuations in air. We first identify the
frequency at which our system transitions from viscous to molecular flow. For
this system, the transition frequency $\omega_{c}/2\pi$ can be found using
${\omega_{c}\tau}+{\lambda\over b}\approx 1$, where $\tau$ is the relaxation
time and $\lambda$ is the mean free path in the fluid [41, 43]. In air,
$\tau\approx 1~{}\rm ns$ and $\lambda\approx 68~{}\rm nm$ [41, 43]. The
relevant length scale, $b=900~{}\rm nm$, is the same for all our devices. We
thus find $\omega_{c}/2\pi\approx 100~{}\rm MHz$ in air. Therefore, our
measurements are mostly within the viscous regime.

Figs. 2, 3, and 4 respectively show the PSDs of the first twelve modes of
beams with $L=30,40$ and $50$ $\rm\mu m$. The low dissipation in air results
in distinctly separated peaks in the PSDs. Each PSD is measured at an antinode
($x=x_{n}$) of mode $n$ as a function of frequency near the peak frequency,
$\omega_{n}^{(a)}/2\pi$, in air. The relatively large modal quality factors,
$17\leq Q_{n}^{(a)}\leq 200$, allow us to determine the eigen-frequencies
$\omega_{n}/2\pi$, the effective spring constants $k_{n}(x_{n})$, and the
theoretical $G_{W}(x_{n},\omega)$ curves in a self-consistent manner. To this
end, we first use the equipartition of energy to find $k_{n}(x_{n})$ from
$\frac{1}{2}k_{n}(x_{n})\left<W_{n}^{2}(x_{n})\right>\\!=\\!\frac{1}{2}k_{B}T$,
where the mean-squared fluctuation amplitude, $\left<W_{n}^{2}(x_{n})\right>$,
is the numerical integral of the experimental $G_{W}(x_{n},\omega)$ data over
frequency. To find the theoretical $G_{W}(x_{n},\omega)$ curve, we calculate
$T_{0}$ and $\Gamma(\omega)$ using the density and viscosity of air at room
temperature. We then insert the experimental $k_{n}(x_{n})$ values along with
$T_{0}$ and $\Gamma(\omega)$ into Eq. (18) and calculate
$G_{W}(x_{n},\omega)$, treating $\omega_{n}$ as a fit parameter. The best fits
are shown as continuous lines in Figs. 2, 3, and 4.

We re-emphasize that, except for $\omega_{n}$, the fits (continuous lines) in
Figs. 2, 3, and 4 are solely determined using known or measured quantities of
the beam and the surrounding fluid. The values of $\omega_{n}/2\pi$ found by
fitting are typically slightly higher than $\omega_{n}^{(a)}/2\pi$ (Table 1).
This is expected because $\omega_{n}/2\pi$ are the eigen-frequencies without
any fluid loading. Our vacuum measurements support this observation (Table 1):
the vacuum frequencies $\omega_{n}^{(v)}/2\pi$ are all slightly larger than
$\omega_{n}^{(a)}/2\pi$. There are very small discrepancies (typically
$<1.5\%$) between the eigen-frequencies determined by fitting, i.e.,
$\omega_{n}/2\pi$, and those obtained from vacuum measurements,
$\omega_{n}^{(v)}/2\pi$. We attribute these small discrepancies to the fact
that the eigen-frequencies of nanomechanical resonators are easily perturbed
by external factors, e.g., the accumulation of adsorbates or changes in the
temperature and humidity of the environment. Since the resonance is very
sharply peaked in air ($Q_{n}^{(a)}>10$) and in vacuum ($Q_{n}^{(v)}\sim
10^{3}$), these perturbations result in small but noticeable frequency shifts
from measurement to measurement.

Table 1: Eigen-frequencies $\omega_{n}/2\pi$ are from fitting our air measurements to our model; peak frequencies, $\omega_{n}^{(v)}/2\pi$ and $\omega_{n}^{(a)}/2\pi$, are measured in vacuum and air, respectively. The effective stiffness $k_{n}(x_{n})$ is found at an antinode, and quality factors $Q_{n}^{(a)}$ in air are found from Lorentzian fits. All frequency values are in MHz. The modal peak for $n=12$ is barely resolvable in the PSD [Fig. 2]. Mode | $L=30~{}\rm\mu m$ | $L=40~{}\rm\mu m$ | $L=50~{}\rm\mu m$
---|---|---|---
|  $\omega_{n}\over 2\pi$ |  $\omega_{n}^{(v)}\over 2\pi$ |  $\omega_{n}^{(a)}\over 2\pi$ |  $k_{n}(x_{n})$ |  $Q_{n}^{(a)}$ |  $\omega_{n}\over 2\pi$ |  $\omega_{n}^{(v)}\over 2\pi$ |  $\omega_{n}^{(a)}\over 2\pi$ |  $k_{n}(x_{n})$ |  $Q^{(a)}$ |  $\omega_{n}\over 2\pi$ |  $\omega_{n}^{(v)}\over 2\pi$ |  $\omega_{n}^{(a)}\over 2\pi$ |  $k_{n}(x_{n})$ |  $Q^{(a)}$
1 |  3.531 |  3.538 |  3.501 |  1.42 |  28 |  2.577 |  2.596 |  2.553 |  1.16 |  22 |  2.017 |  2.019 |  1.995 |  1.04 |  17
2 |  7.158 |  7.224 |  7.111 |  6.05 |  50 |  5.182 |  5.310 |  5.145 |  3.86 |  40 |  4.065 |  4.073 |  4.030 |  3.81 |  35
3 |  10.998 |  11.009 |  10.934 |  21.42 |  69 |  7.906 |  7.997 |  7.856 |  10.94 |  56 |  6.194 |  6.168 |  6.151 |  7.99 |  46
4 |  15.260 |  15.293 |  15.180 |  34.65 |  86 |  10.817 |  10.968 |  10.755 |  25.46 |  70 |  8.356 |  8.326 |  8.303 |  19.11 |  58
5 |  20.159 |  20.141 |  20.067 |  63.28 |  104 |  13.924 |  14.073 |  13.851 |  39.62 |  83 |  10.649 |  10.598 |  10.585 |  29.64 |  68
6 |  25.532 |  25.605 |  25.420 |  118.67 |  121 |  17.373 |  17.537 |  17.289 |  50.77 |  98 |  13.166 |  13.227 |  13.094 |  54.49 |  81
7 |  31.556 |  31.461 |  31.423 |  183.02 |  125 |  21.037 |  21.121 |  20.938 |  93.52 |  110 |  15.717 |  15.725 |  15.636 |  64.15 |  92
8 |  37.940 |  37.784 |  37.789 |  405.19 |  141 |  25.017 |  25.175 |  24.910 |  145.15 |  120 |  18.579 |  18.598 |  18.487 |  83.14 |  102
9 |  45.306 |  45.086 |  45.132 |  495.15 |  150 |  29.275 |  29.354 |  29.156 |  173.99 |  131 |  21.562 |  21.564 |  21.462 |  153.41 |  112
10 |  52.195 |  52.392 |  52.004 |  699.80 |  155 |  34.043 |  34.265 |  33.904 |  375.16 |  137 |  24.682 |  24.695 |  24.573 |  220.33 |  119
11 |  60.014 |  59.904 |  59.802 |  1077.36 |  167 |  38.970 |  38.879 |  38.817 |  375.16 |  144 |  28.112 |  27.988 |  27.991 |  220.33 |  128
12 |  – |  68.342 |  68.293 |  – |  – |  44.291 |  44.234 |  44.120 |  564.47 |  153 |  31.678 |  31.565 |  31.548 |  349.46 |  136

###  IV.3  Eigen-Frequencies and Spring Constants

We next compare the experimental values for $\omega_{n}/2\pi$ and
$k_{n}(x_{n})$ with theoretical predictions of Euler-Bernoulli beam theory
with tension [30]. We estimate the magnitude of the tension force $F_{T}$ by
comparing the eigen-frequencies obtained from experiments, $\omega_{n}/2\pi$,
with those from theory, $\omega_{n}^{(t)}/2\pi$. To determine
$\omega_{n}^{(t)}/2\pi$, we turn to the characteristic equation, which relates
$\omega_{n}^{(t)}/2\pi$ to the unknown tension [35] in terms of the non-
dimensional parameters $\Omega_{n}$ and $U$. In our calculations, we use
nominal beam dimensions as well as $E$ and $\rho_{s}$ values for SiN. Since
reported values for $E$ have a large uncertainty, $200~{}{\rm GPa}\lesssim
E\lesssim 380~{}\rm GPa$ [44, 45, 46], we take the Young’s modulus as
$E=300~{}\rm GPa$. The density has been measured as $\rho_{s}=2960\pm 30~{}\rm
kg/m^{3}$ [29]. For each beam, we sweep the value of $F_{T}$, solve for
$\omega_{n}^{(t)}/2\pi$ for each $F_{T}$, and then compute the error

$\varepsilon=\sum_{n=1}\frac{|\omega_{n}-\omega_{n}^{(t)}|^{2}}{|\omega_{n}|^{2}}$
(19)

between $\omega_{n}$ and $\omega_{n}^{(t)}$, where $n$ encompasses the first
twelve modes. The minimum $\varepsilon$ provides the experimental value of
$F_{T}$. Because all beams are on the same chip, we average the values of
$F_{T}$ for each of our three beams. We thus find $F_{T}=7.43~{}\rm\mu N$. In
Fig. 5(a), we show experimental $\omega_{n}/2\pi$ data (symbols) and
theoretical predictions (continuous lines) using $F_{T}=7.43~{}\rm\mu N$ on
semi-logarithmic axes for all three beams.

To find the theoretical $k_{n}(x_{n})$ at an antinode, we use Eq. (16). To
this end, we calculate $\phi_{n}(x)$ for each beam [30] using $F_{T}$ and
$\omega_{n}/2\pi$; we use nominal $m$, measured $\omega_{n}$, and calculated
$\phi_{n}(x_{n})$ to determine $k_{n}(x_{n})$. We show the experimental and
theoretical values for $k_{n}(x_{n})$ in Fig. 5(b). The experimental spring
constants match predictions closely over two orders of magnitude for
$n\lesssim 7$. With our knowledge of $\phi_{n}(x)$ and $k_{n}(x_{n})$, we can
determine $k_{n}(x_{0})$ for any position $x_{0}$ along the beam.

The differences between experimental and theoretical values of $\omega_{n}$
and $k_{n}(x_{n})$ are most likely due to imperfections, such as the presence
of the gold layer and the undercuts beneath the anchors. There is also an
estimated 10% error in $F_{T}$ due to the fact that we do not know the exact
value of $E$ [28]. Thus, it is more justifiable to use the $\omega_{n}$ and
$k_{n}(x_{n})$ values directly obtained from experiments rather than those
calculated from elasticity theory.

###  IV.4 Water

Our beams are then immersed in water, where we measure the PSDs of
displacement fluctuations, $G_{W}(x_{0},\omega)$, over the continuous
frequency range of $50~{}{\rm kHz}$ to $15~{}\rm MHz$ at two positions on the
beam, $x_{0}=0.25$ and $x_{0}=0.50$. Fig. 6 shows $G_{W}(x_{0},\omega)$ as a
function of frequency at these two positions for all three beams. The first
position, $x_{0}=0.50$, corresponds to an antinode of all odd modes and a node
of all even modes [Fig. 6(a-c)]; the second position, $x_{0}=0.25$, is close
to the antinode of the second mode [Fig. 6(d-f)]. Low quality factors
($Q\gtrsim 1$) in water result in broad and overlapping peaks; the peak
frequencies are significantly lower than those in air.

Next, we compare our measurements in water to our theoretical expression for
$G_{W}(x_{0},\omega)$. The continuous curves in Fig. 6 show predictions based
on Eq. (18). Here, we use the $\omega_{n}$ values found above and determine
$k_{n}(x_{0})$ from $k_{n}(x_{n})$ after correcting for the position
dependence via Eq. (16). We then compute $T_{0}$ and $\Gamma(\omega)$ using
the density and viscosity of water, and combine all factors in Eq. (18) to
generate the curves. The dotted curves show the PSDs of individual modes; the
continuous curve is the sum of the first $12$ modes. A strong agreement
between experiment and theory is evident for $n\lesssim 7$; for $n\gtrsim 7$,
the beam fluctuations remain below our resolution limit. The theory predicts
the peak frequencies and the noise power levels accurately.

The positioning error mentioned in the third paragraph in Section III affects
both $k_{n}(x_{0})$ and the measured $G_{W}(x_{0},\omega)$. This error becomes
more pronounced at higher frequencies. In principle, the theory curves can be
further improved by treating the measurement position as a fit parameter.

##  V Discussion

Eq. (18) describes the Brownian dynamics of a nanomechanical beam in a viscous
fluid and indicates that, at a given frequency, the total noise is found by
adding the noise PSDs in different eigenmodes. The underlying assumption is
that the Brownian force noise is delta-function correlated in space [47, 33].
The form of Eq. (18), i.e., the summation over uncorrelated eigenmodes, should
remain unchanged for any mechanical system as long as the dissipation is
uniform in space. For our system, the dissipation in Eq. (5) from the cylinder
model is indeed spatially homogeneous. The remarkable agreement between
experiment and theory in Fig. 6 for the first seven modes suggests that the
cylinder model remains accurate—to within our experimental resolution. In
other words, the frequency dependence and the spatial homogeneity of the
dissipation in the model are both validated by our experiments.

Hydrodynamic fluctuations in a simple fluid are typically assumed to be delta-
function correlated in space [48]. However, the situation is different for the
nanomechanical beam immersed in a fluid: the flow around the structure and the
fluid-structure interactions are expected to result in spatial correlations in
the force noise, eventually leading to observable deviations from Eq. (18) for
higher modes. This expectation is consistent with the fact that the viscous
dissipation of the oscillating cylinder model is not accurate for higher
modes. As the flow in the axial direction becomes more appreciable with
increasing mode number, the dissipation becomes mode dependent and non-uniform
[49, 39]. However, the agreement between predictions and experiments suggests
that the axial flow is negligible in our parameter space. The smallest length
scale probed in our beams is comparable to the smallest resolved modal
wavelength of $\sim 9~{}\rm\mu m$.

$G_{W}(x_{0},\omega)$ should also be affected by the presence of the nearby
substrate due to the squeeze flow between the beam and the substrate. In our
theory, we have neglected the presence of the substrate. This can be corrected
using numerical simulations [49]: we estimate a decrease in
$\Gamma^{\prime}(\omega/2\pi=0.1~{}\rm MHz)$ by $\sim 30\%$ and an increase in
$\Gamma^{\prime\prime}(\omega/2\pi\\!=\\!0.1~{}\rm MHz)$ by $\sim 30\%$, which
leads to an increase in $G_{W}(\omega/2\pi\\!=\\!0.1~{}\mathrm{MHz},x_{0})$ of
$\lesssim 30\%$ depending on $L$ of the beam. This additional damping should
further broaden the fundamental mode and induce an additional redshift of the
peak [50, 25], which we observe in measurements. At higher frequencies, as the
viscous boundary layer thickness $\delta$ becomes $\delta/g\\!\ll\\!1$, this
effect disappears [51, 27].

To experimentally observe the noise correlation effects due to fluid-structure
interaction, it would be necessary to resolve the fluctuations in higher
modes. To this end, one should first determine the viscous dissipation from
fluid dynamics and asses the regime where the dissipation becomes
nonhomogenous, e.g., due to axial flows. Unfortunately, analytical solutions
cannot be found for most experimental geometries, making numerical simulations
necessary. Once the flow regimes around the structure are determined, one
could design and fabricate structures with softer spring constants $k_{n}$
such that the Brownian motions in higher modes can be resolved.

Nonlinearities in the elastic potential can modify the Brownian dynamics of
the nanomechanical beam [52, 53]. Achieving the nonlinear limit for thermal
fluctuations of a mode requires a very high modal $Q_{n}$ and a very soft
modal spring [53]. We estimate that the very low $Q_{n}$ in fluids along with
the large $k_{n}$ of the beams makes nonlinear effects negligible in our
experiments. As an order of magnitude comparison, the fundamental mode of our
50-$\rm\mu m$-long beam displays nonlinear behavior at amplitudes $\gtrsim
20~{}\rm nm$ [54], while its thermal amplitude remains around $60~{}\rm pm$ in
air. The emergence of nonlinear behavior in higher modes of elastic structures
and Brownian motion in a nonlinear potential are both interesting questions
requiring further study.

###### Acknowledgements.

We acknowledge support from US NSF through grants CMMI-2001559, CMMI-1934271,
CMMI-1934370, and CMMI-2001403. …………. …………. ………….. ………. ……….. ……… ……….

## References

  * Uhlenbeck and Ornstein [1930] G. E. Uhlenbeck and L. S. Ornstein, On the theory of the Brownian motion, Phys. Rev. 36, 823 (1930).
  * Kappler [1931] E. Kappler, Versuche zur Messung der Avogadro-Loschmidtschen Zahl aus der Brownschen Bewegung einer Drehwaage, Ann. Phys. (Leipzig) 403, 233 (1931).
  * Parson [1915] A. Parson, A highly sensitive electrometer, Phys. Rev. 6, 390 (1915).
  * Moll and Burger [1925] W. Moll and H. Burger, LXV. The sensitivity of a galvanometer and its amplification, Lond. Edinb. Dublin Philos. Mag. J. Sci. 50, 626 (1925).
  * Barnes and Silverman [1934] R. B. Barnes and S. Silverman, Brownian motion as a natural limit to all measuring processes, Rev. Mod. Phys. 6, 162 (1934).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, Mesoscopic physics of nanomechanical systems, Rev. Mod. Phys. 94, 045005 (2022).
  * Butt and Jaschke [1995] H.-J. Butt and M. Jaschke, Calculation of thermal noise in atomic force microscopy, Nanotechnology 6, 1 (1995).
  * Saulson [1990] P. R. Saulson, Thermal noise in mechanical experiments, Phys. Rev. D 42, 2437 (1990).
  * Cleland [2013] A. N. Cleland,  _Foundations of Nanomechanics_ (Springer, Berlin, 2013).
  * Cleland and Roukes [2002] A. N. Cleland and M. L. Roukes, Noise processes in nanomechanical resonators, J. Appl. Phys. 92, 2758 (2002).
  * Paolino  _et al._ [2009] P. Paolino, B. Tiribilli, and L. Bellon, Direct measurement of spatial modes of a microcantilever from thermal noise, J. Appl. Phys. 106, 094313 (2009).
  * Wang  _et al._ [2014] Z. Wang, J. Lee, and P. X.-L. Feng, Spatial mapping of multimode Brownian motions in high-frequency silicon carbide microdisk resonators, Nat. Commun. 5, 1 (2014).
  * McRae  _et al._ [2010] T. G. McRae, K. H. Lee, G. I. Harris, J. Knittel, and W. P. Bowen, Cavity optoelectromechanical system combining strong electrical actuation with ultrasensitive transduction, Phys. Rev. A 82, 023825 (2010).
  * Gloppe  _et al._ [2014] A. Gloppe, P. Verlot, E. Dupont-Ferrier, A. Siria, P. Poncharal, G. Bachelier, P. Vincent, and O. Arcizet, Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field, Nat. Nanotechnol. 9, 920 (2014).
  * Arcizet  _et al._ [2006] O. Arcizet, P.-F. Cohadon, T. Briant, M. Pinard, A. Heidmann, J.-M. Mackowski, C. Michel, L. Pinard, O. Français, and L. Rousseau, High-sensitivity optical monitoring of a micromechanical resonator with a quantum-limited optomechanical sensor, Phys. Rev. Lett. 97, 133601 (2006).
  * Levin [1998] Y. Levin, Internal thermal noise in the ligo test masses: A direct approach, Phys. Rev. D 57, 659 (1998).
  * Liu and Thorne [2000] Y. T. Liu and K. S. Thorne, Thermoelastic noise and homogeneous thermal noise in finite sized gravitational-wave test masses, Phys. Rev. D 62, 122002 (2000).
  * Yamamoto  _et al._ [2001] K. Yamamoto, S. Otsuka, M. Ando, K. Kawabe, and K. Tsubono, Experimental study of thermal noise caused by an inhomogeneously distributed loss, Phys. Lett. A 280, 289 (2001).
  * Yamamoto  _et al._ [2002] K. Yamamoto, S. Otsuka, M. Ando, K. Kawabe, and K. Tsubono, Study of the thermal noise caused by inhomogeneously distributed loss, Class. Quantum Gravity 19, 1689 (2002).
  * Schwarz  _et al._ [2016] C. Schwarz, B. Pigeau, L. M. De Lépinay, A. G. Kuhn, D. Kalita, N. Bendiab, L. Marty, V. Bouchiat, and O. Arcizet, Deviation from the normal mode expansion in a coupled graphene-nanomechanical system, Phys. Rev. Appl. 6, 064021 (2016).
  * Franosch  _et al._ [2011] T. Franosch, M. Grimm, M. Belushkin, F. M. Mor, G. Foffi, L. Forró, and S. Jeney, Resonances arising from hydrodynamic memory in Brownian motion, Nature 478, 85 (2011).
  * Paul and Cross [2004] M. R. Paul and M. C. Cross, Stochastic dynamics of nanoscale mechanical oscillators immersed in a viscous fluid, Phys. Rev. Lett. 92, 235501 (2004).
  * Paul  _et al._ [2006] M. R. Paul, M. T. Clark, and M. C. Cross, The stochastic dynamics of micron and nanoscale elastic cantilevers in fluid: fluctuations from dissipation, Nanotechnology 17, 4502 (2006).
  * Clark  _et al._ [2010] M. T. Clark, J. E. Sader, J. P. Cleveland, and M. R. Paul, Spectral properties of microcantilevers in viscous fluid, Phys. Rev. E 81, 046306 (2010).
  * Clarke  _et al._ [2006] R. J. Clarke, O. E. Jensen, J. Billingham, A. P. Pearson, and P. M. Williams, Stochastic elastohydrodynamics of a microcantilever oscillating near a wall, Phys. Rev. Lett. 96, 050801 (2006).
  * Bellon [2008] L. Bellon, Thermal noise of microcantilevers in viscous fluids, J. Appl. Phys. 104, 104906 (2008).
  * Pierro  _et al._ [2019] E. Pierro, F. Bottiglione, and G. Carbone, Thermal fluctuations and dynamic modeling of a dAFM cantilever, Adv. Theory Simul. 2, 1900004 (2019).
  * Ari  _et al._ [2020] A. B. Ari, M. S. Hanay, M. R. Paul, and K. L. Ekinci, Nanomechanical measurement of the Brownian force noise in a viscous liquid, Nano Lett. 21, 375 (2020).
  * Ti  _et al._ [2021] C. Ti, A. B. Ari, M. C. Karakan, C. Yanik, I. I. Kaya, M. S. Hanay, O. Svitelskiy, M. González, H. Seren, and K. L. Ekinci, Frequency-dependent piezoresistive effect in top-down fabricated gold nanoresistors, Nano Lett. 21, 6533 (2021).
  * Barbish  _et al._ [2022] J. Barbish, C. Ti, K. L. Ekinci, and M. R. Paul, The dynamics of an externally driven nanoscale beam that is under high tension and immersed in a viscous fluid, J. Appl. Phys. 132, 034501 (2022).
  * Schlichting [1932] H. Schlichting, Berechnung ebener periodischer Grenzschichtströmungen, Phys. Z. 33, 327 (1932).
  * Retsina  _et al._ [1987] T. Retsina, S. Richardson, and W. Wakeham, The theory of a vibrating-rod viscometer, Appl. Sci. Res. 43, 325 (1987).
  * Sader [1998] J. E. Sader, Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope, J. Appl. Phys. 84, 64 (1998).
  * Rosenhead [1963] L. Rosenhead,  _Laminar Boundary Layers_ (Clarendon, Oxford, 1963).
  * Bokaian [1990] A. Bokaian, Natural frequencies of beams under tensile axial loads, J. Sound Vib. 142, 481 (1990).
  * Stachiv [2014] I. Stachiv, Impact of surface and residual stresses and electro-/magnetostatic axial loading on the suspended nanomechanical based mass sensors: A theoretical study, J. Appl. Phys. 115, 214310 (2014).
  * Callen and Welton [1951] H. B. Callen and T. A. Welton, Irreversibility and generalized noise, Phys. Rev. 83, 34 (1951).
  * Callen and Greene [1952] H. B. Callen and R. F. Greene, On a theorem of irreversible thermodynamics, Phys. Rev. 86, 702 (1952).
  * Van Eysden and Sader [2007] C. Van Eysden and J. E. Sader, Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope: Arbitrary mode number, J. Appl. Phys. 101, 044908 (2007).
  * Wagner [1990] J. W. Wagner, Optical detection of ultrasound, Phys. Acoust. 19, 201 (1990).
  * Kara  _et al._ [2015] V. Kara, Y.-I. Sohn, H. Atikian, V. Yakhot, M. Loncar, and K. L. Ekinci, Nanofluidics of single-crystal diamond nanomechanical resonators, Nano Lett. 15, 8070 (2015).
  * Kouh  _et al._ [2005] T. Kouh, D. Karabacak, D. Kim, and K. Ekinci, Diffraction effects in optical interferometric displacement detection in nanoelectromechanical systems, Appl. Phys. Lett. 86, 013106 (2005).
  * Kara  _et al._ [2017] V. Kara, V. Yakhot, and K. L. Ekinci, Generalized Knudsen number for unsteady fluid flow, Phys. Rev. Lett. 118, 074505 (2017).
  * Tai and Muller [1990] Y.-C. Tai and R. Muller, Measurement of Young’s modulus on microfabricated structures using a surface profiler, in  _Proceedings on Micro Electro Mechanical Systems, An Investigation of Micro Structures, Sensors, Actuators, Machines and Robots._ (IEEE, Napa Valley, 1990) pp. 147–152.
  * Zhang  _et al._ [2000] T.-Y. Zhang, Y.-J. Su, C.-F. Qian, M.-H. Zhao, and L.-Q. Chen, Microbridge testing of silicon nitride thin films deposited on silicon wafers, Acta Mater. 48, 2843 (2000).
  * Kuhn  _et al._ [2000] J. L. Kuhn, R. K. Fettig, S. H. Moseley Jr, A. S. Kutyrev, and J. Orloff, Fracture tests of etched components using a focused ion beam machine, in  _MEMS Reliability for Critical Applications_ , Vol. 4180 (SPIE, Santa Clara, 2000) pp. 40–48.
  * Newland [2012] D. E. Newland,  _An Introduction to Random Vibrations, Spectral & Wavelet Analysis_ (Dover Publications, New York, 2012).
  * Landau  _et al._ [1980] L. D. Landau, E. Lifshitz, and L. Pitaevskii,  _Course of Theoretical Physics: Statistical Physics, Part 2_ (Pergamon, Oxford, 1980).
  * Liem  _et al._ [2021] A. T. Liem, A. B. Ari, C. Ti, M. J. Cops, J. G. McDaniel, and K. L. Ekinci, Nanoflows induced by MEMS and NEMS: Limits of two-dimensional models, Phys. Rev. Fluids 6, 024201 (2021).
  * Green and Sader [2005] C. P. Green and J. E. Sader, Frequency response of cantilever beams immersed in viscous fluids near a solid surface with applications to the atomic force microscope, J. Appl. Phys. 98, 114913 (2005).
  * Clarke  _et al._ [2005] R. J. Clarke, S. M. Cox, P. Williams, and O. Jensen, The drag on a microcantilever oscillating near a wall, J. Fluid Mech. 545, 397 (2005).
  * Dykman and Krivoglaz [1971] M. Dykman and M. Krivoglaz, Classical theory of nonlinear oscillators interacting with a medium, Phys. Status Solidi (B) 48, 497 (1971).
  * Gieseler  _et al._ [2013] J. Gieseler, L. Novotny, and R. Quidant, Thermal nonlinearities in a nanomechanical oscillator, Nat. Phys. 9, 806 (2013).
  * Postma  _et al._ [2005] H. Postma, I. Kozinsky, A. Husain, and M. Roukes, Dynamic range of nanotube-and nanowire-based electromechanical systems, Appl. Phys. Lett. 86, 223105 (2005).


# Simple U-net Based Synthetic Polyp Image Generation: Polyp to Negative and
Negative to Polyp

Hemin Ali Qadir  Ilangko Balasingham  Younghak Shin  The Intervention Centre,
Oslo University Hospital (OUS), Sognsvannsveien 20, 0372 Oslo, Norway
Department of Electronic Systems, Norwegian University of Science and
Technology (NTNU), Høgskoleringen 1, 7491 Trondheim, Norway  Department of
Computer Engineering, Mokpo National University, Jeollanam-do, Muan-gun,
Cheonggye-myeon, Yeongsan-ro, 1666 KR, South Korea

###### Abstract

Synthetic polyp generation is a good alternative to overcome the privacy
problem of medical data and the lack of various polyp samples. In this study,
we propose a deep learning-based polyp image generation framework that
generates synthetic polyp images that are similar to real ones. We suggest a
framework that converts a given polyp image into a negative image (image
without a polyp) using a simple conditional GAN architecture and then converts
the negative image into a new-looking polyp image using the same network. In
addition, by using the controllable polyp masks, polyps with various
characteristics can be generated from one input condition. The generated polyp
images can be used directly as training images for polyp detection and
segmentation without additional labeling. To quantitatively assess the quality
of generated synthetic polyps, we use public polyp image and video datasets
combined with the generated synthetic images to examine the performance
improvement of several detection and segmentation models. Experimental results
show that we obtain performance gains when the generated polyp images are
added to the training set.

###### keywords:

Colonoscopy , polyp detection , image synthesis , convolutional neural network
, generative adversarial networks

††journal: Biomedical Signal Processing and Control

##  1 Introduction

Colorectal cancer (CRC) is the third most frequent cause of cancer mortality
and the second leading cause of death in men and women globally [1]. Polyps
are precursors to colorectal cancer, and resection of polyps can effectively
prevent the development of colorectal cancer [2]. Colonoscopy is an important
standard tool for colon screening globally. It has the advantage of being able
to visualize and remove lesions simultaneously during screening. However, the
effectiveness of a colonoscopy depends on the skills and expertise of the
endoscopy specialist. Therefore, the polyp detection rate varies widely, and
the polyp miss rate is known to be between 6 and 27% in clinical trials [3].

Over the past two decades, many researchers have been working on automatic
polyp detection using efficient computer-aided detection systems. Early
studies attempted to detect polyp by extracting features such as color, shape,
and texture of polyp [4, 5, 6]. Recently, attempts have been made to
automatically detect and segment polyps using deep learning models which
improved performance over traditional feature extraction-based methods [7, 8,
9, 10, 11, 12].

Most deep learning-based studies use several public polyp datasets [7, 9, 10,
11, 12] for research purposes and focus on developing deep learning model
architectures with proper pre- and post-processing methods to improve polyp
detection and segmentation performance. One of the difficulties in polyp
detection and segmentation is the diversity of their features. Polyps vary
widely in color, shape, size, and texture. However, widely used polyp datasets
only have dozens of unique polyps used for deep learning-based detection and
segmentation model training [13, 14, 15]. Because in medical fields, abnormal
data is inherently more difficult to collect than normal data. To obtain the
location information required for polyp detection and segmentation, manual
annotations by colonoscopy experts are essential and are usually costly and
logistically difficult.

In order to overcome the limitation of availability of large polyp image data,
many deep learning-based studies are actively using transfer learning and
image augmentation techniques [9, 5]. However, to effectively train a deep
learning model with a large number of learnable parameters, it is necessary to
increase the diversity of polyp appearances within the training dataset.

Recently, a GAN(generative adversarial network)-based image generation
approach has also been applied to polyp data to overcome insufficient training
datasets [16, 17]. In [16], GAN based domain adaptation approach has been
applied to the Endoscopy dataset. They focused on adapting real data to
synthetic-like data to remove patient-specific details from real images. This
is the difference from our study, which aims to generate new synthetic polyp
images to improve polyp detection and segmentation performance.

In [18], to generate a synthetic polyp, a polyp mask is first generated using
a GAN framework, and the polyp part is generated when the polyp mask is given
using the conditional GAN approach. Then, the polyp-free part is found in the
real polyp image and the generated polyp is synthesized. The difference from
the proposed approach is that two GAN frameworks should be used for polyp
generation. In addition, a separate algorithm for each image is required to
find the polyp-free part, and another separate algorithm is used to naturally
synthesize the generated polyp part into the existing image. On the other
hand, the proposed method uses a single conditional GAN framework to convert a
polyp image into a negative image and convert the negative image back to a
polyp image. Therefore, there is no need to use a separate algorithm to find a
polyp-free part or perform natural synthesis.

In [17], a conditional GAN was used to generate synthetic polyps. The authors
proposed to use a conditional mask to control the position, size, and shape of
generated polyps and showed improved polyps detection performance when the
generated polyps were added as additional training data. However, negative
images must be used to generate the polyp image, and additional processing
(i.e., edge filtering) was required to create the input condition image.
Furthermore, it is not possible to generate polyps with various
characteristics from the same input condition. From a performance verification
point of view, the performance improvement was weakly evaluated, i.e., only
one deep learning model was evaluated for the polyp detection task.

In this paper, to overcome the shortcomings of the previous polyp image
generation work, we conduct a study that generates synthetic polyp images
similar to real ones using only a given polyp dataset without the use of
additional datasets and processes. We propose a framework that uses the same
GAN architectures to convert the polyp images into a negative image and the
negative image back into a novel polyp image. In addition, we attempt to
generate various polyp features by introducing a controllable polyp mask with
one input condition image. The experimental results show performance
improvement when synthetic polyp images generated by the proposed method are
used for training various deep learning models for object detection and
segmentation.

The remainder of this paper is organized as follows. In Section 2, the
proposed polyp to negative and negative to polyp image generation framework
including network architecture and preparation of input conditioned images are
introduced. We also explain the polyp detection and segmentation models used
for the qualitative evaluation of generated polyp images. In Section 3,
experimental datasets used in this study are described. In Section 4,
experimental results and discussions are described. Finally, we conclude this
study in Section 5.

Fig. 1: Proposed conditional GANs based polyp image generation framework. Fig.
2: Proposed U-Net based generator and discriminator architectures. The same
architecture is used for both polyp to negative and negative to polyp tasks.

##  2 Methods

Figure 1 shows the concept of the proposed synthetic polyp generation
procedure. Using a given polyp image and its corresponding label mask, we
first transform the polyp part to the negative part. Then, we transform the
generated negative image into a new synthetic polyp image using a controllable
input condition mask. The controllable mask is used to force the generator
network to produce polyps with different characteristics depending on the
color of the polyp mask in the input condition image (see Section 2.3). For
both transformations, we use the same network architectures. We explain each
method in the following subsections in detail.

###  2.1 Network Architectures

Figure 2 shows the structure of the generator and discriminator networks used
for the conditional GAN framework in this study. Recent GAN studies have
focused on proposing new network architectures to obtain task-specific high
performance image generation [17], image inpainting [19] and image detection
and segmentation [12, 20]. On the other hand, we adopt a simple U-net based
generator architecture used in [21] to make it universally work for two
different tasks (e.g., converting from polyp to negative and from negative to
polyp). The difference from the paper [21] is that a multi-patch discriminator
is used in the discriminator and WGAN-GP loss [22] is adopted as the GAN loss
function. Detailed explanations are provided in the following paragraphs.

As shown in Figure 2, for the input image size of 256×256, 4×4 convolution
with stride 2 is applied in each encoding layer, so that the feature map size
is reduced to 1/2 and then progresses to the final 1×1 feature map. In each
subsequent decoding layer, transposed convolution is applied to increase the
feature map size by 2 and a final 256×256 image is generated. In the encoder
part of the generator, we use Convolution-BatchNorm-leakyReLU module in each
layer while in the decoder part, we use transposed Convolution-BatchNorm-ReLU
module. Dropout is used in the first 3 layers of the decoder part with a
dropout rate of 50%.

In the case of the discriminator, a patch-based discriminator which is widely
used for image-to-image transformation tasks is applied [17, 21].
Additionally, to reflect the diversity of polyp sizes well, a multi-patch
discriminator was used that utilizes two output patches at the same time as
shown in Figure 2. In the discriminator, the Convolution-BatchNorm-leakyReLU
module was applied to the input image in the same way as the encoding layer of
the generator, and the final feature map size is 16×16 as shown in Figure 2.
For the multi-patch discriminator, GAN loss for 64×64 and 16×16 feature maps
is added.

Our training loss is comprised of a GAN term and a reconstruction (Reconst)
term which are commonly used in the existing methods [17, 21]:

$arg\>\underset{G}{min}\>\underset{D}{max}\>L_{GAN}(G.D)+\lambda
L_{Reconst}(G).$ (1)

where $G$ refers to the generator network, $D$ refers to the discriminator
network, and $\lambda$ is a parameter to control the balance between two
different loss terms. In the case of the GAN loss term, we use the WGAN-GP
loss [22], which is known to improve training stability and has been recently
used in many GAN studies [22, 23]. For reconstruction loss, we use the L1 loss
between the generated image and the original target image [21]. L1 loss is
needed to model low-frequency structure (blurry results) because the GAN loss
can only model high-frequency structure and give much sharper results. The
framework can introduce virtual artifacts when the GAN loss is used alone
[21]. When both losses are used together, they lead the framework to generate
realistic images to fool $D$.

###  2.2 Polyp Image to Negative Image Generation

To generate a new polyp image, we first train our conditional GAN model to
transform a given polyp image to a new negative image. To do this, we need to
prepare a pair of training images, i.e., input conditioned image and output
target image as shown in Figure 1.

Each row in Figure 3 is an example of a pair of images used for training polyp
to negative image transformation. An input condition image (shown on the left
side in Figure 3) is created from combining a polyp training image and a
generated mask. The original polyp image (shown on the right side in Figure 3)
is used as an output target image. Thus, the proposed conditional GAN trains
an image inpainting task that targets generating the negative parts in the
right images by using the masks in the left images in Figure 3 as its
condition.

To create the masks shown in the left images in Figure 3, we apply basic image
augmentation such as rotation, scaling, position translation, and perspective
transform with randomly selected parameters to the original polyp masks in the
training dataset. Since the network will be trained to reconstruct the
negative image portion only, the generated masks are randomly placed in the
input condition images, avoiding overlapping with the polyp portion.

Fig. 3: Training image pairs for polyp to negative image generation. Left:
input condition image, Right: output target image.

In the inference phase, the goal is to transform the polyp region in a given
polyp image to a negative region. The polyp image is combined with its
corresponding polyp mask to create the input condition image. However, if we
use the original polyp mask, shadows or connections around the polyp will
remain as shown in Figure 4(b). This shadow will appear as an artifact and
prevent generating a realistic negative image. To avoid this issue, we dilate
the original polyp mask by 10 pixels (as shown in Figure 4(c),) as the final
mask for the input condition in the inference.

###  2.3 Negative Image to Polyp Image Generation

Next, we prepare training data for the final objective, the negative image to
polyp image transformation task. Figure 5 shows two pairs of training images
used for the negative image to polyp image transformation. The input condition
image is generated from combining the polyp images with their corresponding
polyp masks provided with the datasets. Figure 5 shows two pairs of training
images: the input condition image is shown on the left side, and the output
target shown on the right is the polyp image. We aim to train the generator
network to regenerate the same form of polyp in the given mask region.

In the previous image to image translation method [21], there is a shortcoming
in that the same image is always generated for the same input condition. Also,
in [17], there is a shortcoming in that the characteristics of the generated
polyps for various input conditions are very similar. The reason is that a
single binary polyp mask was used for various polyps. To overcome this
problem, we try to use a controllable input condition. It is intended to
generate a new polyp by adjusting the mask values even if the the shape of the
mask in the input condition is the same shape (see the left side in Figure 5).

Fig. 4: Example of input conditioned image for negative image generation in
inference. (a) original polyp image (b) original polyp image combined with
original polyp mask (c) original polyp image combined with extended polyp
mask.

The dataset of polyp images used for training comprises 34 different unique
polyps (see the information of CVC-ClinicDB dataset in section 3). We allocate
34 different values evenly distributed between 0 and 255 rather than
traditional binary masks for the polyp regions in the training images.
Therefore, the value of the polyp masks varies depending on the polyp in the
given image as shown on the left side of Figure 5, meaning a different gray-
scale value is assigned for each unique polyp in the dataset.

Fig. 5: Training image pairs for negative to polyps image generation. Left:
input condition images, Right: output target images (original polyp images).

By controlling the value of the polyp masks in the input condition images, we
can generate a polyp with different features and characteristics from the same
input condition of the same shape at the inference phase, see Section 4 for
detailed generation results. When inferring the final polyp generation, the
input condition image is produced by combining a label mask, which is created
from a randomly generated shape and selected value, with a negative image,
which is obtained by the polyp to negative image transformation model.

###  2.4 Polyp Detection and Segmentation Models

For quantitative evaluation of the generated polyp images, we assess the
performance of several polyp detection and segmentation models when trained
with and without generated synthetic polyps. In more detail, performance
evaluation is performed on several detection and segmentation models trained
only with the data given as the original training set and the models trained
with the generated polyp images added to the original training set.

For polyp detection, we select three well-known detection models: Faster R-CNN
ResNet101, Faster R-CNN Inception Resnet v2 [24] and R-FCN ResNet101 [25].
Faster R-CNN utilizes a pre-trained deep CNN as a feature extractor and uses
the last convolution layer as a region proposal network (RPN) [24]. Then, the
classification and regression are performed on the RoIs (Region of Interest)
proposed by the RPN. We use two deep CNN models i.e. ResNet101 and Inception
Resnet v2, as the feature extractors pre-trained on Microsoft’s (MS) Common
Objects in Context (COCO) datasets.

R-FCN has been proposed to effectively overcome the slow training and
inference speeds occurring during RoI-specific processing of the existing
Faster R-CNN [25]. R-FCN implements the position-sensitive score maps
technique that can include information about the location of objects while
sharing the operation of RoIs extracted from RPN using a fully convolutional
layer. Resnet101 pre-trained on MS COCO dataset is used as the feature
extractor of R-FCN.

We also use three segmentation models to qualitatively evaluate the generated
synthetic polyps when used as additional training data. All three models are
based on encoder-decoder architecture of the U-Net family: TernausNet-16 [20,
26], AlbuNet-34 [20], and MDeNetplus [12]. U-Net is a widely used segmentation
model in the field of medical imaging where training data is limited [27]. A
U-Net architecture consists of an encoder-decoder structure with skip
connections. The skip connections enable feature combinations of the encoder
layers and decoder layers, which allows for sophisticated localization. The
TernausNet-16 model uses ImageNet pre-trained VGG-16 for the encoder network,
while the AlbuNet-34 uses ImageNet pre-trained ResNet-34 as an encoder and
improves the method of skip-connections from U-Net [20]. MDeNetplus has not
only skip connections from the encoder layers to decoder layers, but also
feedback connections. The feedback connections sum the activation maps of
similar layers of different decoders.

###  2.5 Model Training

We train the GAN network from scratch and the Adam optimizer with a momentum
of 0.5 at a learning rate of 0.0002. Before training, the input images of
256×256 are resized to 312×312 and then randomly cropped back to 256×256 for
applying random jittering [21]. For other parameters in model training such as
learning rate, batch size, and weight for reconstruction loss, we refer to the
values which were used in the image to image translation work [21].

For training of detection models, i.e., Faster R-CNN ResNet101 and Inception
Resnet v2 and R-FCN ResNet101, we use the Tensorflow object detection API
[28]. The stochastic gradient descent (SGD) method with a momentum of 0.9 is
used as an optimizer. In each iteration of the RPN training, 256 training
samples are randomly selected from each training image and the ratio between
positive (‘polyp’) and negative (‘background’) samples is 1:1. For all other
parameters such as learning rate, non-maximum suppression (NMS), and the
maximum number of proposals, we use the default values which were used in the
original Faster R-CNN work [24]. Each segmentation model is trained with an
Adam optimizer with a learning rate of 0.0001 and batch size of 5 for a
different number of epochs to achieve the maximum performance for each model.
During training, the training dataset was randomly split into training and
validation subsets using 5-fold cross-validation. The 5-fold cross-validation
method was performed to choose the best hyper-parameters for the models (e.g.
the learning rate relative to the batch size, best epoch, etc) and thus avoid
overfitting. The size of the images was changed to 512×512. Finally, we use
Jaccard loss combined with a pixel-wise binary cross-entropy to optimize the
parameters of the models

Fig. 6: Some results of the polyp to negative translation. Each column
represents different image generation result. (a) is polyp image used for
preparing the corresponding input conditioned image (b). (c) is generated
negative image from the corresponding input condition (b).

##  3 Experimental Datasets

We use three publicly available polyp datasets, CVC-ClinicDB [13], ETIS-LARIB
[14] and CVC-ClinicVideoDB [15]. The CVC-ClinicDB image dataset is used for
the whole procedure of the polyp image generation task, i.e., model training
for negative and polyp image generation. It is also utilized for training
detection and segmentation models as a baseline original image set which is
compared to those of the original plus generated synthetic image set. The CVC-
ClinicVideoDB and ETIS-LARIB image dataset are used for detection and
segmentation model tests, respectively.

The CVC-ClinicDB image dataset consists of 612 polyp image frames with
corresponding 612 binary polyp mask images. The mask images were annotated by
skilled clinicians. The 612 polyp image frames have a pixel resolution of
388×284 pixels in SD (standard definition). These images are extracted from 31
colonoscopy videos which contain 34 unique polyps. The CVC-ClinicVideoDB
dataset contains 18 SD (384×288 pixels) videos. In this dataset, a total of
10025 frames out of 11954 frames contain a polyp. For each image in the video
database, an approximate circle-shaped binary polyp mask is given which makes
the dataset suitable for the evaluation of performance improvement the polyp
detection models and benchmark our results with the existing studies. For the
evaluation of detection models, we use the whole 11954 frames as a test set.

The ETIS-LARIB dataset comprises 196 still images extracted from 34
colonoscopy videos. The images have an HD (high definition) 1255×966 pixels
resolution. In this dataset, 44 different polyps are presented in 196 images.
A binary polyp segmentation mask annotated by an experienced clinician is
provided for each image. This dataset is suitable for the evaluation of
performance improvement of the segmentation models because the polyp masks are
polygon shapes drawn around the polyp boundaries.

##  4 Results and Discussion

###  4.1 Generated Negative Images

In this work, we first transform the given polyp image into a negative image
and then create a new polyp in the generated negative image using the polyp
condition. In both image translation tasks, i.e., polyp to negative and
negative to polyp, we use the same proposed conditional GAN architectures.

Figure 6 shows some results of the generated negative images from polyp to
negative translation task. In each column of Figure 6, a different generated
negative image is represented in (c) which is created from images presented in
(b) as the input condition image and (a) is the corresponding original polyp
image. Thus, from this task, we transform polyp image (a) to negative image
(c). As we can see in Figure 6, the generated negative images not only clearly
remove the corresponding polyp in the polyp portion given by input condition
(b), but also naturally harmonize the removed and surrounding parts. We can
see that it produces natural negative images even for the third and fifth
columns, which are relatively large polyps.

Figure 7 is an example of generated negative images with large texture changes
in the part of the binary condition mask. As we can see in the generated
negative images (c) of Figure 7, our trained model adaptively creates a
realistic negative colon image by applying the shape and texture of the
surrounding parts of the polyp mask. As a result, we can see that capillaries
or folds of colon in the generated negative images in the blue circle part of
(c).

Fig. 7: Example of images with large texture changes in the corresponding
polyp mask part among the generated negative images. (a), (b) and (c)
represents original polyp image, input condition and generated negative image.
Fig. 8: Some results of the negative to polyp translation. Each column
represents a different image generation result. (a) is generated negative
image from polyp to negative translation. (b) and (c) is the input conditioned
image and the corresponding generated polyp image.

###  4.2 Generated Polyp Images

Our main goal is to create new polyp-looking images. For this purpose, the
negative to polyp translation is performed using the previously mentioned
method. Figure 8 shows several examples of negative to polyp translation and
each column shows a different image result. In this experiment, we use the
generated negative image (a) for the preparation of the input condition image
(b). Then, our trained model, i.e., negative to polyp model generated output
polyp image (c) using the (b) as the input condition.

As mentioned in Section 2.3, in the input condition image (b) of Figure 8, the
value of the polyp mask is randomly selected among 34 different pixel values
representing different gray-scale colors. The resulting images in Figure 8 (c)
show that the position and shape of the generated polyp are controlled by the
polyp mask of (b). In addition, it can be seen that various gray-scale colors
and textures are generated rather than similar characteristics, which was a
challenge in the existing synthetic polyp generation study [17]. This is
because different pixel values are used as input conditions for the polyp mask
part.

To analyze this in detail, we have generated synthetic polyp images by varying
the pixel value of the polyp mask part in the same input condition image in
Figure 9. Each row in Figure 9 is an example of a polyp image generation
result when only the pixel value of the polyp mask part varies for the same
input condition. In other words, each column is generated from different input
conditions with different polyp mask values.

As it can be seen from the results of each row in Figure 9, even with the same
shape of polyp mask and background image, it is possible to create polyps with
various characteristics in terms of color or texture by controlling the pixel
value of the polyp mask. If we look at each column in Figure 9, we can confirm
that the model learns to generate synthetic polyps depending on the polyp mask
value rather than the shape or the background image. Because even if different
background images and shapes are used, the model can generate polyps with
similar characteristics. These two experiments confirm that diverse image
generation using the same input, a limitation of the existing image to image
transform study [17, 21], can easily be achieved by adjusting the mask pixel
value of the input condition rather than changing the network structure or
training loss. This methodology makes it possible to control the generation of
synthetic polyps with characteristics that are difficult for the model to
detect, or polyps with clinically malignant characteristics that should not be
overlooked.

Furthermore, we measured the inference speed of the generator network after
being trained within the Conditional GAN framework used in this paper. As a
result, the average inference time was measured to be 51.33 msec on an NVIDIA
RTX 2080Ti GPU.

###  4.3 Evaluation of Polyp Detection Performance

The advantage of synthetic polyp generation based on the proposed input
condition mask is that the generated polyp image can be used directly as a
training image for polyp detection or segmentation task without mask
annotations. This is significant because it can reduce the manual labeling
cost of experts in mask annotation of medical data.

For evaluation of generated polyp images, we first assess polyp detection
performance when generated polyp images are used as additional training
samples. As polyp detection models, the Faster R-CNN ResNet101, Faster R-CNN
Inception ResNet v2, and R-FCN ResNet101 models mentioned in section 2.4 are
used. To define evaluation metrics for polyp detection, it is important to
specify and compute the following medical terminologies: true positive (TP),
false positive (FP), false negative (FN), and true negative (TN) where:

TP = detection output within the polyp ground truth.

FP = any detection output outside the polyp ground truth.

FN = polyp not detected for positive (with polyp) image.

TN = no detection output for negative (without polyp) image.

If there is more than one detection output, only one TP is counted per polyp.
Based on the above parameters, we define three performance evaluation metrics,
precision (Pre), recall (Rec) and, f1-score (f1):

$Pre=\frac{TP}{TP+FP}\times 100.\quad Rec=\frac{TP}{TP+FN}\times 100.\\\
f1=\frac{2\times{Rec\times Pre}}{Rec\times Pre}\times
100.\quad\quad\quad\quad\quad\quad\quad\quad\quad$ (2)

For evaluation of polyp detection performance, the CVC-ClinicVideoDB dataset
(see Section 3 for details) is used as the test dataset.

Table 1 lists the polyp detection performance for three different models. The
original columns show the performance when the models are trained with The
CVC-ClinicDB image dataset (612 polyp images), which is the originally given
training dataset. The paper [17] columns show the detection results when the
models are trained on the original training dataset plus additional 372
synthetic polyp images generated by the method proposed in paper [17]. The
combined columns show the performance when the models are trained by adding
350 polyp images generated by the proposed method in this paper to the
original training dataset.

Table 1: Comparison of polyp detection performance for the three detection
models using the original training set and combined training set by generated
polyp image. Number of samples is: Original (612); [17] (612+372); Combined
(612+350)  Faster R-CNN ResNet101
---
Models | TP | TN | FP | FN | Rec | Pre | F1
Original | 6047 | 1431 | 1513 | 3978 | 60.32 | 79.99 | 68.76
[17] | 5370 | 1603 | 1049 | 4655 | 53.57 | 83.66 | 65.31
Combined | 6263 | 1508 | 991 | 3762 | 62.47 | 86.34 | 72.49
Faster R-CNN Inception ResNet v2
Models | TP | TN | FP | FN | Rec | Pre | F1
Original | 6011 | 1496 | 1333 | 4014 | 60 | 81.9 | 69.22
[17] | 6831 | 1399 | 1177 | 3194 | 68.1 | 85.3 | 75.74
Combined | 7056 | 1351 | 1212 | 2969 | 70.38 | 85.34 | 77.14
R-FCN ResNet101
Models | TP | TN | FP | FN | Rec | Pre | F1
Original | 5762 | 1304 | 2062 | 4263 | 57.48 | 73.65 | 64.56
[17] | 5554 | 1653 | 809 | 4471 | 55.4 | 87.29 | 67.78
Combined | 6555 | 1596 | 1032 | 3470 | 65.38 | 86.39 | 74.43

The results in Table 1 show that the use of generated images (Combined)
demonstrates better polyp detection performance for all three models than the
use of the original dataset only. For the three models, both recall and
precision increased. This means that polyps missed by the model trained with
the existing training dataset are improved when the combined dataset is used.
For example, when training using the combined set based on the Faster R-CNN
Inception ResNet v2 model, 1045 additional polyps are detected compared to the
polyps detected by the original set training. In addition, when the polyp
images generated by the proposed method are used for training, improved
performance is demonstrated in all three models based on f1-score compared to
using the polyp image generated by the existing method [17].

Fig. 9: Example of generated polyps. The images in each row are generated when
only the pixel value of the polyp mask is changed under the same input
condition..

In [18], the same dataset as of this study was used for the evaluation of
polyp detection performance, i.e., the detection model was trained on the
synthetic polyp generated by their proposed method together with the CVC-
ClinicDB training data. They evaluated the detection performance of Faster
R-CNN ResNet50 on the CVC-ClinicVideoDB. For comparison, we measure the
detection performance for the same model using a dataset generated by our
proposed method (i.e., Combined). Table 3 summarizes the detection performance
of Faster R-CNN ResNet50 for both methods and shows that our method obtains
better results in terms of Precision and F1-score while their model obtained
better Recall.

Table 2: Comparison of polyp detection performance of Faster R-CNN ResNet50 between our results and the results shown in [18] Models | Faster R-CNN ResNet50 |
---|---|---
TP | TN | FP | FN | Rec | Pre | F1
[18] | 7378 | - | 6313 | 2647 | 73.5 | 53.8 | 62.1
Ours | 6597 | 1692 | 859 | 3428 | 65.8 | 88 | 75.3

###  4.4 Evaluation of Polyp Segmentation Performance

To quantitatively evaluate the effectiveness of the generated synthetic polyps
in improving the performance of the segmentation models, we apply commonly
used evaluation metrics: Jaccard index (J) also known as intersection over
union (IoU) and Dice similarity score (D) as follows:

$J(A.B)=\frac{\mid A\cap B\mid}{\mid A\cup B\mid}=\frac{\mid A\cap B\mid}{\mid
A\mid+\mid B\mid-\mid A\cap B\mid}.\\\ Dice(A.B)=\frac{2\mid A\cap B\mid}{\mid
A\mid+\mid B\mid}.\quad\quad\quad\quad\quad\quad\quad\quad\quad\>\>$ (3)

where $A$ refers to a generated output by the models and $B$ is the
corresponding ground-truth of $A$.

Table 3 demonstrates performance evaluation of the three segmentation models
where two scenarios in terms of training data, i.e., original and original
plus 350 synthetic polyp images (combined training), are used to train the
models. As it can be seen from the table, the performance of the three
segmentation models improves when the generated synthetic polyps are added to
the original training data as follows: The performance of TernausNet-16
improved by 6.44% Jaccard and 9.66% Dice, and AlbuNet-34 by 4.11% Jaccard and
5.87% Dice, and MDeNetplus by 5.27% Jaccard, 6.47% Dice.

Table 3: Comparison of polyp segmentation performance for three different models using original training set and combined training set by generated polyp image. Models | TernausNet-16 | AlbuNet-34 | MDeNetplus
---|---|---|---
Jaccard | Dice | Jaccard | Dice | Jaccard | Dice
Original | 35.47 | 43.18 | 46.75 | 56.98 | 45.77 | 56.53
[17] | 37.86 | 47.57 | 50.48 | 60.31 | 50.49 | 60.77
Combined | 41.91 | 52.84 | 50.86 | 62.85 | 51.04 | 63

In Table 2, we use the same segmentation models to compare the improvement
capability of the synthetic polyps generated by this work and the synthetic
polyps generated by the method presented in [17]. For all three models, the
synthetic polyps of the current study show a better result improvement than
the synthetic polyps of [17] study. This might be because we are capable of
controlling the generation of synthetic polyps with various shapes and
features from the same input condition image. As a result, with the proposed
method, polyps with various characteristics can be generated. This is a
difference from [17] in which only the shape of the polyp can be controlled.

###  4.5 Limitations and Future Work

One natural question might be how performance changes as the number of
generated polyps increases. In Figure 10, the detection performance of the
RFCN ResNet101 model is evaluated when different numbers of synthetic polyp
images are added to the original training dataset. As it can be seen that the
model starts gradually improving its performance when it is exposed to
additional synthetic polyp images. However, the model reaches a saturation
performance improvement after a certain point even if more synthetic polyp
images are added to the training set. This is due to the limitation of this
method that is unable to introduce new unseen features. This method only
manipulates the existing features in the training dataset used to train the
GAN model and tries to reuse the same set of features to generate new-looking
synthetic polyps. In other words, the proposed framework cannot improve data
distribution and is unable to add new features to the training data, it only
leads to an image-level transformation and feature-level manipulation. This is
a common limitation of GAN-based image generation models.

For future work, studies that generate polyps with new characteristics not
included in the training data will be needed. Performance can further be
improved if characteristics of two/multiple polyps can be fused into a single
synthetic polyp. This could be done by fusing obtained abstract features from
different polyps or by manipulating the features by adding random noise. Both
approaches could be performed in the latent space. Alternatively, we can apply
a method such as CAN (Creative Adversarial Networks) [29] that adjusts the
loss term of GAN to generate a creative image by maximizing deviation from
established styles and minimizing deviation from the distribution of several
datasets.

$Original$$150$$350$$550$$750$$55$$60$$65$$70$$75$$80$$85$$90$$95$Number of
imagesPerformance (%) SensitivityPrecisionf1-score Fig. 10: Performance
improvement with different number of synthetic polyp images added to the
training data

##  5 Conclusion

It is expensive to acquire images of various types of polyps that can be used
for training deep learning based automatic polyp detection and segmentation
models. We proposed a synthetic polyp image generation framework based on the
conditional GAN architectures. In the proposed framework, we converted a given
polyp image into a negative image and then the negative image back into a new-
looking synthetic polyp image using the same networks for both tasks. In the
previous polyp generation studies, a single value was assigned to the polyp
mask in the condition input image. This led the models to face difficulties to
generate various polyps with different characteristics. In this study, we
attempted to overcome this shortcoming. We developed a framework to generate
various polyps with different features by controlling the value of the polyp
masks in the input condition images. The experimental results showed that the
proposed framework could generate polyps with various characteristics similar
to the real ones. In addition, it was confirmed that the polyp detection and
segmentation performance could be improved when the generated synthetic polyp
images were used to train several detection and segmentation models.

## Acknowledgements

This research was supported by the Research Funds of Mokpo National University
in 2020. The authors would like to thank Jacob Bergsland at Oslo University
Hospital for his valuable comments.

## References

  * [1] F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, A. Jemal, Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA: a cancer journal for clinicians 68 (6) (2018) 394–424.
  * [2] A. G. Zauber, S. J. Winawer, M. J. O’Brien, I. Lansdorp-Vogelaar, M. van Ballegooijen, B. F. Hankey, W. Shi, J. H. Bond, M. Schapiro, J. F. Panish, et al., Colonoscopic polypectomy and long-term prevention of colorectal-cancer deaths, N Engl J Med 366 (2012) 687–696.
  * [3] S. B. Ahn, D. S. Han, J. H. Bae, T. J. Byun, J. P. Kim, C. S. Eun, The miss rate for colorectal adenoma determined by quality-adjusted, back-to-back colonoscopies, Gut and liver 6 (1) (2012) 64.
  * [4] S. A. Karkanis, D. K. Iakovidis, D. E. Maroulis, D. A. Karras, M. Tzivras, Computer-aided tumor detection in endoscopic video using color wavelet features, IEEE transactions on information technology in biomedicine 7 (3) (2003) 141–152.
  * [5] N. Tajbakhsh, J. Y. Shin, S. R. Gurudu, R. T. Hurst, C. B. Kendall, M. B. Gotway, J. Liang, Convolutional neural networks for medical image analysis: Full training or fine tuning?, IEEE transactions on medical imaging 35 (5) (2016) 1299–1312.
  * [6] Y. Shin, I. Balasingham, Automatic polyp frame screening using patch based combined feature and dictionary learning, Computerized Medical Imaging and Graphics 69 (2018) 33–42.
  * [7] J. Bernal, N. Tajkbaksh, F. J. Sanchez, B. J. Matuszewski, H. Chen, L. Yu, Q. Angermann, O. Romain, B. Rustad, I. Balasingham, et al., Comparative validation of polyp detection methods in video colonoscopy: results from the miccai 2015 endoscopic vision challenge, IEEE transactions on medical imaging 36 (6) (2017) 1231–1249.
  * [8] L. Yu, H. Chen, Q. Dou, J. Qin, P. A. Heng, Integrating online and offline three-dimensional deep learning for automated polyp detection in colonoscopy videos, IEEE journal of biomedical and health informatics 21 (1) (2016) 65–75.
  * [9] Y. Shin, H. A. Qadir, L. Aabakken, J. Bergsland, I. Balasingham, Automatic colon polyp detection using region based deep cnn and post learning approaches, IEEE Access 6 (2018) 40950–40962.
  * [10] H. A. Qadir, I. Balasingham, J. Solhusvik, J. Bergsland, L. Aabakken, Y. Shin, Improving automatic polyp detection using cnn by exploiting temporal dependency in colonoscopy video, IEEE journal of biomedical and health informatics 24 (1) (2019) 180–193.
  * [11] H. A. Qadir, Y. Shin, J. Solhusvik, J. Bergsland, L. Aabakken, I. Balasingham, Polyp detection and segmentation using mask r-cnn: Does a deeper feature extractor cnn always perform better?, in: 2019 13th International Symposium on Medical Information and Communication Technology (ISMICT), IEEE, 2019, pp. 1–6.
  * [12] H. A. Qadir, Y. Shin, J. Solhusvik, J. Bergsland, L. Aabakken, I. Balasingham, Toward real-time polyp detection using fully cnns for 2d gaussian shapes prediction, Medical Image Analysis 68 (2021) 101897.
  * [13] J. Bernal, F. J. Sánchez, G. Fernández-Esparrach, D. Gil, C. Rodríguez, F. Vilariño, Wm-dova maps for accurate polyp highlighting in colonoscopy: Validation vs. saliency maps from physicians, Computerized Medical Imaging and Graphics 43 (2015) 99–111.
  * [14] J. Silva, A. Histace, O. Romain, X. Dray, B. Granado, Toward embedded detection of polyps in wce images for early diagnosis of colorectal cancer, International journal of computer assisted radiology and surgery 9 (2) (2014) 283–293.
  * [15] Q. Angermann, J. Bernal, C. Sánchez-Montes, M. Hammami, G. Fernández-Esparrach, X. Dray, O. Romain, F. J. Sánchez, A. Histace, Towards real-time polyp detection in colonoscopy videos: Adapting still frame-based methodologies for video sequences analysis, in: Computer Assisted and Robotic Endoscopy and Clinical Image-Based Procedures, Springer, 2017, pp. 29–41.
  * [16] F. Mahmood, R. Chen, N. J. Durr, Unsupervised reverse domain adaptation for synthetic medical images via adversarial training, IEEE transactions on medical imaging 37 (12) (2018) 2572–2581.
  * [17] Y. Shin, H. A. Qadir, I. Balasingham, Abnormal colon polyp image synthesis using conditional adversarial networks for improved detection performance, IEEE Access 6 (2018) 56007–56017.
  * [18] V. de Almeida Thomaz, C. A. Sierra-Franco, A. B. Raposo, Training data enhancements for improving colonic polyp detection using deep convolutional neural networks, Artificial Intelligence in Medicine 111 (2021) 101988.
  * [19] L. Yuan, C. Ruan, H. Hu, D. Chen, Image inpainting based on patch-gans, IEEE Access 7 (2019) 46411–46421.
  * [20] A. A. Shvets, V. I. Iglovikov, A. Rakhlin, A. A. Kalinin, Angiodysplasia detection and localization using deep convolutional neural networks, in: 2018 17th IEEE international conference on machine learning and applications (ICMLA), IEEE, 2018, pp. 612–617.
  * [21] P. Isola, J.-Y. Zhu, T. Zhou, A. A. Efros, Image-to-image translation with conditional adversarial networks, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 1125–1134.
  * [22] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, A. Courville, Improved training of wasserstein gans. in advances in neural information processing systems.
  * [23] T. R. Shaham, T. Dekel, T. Michaeli, Singan: Learning a generative model from a single natural image, in: Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 4570–4580.
  * [24] S. Ren, K. He, R. Girshick, J. Sun, Faster r-cnn: Towards real-time object detection with region proposal networks, Advances in neural information processing systems 28 (2015) 91–99.
  * [25] J. Dai, Y. Li, K. He, J. Sun, R-fcn: Object detection via region-based fully convolutional networks, in: Advances in neural information processing systems, 2016, pp. 379–387.
  * [26] V. Iglovikov, A. Shvets, Ternausnet: U-net with vgg11 encoder pre-trained on imagenet for image segmentation, arXiv preprint arXiv:1801.05746.
  * [27] O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, in: International Conference on Medical image computing and computer-assisted intervention, Springer, 2015, pp. 234–241.
  * [28] J. Huang, V. Rathod, C. Sun, M. Zhu, A. Korattikara, A. Fathi, I. Fischer, Z. Wojna, Y. Song, S. Guadarrama, et al., Speed/accuracy trade-offs for modern convolutional object detectors, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 7310–7311.
  * [29] A. Elgammal, B. Liu, M. Elhoseiny, M. Mazzone, Can: Creative adversarial networks, generating" art" by learning about styles and deviating from style norms, arXiv preprint arXiv:1706.07068.