# Fabrication, characterization, and simulation of glass devices with AlN-
thin-film-transducers for excitation of ultrasound resonances

André G. Steckel<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark, DTU Physics Building 309, DK-2800 Kongens Lyngby,
Denmark  Henrik Bruus<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark, DTU Physics Building 309, DK-2800 Kongens Lyngby,
Denmark  Paul Muralt<EMAIL_ADDRESS>Materials Science, Station 12, EPFL,
1015 Lausanne, Switzerland  PIEMACS Sàrl, EPFL Innovation Parc, Bâtiment C,
1015-Lausanne Switzerland  Ramin Matloub<EMAIL_ADDRESS>PIEMACS
Sàrl, EPFL Innovation Parc, Bâtiment C, 1015-Lausanne Switzerland

(16 November 2020)

###### Abstract

We present fabrication of 570-µm-thick, millimetric-sized soda-lime-silicate
float glass blocks with a 1-µm-thick AlN-thin-film piezoelectric transducer
sandwiched between thin metallic electrodes and deposited on the top surface.
The electromechanical properties are characterized by electrical impedance
measurements in the frequency range from 0.1 to 10 MHz with a peak-to-peak
voltage of 0.5 V applied to the electrodes. We measured the electrical
impedance spectra of 35 devices, all of width 2 mm, but with 9 different
lengths ranging from 2 to 6 mm and with 2-7 copies of each individual
geometry. Each impedance spectrum exhibits many resonance peaks, of which we
carefully measured the 5 most prominent ones in each spectrum. We compare the
resulting 173 experimental resonance frequencies with the simulation result of
a finite-element-method model that we have developed. When using material
parameters from the manufacturer, we obtain an average relative deviation of
the 173 simulated resonance frequencies from the experimental ones of
$(-4.2\pm 0.04)$%. When optimizing the values of the Young’s modulus and the
Poisson ratio of the float glass in the simulation, this relative deviation
decreased to $(-0.03\pm 0.04)$%. Our results suggest a method for an accurate
_in-situ_ determination of the acoustic parameters at ultrasound frequencies
of any elastic solid onto which a thin-film transducer can be attached.

##  I Introduction

Aluminum nitride (AlN) is one of the most commonly used materials in
integrated thin-film piezoelectric transducers for actuating electromechanical
systems (MEMS). Due to their low dielectric and mechanical loss tangent, their
structural and chemical stability, as well as their compatibility with
standard silicon-based CMOS microfabrication techniques, AlN-sputtered thin
films are commercially used in thin film bulk-wave acoustic resonator filters
Ruby (2017). The academic literature reports applications of AlN thin-film
transducers as RF filters Dubois and Muralt (1999), contour mode resonators
Piazza  _et al._ (2006), switches Zaghloul and Piazza (2014); Sinha  _et al._
(2009), suspended microchannel resonators De Pastina  _et al._ (2018), and
accelerometers Olsson, III  _et al._ (2009). Detailed studies on AlN-thin-
film-actuated high-tone bulk acoustic resonators with high quality factors
have been studied in the 0.3 - 3 GHz range on 350-µm-thick substrates of
sapphire, crystal quartz, fused silica, and silicon Zhang  _et al._ (2006),
and on 30-µm-thick Si membranes Masson  _et al._ (2007); Fujikura  _et al._
(2008).

Because AlN is a nonferroelectric polar material, a reorientation of the polar
axis is not possible. The material growth process therefore has to provide for
a textured structure that includes the alignment of the polar directions
Muralt (2008), since otherwise the random orientation of the grains
constituting the AlN thin film would result in a zero average piezoelectric
effect. In this work, we therefore follow the process developed in our
previous work on pure Howell  _et al._ (2019) and Sc-doped Matloub  _et al._
(2011, 2013) AlN films, and select Pt/Ti as the seed layer for the bottom
electrode stack as it promotes nucleation of (002)-oriented AlN grains.

The past few years, piezoelectric thin films made of nonferroelectric polar
materials, such as zinc oxide (ZnO) and AlN, have found a broad range of lab-
on-chip applications such as biosensing, particle/cell concentrating, sorting,
patterning, pumping, mixing, nebulization and jetting. Integrated devices with
acoustic transducers, sensors, and microfluidic channels have been fabricated
by depositing these piezoelectric films onto wide range of materials beyond
the usual silicon substrates such as ceramics, diamond, glass, and more
recently also polymer, metallic foils and flexible glass for making flexible
devices. Such thin-film acoustic wave devices have great potential for
implementing integrated, disposable, or bendable lab-on-a-chip devices into
various sensing and actuating applications Fu  _et al._ (2017). To advance
such applications, we focus in this work on AlN-thin-film transducers
deposited on glass substrates and driven in the 0.1 - 10 MHz range.

The deposition of AlN thin films and their texture control are more difficult
than in the case with ZnO films Lee  _et al._ (2007), and their integration
into MEMS is not a trivial task, as they tend to have high in-plane stress
with a sharp transition around zero stress Dubois and Muralt (2001). However,
the thermal conductivity of AlN is more than 3 times larger than of ZnO, which
enables the operation of AlN films at a higher power level, an important
feature for applications in RF technology. AlN films are normally sputter-
deposited at moderate temperatures (typically 200-400
${}^{\circ}\\!\textrm{C}$) to achieve an optimal performance, although “room
temperature” deposition with a local temperature near the target of around 150
${}^{\circ}\\!\textrm{C}$ is sometimes used Lee  _et al._ (2002). Also
deposition conditions, especially oxygen or moisture in vacuum chambers, have
significant effects on AlN film growth and microstructure. Growing AlN films
thicker than a few microns is particularly challenging because of the
potentially large film stress and the tendency to form microcracks.

In this work, we used our titanium-seed-layer technique Howell  _et al._
(2019) to improve the fabrication of AlN-thin-film transducers deposited on
glass substrates for future use in biosensing and lab-on-a-chip applications.
In Section II, we describe the design of the devices, followed in Section III
by a presentation of the fabrication and characterization processes. In
Section IV, we then characterize the acoustic properties of the devices by
measuring electrical impedance spectra in the frequency range from 0.1 to 10
MHz. The numerical model and its finite-element-method implementation is
described in Section V. The comparison of the simulation results with the
experiments presented in Section VI are in so good agreement that our
numerical handling of the measured impedance spectra can be turned into an
improved _in-situ_ determination of the acoustic properties of the glass
substrate. Finally, concluding remarks are presented in Section VII.

Figure 1:  (a) Design of the electrode mask for the lithography step in the
fabrication process that defines the gold electrodes on the top surface of the
glass wafer. The relevant devices are taken from row A - G and column 1 - 9.
The inset show row B column 1 - 9. Devices with the same column index are
nominally identical. (b) A photograph of the fabricated wafer before dicing it
in individual devices. (c) Microscope images of the actual device B2 and B9.
(d) Rendering of the two 2-mm-wide and 0.58-mm-high devices B2 (2.5 mm long)
and B9 (6 mm long), showing the outline of the 1-mm-wide top electrode split
in two halves by a 40-µm-wide gap.

##  II Device design

In Fig. 1 is shown the design of the devices under study. This design is
chosen to ease the microfabrication, the experimental characterization, and
the subsequent validation of a finite-element method model of the system. All
devices consist of a simple box-shaped 570-µm-thick and 2-mm-wide block of
soda-lime-silicate (SLS) float glass Sch  of length $L$ diced from a wafer.
$L$ is in the range from 2 to 6 mm in steps of 0.5 mm, thus resulting in 9
groups of nominally identical devices. A 1-µm-thick AlN-thin-film
piezoelectric transducer, sandwiched between sub-µm-thick metal electrodes, is
deposited on top surface of the glass block. To ease the dicing process and to
avoid short circuiting between the two electrodes, the top electrode of the
transducer has its width reduced to 1 mm and its length to
$L_{\mathrm{te}}=L-1$ mm. Moreover, to allow for anti-symmetric voltage
actuation, the top electrode is split in half by a 40-µm-wide gap
perpendicular to the length direction. As listed in Table 1, there are up to
seven nominally identical copies A, B, C, D, E, F, and G in each group of
devices.

Table 1:  A list of the 35 devices under study. The rectangular glass blocks are all 0.58 mm thick, 2 mm wide, and with a length $L$ in the range from 2 to 6 mm. The deposited AlN-thin-film piezoelectric transducer is 1 µm thick and sandwiched between 0.1-µm-thick metal electrodes. The top electrode has the width 1 mm and length $L_{\mathrm{te}}=L-1$ mm. All devices with the same column index are nominally identical. Device | Device | $L$ | $L_{\mathrm{te}}$
---|---|---|---
column | row | [mm] | [mm]
1 | C, D, G | 2.0 | 1.0
2 | C, D, E, F, G | 2.5 | 1.5
3 | A, B, C, G | 3.0 | 2.0
4 | A, B, D | 3.5 | 2.5
5 | A, B, C | 4.0 | 3.0
6 | B, C | 4.5 | 3.5
7 | A, B, C, D, E | 5.0 | 4.0
8 | A, B, C, D, E, F, G | 5.5 | 4.5
9 | C, D, E | 6.0 | 5.0

##  III Fabrication and characterization of the piezoelectric thin-film
transducer

The device fabrication is based on a 4-inch SLS float glass wafer (SCHOTT
Suisse SA, Yverdon). This glass has a low thermal expansion, a high thermal
shock resistance, and an ability to withstand temperatures up to
$490~{}^{\circ}\\!\textrm{C}$ for long periods, so by not raising the process
temperatures above $300~{}^{\circ}\\!\textrm{C}$, we ensure a high quality and
reproducibility in our fabrication process. Following a microfabrication
process similar to the one described in Ref. Howell  _et al._ (2019) and
sketched in Fig. 2, the first step is cleaning of the glass substrate with
isopropyl alcohol and acetone to remove organic contaminations. Additional
plasma cleaning was performed right before the deposition of the bottom
electrode in the deposition tool. We prepared textured Ti, Pt, AlN thin films
using a dedicated cluster sputtering chamber (Pfeiffer Vacuum), and then
sputter-deposited 10 nm Ti at $300~{}^{\circ}\\!\textrm{C}$ serving as an
adhesion layer between the glass substrate and the Pt bottom electrode. Next,
the 100-nm thick Pt electrode was deposited with the same tool in a different
chamber also at $300~{}^{\circ}\\!\textrm{C}$ and without breaking the vacuum.
Successively, $1~{}\textrm{\textmu{}m}$ (0001)-textured AlN was deposited
using a dedicated reactive, pulsed, direct current magnetron sputtering
chamber at $300~{}^{\circ}\\!\textrm{C}$, and on top of the AlN film 10 nm Cr
and 150 nm Au were sputter deposited as top electrode with Cr being the
adhesion layer between the AlN and Au. As the next process step, standard
photolithography was exploited to transfer the intended pattern on to the top
electrode, and this pattern was realized by ion-beam etching (Veeco). Finally,
the wafer was diced into the smaller chips ready for use, as shown in Fig.
1(c).

After fabrication and dicing, the quality of the resulting devices was checked
by visual inspection in a microscope, and their dimensions were checked by
measurements using a digital caliper (Vernier Caliper, Cocraft) with 0.03-mm
accuracy and 0.01-mm repeatability. 29 representative devices were selected.
Their lengths were found to range from $(2034\pm 5)~{}\textrm{\textmu{}m}$ to
$(6034\pm 5)~{}\textrm{\textmu{}m}$ in steps of 500 µm, the widths were
$(2025\pm 5)~{}\textrm{\textmu{}m}$, and the heights were $(570\pm
3)~{}\textrm{\textmu{}m}$. The uncertainties represent one standard deviation
of several repeated measurements, and the relative uncertainty ranges from
0.6% for the height, down to 0.1% for the longest length.

Figure 2:  Sketch of the fabrication process, where a stack consisting of a
bottom electrode, the AlN thin film, and the top electrode is deposited on the
glass block. (a) Deposition of a 10-nm-thick Ti adhesion layer (not shown)
followed by a 100-nm-thick platinum layer. (b) Deposition of the 1-µm-thick
AlN film. (c) Deposition of a 10-nm-thick Cr adhesion layer (not shown)
followed by a 150-nm-thick gold layer. (d) Deposition and development of
photoresist. (e) Patterning of the top electrode by ion-beam etching of the
exposed gold followed by removal of the photoresist.

##  IV Measured electrical impedance

The electrical impedance of the 35 devices listed in Table 1 were measured as
a function of frequency from from 0.1 to 3 MHz in steps of 500 Hz using an
Agilent 4294A Precision Impedance Analyzer (Agilent Technologies AG, Basel,
Switzerland) with a 0.5-Vpp ac-voltage applied in anti-phase between the two
halves of the split top electrode. The bottom electrode was left floating.

Figure 3:  (a) Three photographs of the spectrum analyzer and the attached metal box that provides electrical shielding of the device, but here with its lid removed to reveal the device mounted inside. A zoom-in on the interior of the metal box shows the coax connectors and the wires leading to the device. A further zoom-in shows the spring-loaded electrodes attached to the device. (b) An end-view sketch of the spring-loaded electrodes attached to the device consisting of the glass block (beige) and the AlN thin-film transducer (purple) with a split top electrode (red-blue). (c) The impedance spectra showing the magnitude $|Z(f)|$ (red curves) and the phase $\theta(f)$ (blue curves) for the seven nominally identical devices A8, B8, C8, D8, E8, F8, and G8 shown in the inset. The phase curves differ as the devices were not discharged initially. (d) Same as for panel (c) but for the five nominally identical devices C2, D2, E2, F2, and G2. Here, the phase curves are nearly identical as the devices were discharged initially. Table 2:  For each of the nine device geometries, labeled by column index $k=1$ \- 9 in Table 1, five resonance peaks $(k,n)$ of increasing frequency with $n=1$ \- 5 were selected the impedance spectrum for comparing simulations with measurements. Peak $(4,1)$ and $(4,2)$ marked by an asterisk were not measured on device B in column $n=4$. The experimental value $f^{\mathrm{exp}}_{k,n}$ is the average resonance frequency of peak $n$ over all devices (row index A - G in Table 1) with column index $k$. The simulation result $f^{\mathrm{sim,A}}_{k,n}$ is the corresponding resonance frequency obtained by the numerical modeling described in Section V and using the parameter values given in the literature. The simulation result $f^{\mathrm{sim,B}}_{k,n}$ is obtained by using the optimized values listed in Table 4 of Young’s modulus and Poisson’s ratio for the glass block. Finally, for X = A or B, $\Delta^{\mathrm{X}}_{k,n}=(f^{\mathrm{sim,X}}_{k,n}-f^{\mathrm{exp}}_{k,n})/f^{\mathrm{exp}}_{k,n}$ is the relative deviation of the simulation result from the experimental value for resonance peak $n$ in device geometry $k$. Peak | $f^{\mathrm{exp}}_{k,n}\;$ | $f^{\mathrm{sim,A}}_{k,n}$ | $f^{\mathrm{sim,B}}_{k,n}$ | $\;\;\Delta^{\mathrm{A}}_{k,n}$ | $\;\;\Delta^{\mathrm{B}}_{k,n}$ | Peak | $f^{\mathrm{exp}}_{k,n}\;$ | $f^{\mathrm{sim,A}}_{k,n}$ | $f^{\mathrm{sim,B}}_{k,n}$ | $\;\;\Delta^{\mathrm{A}}_{k,n}$ | $\;\;\Delta^{\mathrm{B}}_{k,n}$
---|---|---|---|---|---|---|---|---|---|---|---
$k,n$ | [MHz] | [MHz] | [MHz] | [%] | [%] | $k,n$ | [MHz] | [MHz] | [MHz] | [%] | [%]
1, 1 | $0.943$ | $0.883$ | $0.926$ | $-6.3$ |  $-1.7$ | 5, 4 | $2.712$ | $2.607$ | $2.721$ | $-3.9$ | $0.3\rule{0.0pt}{10.76385pt}$
1, 2 | $1.424$ | $1.369$ | $1.422$ | $-3.9$ |  $-0.2$ | 5, 5 | $2.901$ | $2.780$ | $2.918$ | $-4.2$ | $0.6$
1, 3 | $2.065$ | $1.958$ | $2.054$ | $-5.2$ |  $-0.5$ | 6, 1 | $1.425$ | $1.368$ | $1.422$ | $-4.0$ | $-0.2$
1, 4 | $2.717$ | $2.618$ | $2.720$ | $-3.7$ |  $0.1$ | 6, 2 | $1.537$ | $1.470$ | $1.538$ | $-4.3$ | $0.1$
1, 5 | $2.896$ | $2.766$ | $2.906$ | $-4.5$ |  $0.3$ | 6, 3 | $1.606$ | $1.540$ | $1.601$ | $-4.1$ | $-0.3$
2, 1 | $0.737$ | $0.701$ | $0.736$ | $-4.9$ |  $-0.2$ | 6, 4 | $2.741$ | $2.641$ | $2.739$ | $-3.7$ | $-0.1$
2, 2 | $1.418$ | $1.365$ | $1.419$ | $-3.7$ |  $0.1$ | 6, 5 | $2.910$ | $2.748$ | $2.922$ | $-5.6$ | $0.4$
2, 3 | $1.791$ | $1.706$ | $1.788$ | $-4.8$ |  $-0.2$ | 7, 1 | $0.321$ | $0.306$ | $0.322$ | $-4.8$ | $0.3$
2, 4 | $1.915$ | $1.830$ | $1.906$ | $-4.5$ |  $-0.4$ | 7, 2 | $1.440$ | $1.379$ | $1.431$ | $-4.2$ | $-0.6$
2, 5 | $2.942$ | $2.802$ | $2.939$ | $-4.7$ |  $-0.1$ | 7, 3 | $1.561$ | $1.495$ | $1.555$ | $-4.2$ | $-0.4$
3, 1 | $0.612$ | $0.582$ | $0.611$ | $-4.9$ |  $-0.2$ | 7, 4 | $2.733$ | $2.632$ | $2.731$ | $-3.7$ | $-0.1$
3, 2 | $1.403$ | $1.348$ | $1.409$ | $-3.9$ |  $0.4$ | 7, 5 | $2.905$ | $2.783$ | $2.920$ | $-4.2$ | $0.5$
3, 3 | $1.489$ | $1.424$ | $1.482$ | $-4.3$ |  $-0.5$ | 8, 1 | $1.438$ | $1.375$ | $1.426$ | $-4.4$ | $-0.9$
3, 4 | $1.752$ | $1.681$ | $1.754$ | $-4.0$ |  $0.1$ | 8, 2 | $1.533$ | $1.470$ | $1.531$ | $-4.1$ | $-0.2$
3, 5 | $2.735$ | $2.640$ | $2.737$ | $-3.5$ |  $0.1$ | 8, 3 | $1.834$ | $1.760$ | $1.835$ | $-4.0$ | $0.1$
4, 1∗ | $0.522$ | $0.498$ | $0.527$ | $-4.6$ |  $0.9$ | 8, 4 | $2.718$ | $2.618$ | $2.725$ | $-3.7$ | $0.3$
4, 2∗ | $0.544$ | $0.519$ | $0.541$ | $-4.6$ |  $-0.4$ | 8, 5 | $2.903$ | $2.786$ | $2.922$ | $-4.0$ | $0.7$
4, 3 | $1.431$ | $1.374$ | $1.426$ | $-4.0$ |  $-0.3$ | 9, 1 | $1.429$ | $1.373$ | $1.426$ | $-3.9$ | $-0.2$
4, 4 | $1.667$ | $1.600$ | $1.669$ | $-4.0$ |  $0.1$ | 9, 2 | $1.511$ | $1.451$ | $1.511$ | $-4.0$ | $0.0$
4, 5 | $2.725$ | $2.627$ | $2.728$ | $-3.6$ |  $0.1$ | 9, 3 | $1.765$ | $1.696$ | $1.768$ | $-3.9$ | $0.2$
5, 1 | $1.429$ | $1.372$ | $1.425$ | $-4.0$ |  $-0.3$ | 9, 4 | $2.089$ | $2.018$ | $2.111$ | $-3.4$ | $1.1$
5, 2 | $1.609$ | $1.545$ | $1.611$ | $-4.0$ |  $0.1$ | 9, 5 | $2.744$ | $2.644$ | $2.746$ | $-3.7$ | $0.1$
5, 3 | $2.104$ | $2.015$ | $2.107$ | $-4.2$ |  $0.2$ |  |  |  |  |  |

To reduce the electrical noise, a given device under study was placed in a
grounded metal box, from where it was connected by short coax cables to the
impedance analyzer by pressing a spring-loaded needle electrodes to each half
of the split top electrode, see Fig. 3(a)-(b). The setup was also designed to
minimize undesired parasitic components: the parasitic inductance in the
system was reduced to somewhere around $L_{\mathrm{ext}}=(28\pm 15)$ nH, the
parasitic resistance was roughly $R_{\mathrm{ext}}=(5\pm
3)~{}\mathrm{\Omega}$, and the parasitic capacitance was negligible compared
to the intrinsic capacitance $C_{\mathrm{dev}}$ of the device formed by the
electrodes separated by the thin AlN film, and ranging from 18 pF for the
smallest device to 86 pF for the largest. This design ensured that the
resonance frequency
$f_{\mathrm{ext}}=\frac{1}{2\pi}(L_{\mathrm{ext}}C_{\mathrm{dev}})^{-\frac{1}{2}}$
involving the external inductance was pushed up in value,
$100~{}\textrm{MHz}\lesssim f_{\mathrm{ext}}\lesssim 220~{}\textrm{MHz}$, much
larger than 10 MHz, and therefore not interfering with the intrinsic
resonances of the device under study. Finally, the device under study was
placed on a folded, soft piece of tissue paper to acoustically isolate from
the environment.

Examples of the recorded 35 impedance spectra
$Z(f)=|Z|\>\mathrm{e}^{\mathrm{i}\theta}$ are shown in Fig. 3(c) and (d),
where we plot the magnitude $|Z(f)|$ and the phase $\theta(f)$. We first note
that the intrinsic resonances of the devices are most clearly visible in the
phase plots $\theta(f)$. Due to the parasitic impedances of the system, the
base level of $\theta(f)$ is a linear increase from the $-90^{\circ}$ of an
ideal capacitor in the dc limit to a value between $-89.7^{\circ}$ and
$-87.8^{\circ}$ at 10 MHz. The device resonances appear as sharp peaks in the
phase-curve with a width of around 0.01 times the resonance frequency itself
and with an amplitude ranging from 0.1∘ to 3.5∘. This technique is a variant
of standard resonant ultrasound spectroscopy (RUS) Radovic  _et al._ (2004);
Cachiaras  _et al._ . However, in contrast to conventional RUS, the ultrasound
transducers in our case are an integrated but tiny part of the structures to
be characterized.

Next, we note that position and amplitudes of the resonance peaks are nearly
identical for samples with nominally identical geometry from different parts
of the wafer, see Fig. 1(a), for the about 30 peaks easily identified in the
range from 0.1 to 10 MHz. About 40% of the devices selected from column 1 - 9
and row A - G had a short circuit between one or both of the top electrodes
with the bottom electrode, so they were discarded from the study. This is the
reason that only column 8 in Table 1 all seven fabricated devices. We found
that the slope of the base level of the phase spectra $\theta(f)$ could vary
even for nominally identical devices, as seen in Fig. 3(d). It turned out that
this variance was due to spurious charging effects, perhaps due to the
floating electrode. This effect could be avoided by shortly short-circuiting
the electrodes with a piece of metallic wire, and thus discharging the system.
After discharging, the phase spectra became nearly identical for all devices
with the same geometry, see Fig. 3(e). When taking the mentioned precautions,
the impedance spectra exhibits a large degree of repeatability and
reproducibility.

For each of the nine different device geometries, labeled by column index
$k=1$ \- 9 in Table 1, among our 35 selected devices, we measure the resonance
frequency of the five most prominent resonance peaks, labeled by $n=1$ \- 5 in
the impedance-phase spectra $\theta(f)$. This resulted in $5\times 35-2=173$
resonance frequencies (for device 4B we did not measure all five peaks, but
only the three highest). In Table 2 we list the experimentally determined
values of the resonance frequencies, but to avoid showing all 173 of them, we
only show the average resonance frequency $f^{\mathrm{exp}}_{k,n}$ of peak $n$
over all devices (row index A - G in Table 1) with column index $k$, which
results in $9\times 5=45$ entries.

##  V Modeling and numerical simulation

We want to compare the experimental results for the impedance spectra and
resonance frequencies with numerical simulation. We therefore use the three-
dimensional model for acoustofluidic systems recently developed by Skov et al.
Skov  _et al._ (2019a). The model couples the mechanical displacement field
$\bm{u}$ and the electric potential $\varphi$ in a system driven by a time-
harmonic electric potential,
$\tilde{\varphi}=\varphi_{0}\>\mathrm{e}^{-{\mathrm{i}\omega t}}$ with angular
frequency $\omega=2\pi f$ and frequency $f$, applied to the interfaces between
the metal electrodes and the piezoelectric AlN-thin-film transducer. This
time-harmonic boundary condition excites the time-harmonic fields: the
electric potential $\tilde{\varphi}(\bm{r},t)$ in the AlN-thin-film transducer
and the displacement $\tilde{\bm{u}}(\bm{r},t)$ in all solids of the system,

$\displaystyle\tilde{\varphi}(\bm{r},t)=\varphi(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega
t}},\qquad\tilde{\bm{u}}(\bm{r},t)=\bm{u}(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega
t}}.$ (1)

Fields with a tilde are time-harmonic, whereas those without a tilde are
complex-valued amplitudes.

The motion of an elastic solid with a given density $\rho$ (in the model it is
the glass, the metal electrodes, and the AlN-thin-film transducer) is
described by the displacement field $\bm{u}$, whereas the electrodynamics is
described by the electric displacement field $\bm{D}$. The governing equations
for $\bm{u}$ and $\bm{D}$ are Cauchy’s elastodynamic equation and Maxwell’s
quasi-electrostatic equation, $\bm{\nabla}\cdot\bm{D}=0$,

$-\rho\omega^{2}(1+\mathrm{i}\Gamma_{\mathrm{sl}})\>\bm{u}=\bm{\nabla}\cdot\bm{\sigma},\quad\bm{\nabla}\cdot\big{[}-(1+\mathrm{i}\Gamma_{\varepsilon})\bm{\varepsilon}\cdot\bm{\nabla}\varphi\big{]}=0.$
(2)

Here, $\Gamma_{\mathrm{sl}},\Gamma_{\varepsilon}\ll 1$ are weak damping
coefficients, $\bm{\varepsilon}$ is the dielectric tensor, and $\bm{\sigma}$
is the stress tensor, which is coupled to $\bm{u}$ through a stress-strain
relation depending on the material-dependent elastic moduli $C_{ik}$.

For the isotropic glass and metal electrodes, the electric field is
unimportant, whereas the relation between the stress tensor $\sigma_{ik}$ and
strain components $\frac{1}{2}(\partial_{i}u_{k}+\partial_{k}u_{i})$ is given
in the compact Voigt representation as

$\left(\begin{array}[]{c}\\!\\!\sigma_{xx}\\!\\!\\\\[2.84526pt]
\\!\\!\sigma_{yy}\\!\\!\\\\[2.84526pt] \\!\\!\sigma_{zz}\\!\\!\\\\[2.84526pt]
\hline\cr\\!\\!\sigma_{yz}\\!\\!\\\\[2.84526pt]
\\!\\!\sigma_{xz}\\!\\!\\\\[2.84526pt]
\\!\\!\sigma_{xy}\\!\\!\end{array}\right)\\!=\\!\left(\begin{array}[]{ccc|ccc}\\!C_{11}&C_{12}&C_{12}&0&\\!\\!0&\\!\\!0\\!\\!\\\\[2.84526pt]
\\!C_{12}&C_{11}&C_{12}&0&\\!\\!0&\\!\\!0\\!\\!\\\\[2.84526pt]
\\!C_{12}&C_{12}&C_{11}&0&\\!\\!0&\\!\\!0\\!\\!\\\\[2.84526pt]
\hline\cr\\!0&0&0&\\!C_{44}&\\!0&\\!\\!0\\!\\!\\\\[2.84526pt]
\\!0&0&0&0&\\!\\!\\!C_{44}&\\!\\!0\\!\\!\\\\[2.84526pt]
\\!0&0&0&0&0&\\!\\!\\!C_{44}\\!\\!\\!\end{array}\right)\\!\\!\left(\begin{array}[]{c}\\!\partial_{x}u_{x}\\\\[2.84526pt]
\\!\partial_{y}u_{y}\\\\[2.84526pt] \\!\partial_{z}u_{z}\\\\[2.84526pt]
\hline\cr\\!\\!\partial_{y}u_{z}\\!+\\!\partial_{z}u_{y}\\!\\!\\\\[2.84526pt]
\\!\\!\partial_{x}u_{z}\\!+\\!\partial_{z}u_{x}\\!\\!\\\\[2.84526pt]
\\!\\!\partial_{x}u_{y}\\!+\\!\partial_{y}u_{x}\\!\\!\end{array}\right).$ (3)

Here, $C_{ik}$ are the elastic moduli which are listed for SLS float glass in
Table 3, and $C_{12}=C_{11}-2C_{44}$.

In the piezoelectric AlN film, the stress tensor $\bm{\sigma}$ and the
electric displacement field $\bm{D}$ are given by the strain tensor
$\frac{1}{2}\big{[}\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{T}\big{]}$, the
electrical field $-\bm{\nabla}\varphi$, the electromechanical coupling
coefficients $e_{ik}$, the elastic moduli $C_{ik}$, and the electric
permittivities $\varepsilon_{ik}$ as

$\displaystyle\left(\begin{array}[]{c}\sigma_{xx}\\\ \sigma_{yy}\\\
\sigma_{zz}\\\ \hline\cr\sigma_{yz}\\\ \sigma_{xz}\\\ \sigma_{xy}\\\ \hline\cr
D_{x}\\\ D_{y}\\\ D_{z}\\\ \end{array}\right)=$ (13)
$\displaystyle\left(\begin{array}[]{c@{\:}c@{\:}c@{\:}|c@{\:}c@{\:}c@{\:}|c@{\:}c@{\:}c}\\!\\!C_{11}\hfil\>&C_{12}\hfil\>&C_{13}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!0\hfil\>&\\!\text{-}e_{31}\\!\\!\\\
\\!\\!C_{12}\hfil\>&C_{11}\hfil\>&C_{13}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!0\hfil\>&\\!\text{-}e_{31}\\!\\!\\\
\\!\\!C_{13}\hfil\>&C_{13}\hfil\>&C_{33}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!0\hfil\>&\\!\text{-}e_{33}\\!\\!\\\
\hline\cr\\!\\!0\hfil\>&0\hfil\>&0\hfil\>&\\!C_{44}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!\text{-}e_{15}\hfil\>&\\!0\\!\\!\\\
\\!\\!0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!C_{44}\hfil\>&0\hfil\>&\text{-}e_{15}\hfil\>&\\!0\hfil\>&\\!0\\!\\!\\\
\\!\\!0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!C_{66}\hfil\>&0\hfil\>&\\!0\hfil\>&\\!0\\!\\!\\\
\hline\cr\\!\\!0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!e_{15}\hfil\>&0\hfil\>&\varepsilon_{11}\hfil\>&\\!0\hfil\>&\\!0\\!\\!\\\
\\!\\!0\hfil\>&0\hfil\>&0\hfil\>&e_{15}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!\varepsilon_{11}\hfil\>&\\!0\\!\\!\\\
e_{31}\hfil\>&e_{31}\hfil\>&e_{33}\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&\\!0\hfil\>&\\!\varepsilon_{33}\\!\\!\\\
\end{array}\right)\\!\\!\left(\begin{array}[]{c}\partial_{x}u_{x}\\\
\partial_{y}u_{y}\\\ \partial_{z}u_{z}\\\
\hline\cr\\!\\!\\!\partial_{y}u_{z}\\!+\\!\partial_{z}u_{y}\\!\\!\\!\\\
\\!\\!\\!\partial_{x}u_{z}\\!+\\!\partial_{z}u_{x}\\!\\!\\!\\\
\\!\\!\\!\partial_{x}u_{y}\\!+\\!\partial_{y}u_{x}\\!\\!\\!\\\
\hline\cr-\partial_{x}\varphi\\\ -\partial_{y}\varphi\\\
-\partial_{z}\varphi\\\ \end{array}\right)\\!.$ (32)

Here, the symmetry of AlN leads to $C_{12}=C_{11}-2C_{66}$.

In terms of volume (mass) ratios, the device consists of 99.8% (99.5%) glass,
0.2% (0.3%) AlN film, and 0.0% (0.2%) metal electrodes. We have found that
within the numerical uncertainties, the simulation results are unchanged when
leaving out the metal electrodes. So to reduce the memory requirements and to
speed up computation times, we have left out the tiny metal electrodes in our
modeling. Consequently, the boundary conditions for the electric potential
$\varphi$ are: (1) on the interface between the bottom electrode and the AlN
film, the potential is grounded, $\varphi=0$; (2) on the interface between one
half of the split electrode and the AlN film, $\varphi=+\frac{1}{2}V_{0}$; (3)
on the interface between the other half of the split electrode and the AlN
film it is in anti-phase, $\varphi=-\frac{1}{2}V_{0}$; (4) on the sides walls
with normal vector $\bm{n}$ of the AlN film, the no-free-charge condition
applies, $\bm{n}\cdot\bm{\nabla}\varphi=0$. For the mechanical displacement
filed in the solids, the non-stress condition applied to all external
surfaces, $\bm{n}\cdot\bm{\sigma}=\bm{0}$.

Table 3:  Material parameters at $25~{}^{\circ}\\!\textrm{C}$ for the AlN-thin-film transducer and the SLS float-glass substrate. The frequency-dependent damping coefficients $\Gamma_{\mathrm{sl}}(f)$ and $\Gamma_{\varepsilon}(f)$ are obtained by fitting to the measured width of the resonance peaks. The tiny metal electrodes are not simulated in this study. Parameter | Value |  | Parameter | Value
---|---|---|---|---
_Thin film aluminum nitride_ , AlN Caro  _et al._ (2015)
$\rho$ | 3300 $\textrm{kg}\>\textrm{m}^{-3}$ |  | $\Gamma_{\mathrm{sl}}$ | 0.001
$C_{11}$ | 410.2 GPa |  | $C_{33}$ | 385.0 GPa
$C_{12}$ | 142.4 GPa |  | $C_{44}$ | 122.9 GPa
$C_{13}$ | 110.1 GPa |  | $C_{66}$ | 133.9 GPa
$e_{13}$ | $-1.0475~{}\textrm{C}\>\textrm{m}^{-2}$ |  | $e_{15}$ | $-0.39~{}\textrm{C}\>\textrm{m}^{-2}$
$e_{33}$ | 1.46 $\textrm{C}\>\textrm{m}^{-2}$ |  | $\Gamma_{\varepsilon}$ | $-\frac{1.43\,\textrm{MHz}+0.98f}{1\,\textrm{GHz}}$
$\epsilon_{11}$ | 9.21 $\epsilon_{0}$ |  | $\epsilon_{33}$ | 10.12 $\epsilon_{0}$
_SLS float-glass substrate_ Sch
$\rho$ | (2522 $\pm$ 15) $\textrm{kg}\>\textrm{m}^{-3}$ | (measured at EPFL)
$E$ |  $(72\pm 2)$ GPa |  | $s$ | $0.23\pm 0.01$
$C_{11}$ | 83.5 GPa |  | $C_{44}$ | 29.3 GPa
$C_{12}$ | 24.9 GPa |  | $\Gamma_{\mathrm{sl}}(f)$ |  $f^{-1}$ 0.018 MHz
$c_{\mathrm{lo}}$ | 6118 $\textrm{m}\,\textrm{s}^{-1}$ |  | $c_{\mathrm{tr}}$ | 3623 $\textrm{m}\,\textrm{s}^{-1}$

Once a solution to $\bm{u}$ and $\varphi$ has been obtained, the electrical
impedance $Z$ can be computed in the model as, the ratio of the total voltage
drop $V_{0}$ and the integral over the electrical current density $\bm{J}$
through the surface $\partial\Omega_{\frac{1}{2}}$ one half of the split
electrode

$Z=\frac{V_{0}}{I}=\frac{V_{0}}{\int_{\partial\Omega_{\frac{1}{2}}}\bm{J}\cdot\bm{n}\>\mathrm{d}a}.$
(33)

As there are no free charges in the system, the only contribution to the
current density $\bm{J}$ is the polarization current density
$\bm{n}\cdot\bm{J}=\bm{e}_{z}\cdot\bm{J}_{\mathrm{pol}}=\partial_{t}P_{z}=-\mathrm{i}\omega(D_{z}+\varepsilon_{0}\partial_{z}\varphi)$.
We compute $D_{z}$ from Eq. (13) and obtain, as in Ref. Skov  _et al._
(2019b),

$Z=\frac{\frac{\mathrm{i}}{\omega}V_{0}}{\int_{\partial\Omega_{\frac{1}{2}}}\\!\\!\big{[}e_{31}(\partial_{x}u_{x}\\!+\\!\partial_{y}u_{y})\\!+\\!e_{33}\partial_{z}u_{z}\\!+\\!(\varepsilon_{0}\\!-\\!\varepsilon_{33})\partial_{z}\varphi\big{]}\mathrm{d}a}.$
(34)

We have implemented the model of the AlN-thin-film-actuated glass-block
devices in the finite-element-method software Comsol Multiphysics 5.4 Com
following the procedure described in Ref. Skov  _et al._ (2019a). The coupled
field equations (2) for the elastic-solid displacement $\bm{u}$ and the
electric potential $\varphi$ are implemented together with the constitutive
equations (3) and (13) using the weak form interface “PDE Weak Form”. A COMSOL
script with a PDE-weak-form implementation of acoustofluidics is available as
supplemental material in Ref. Muller and Bruus (2015).

The system geometry described in Section II is meshed fairly coarse, because
we study the low end of the spectrum, mainly $f<3$ MHz, where the wavelengths
are larger than 1.1 mm. The mesh of the glass block is a free tetrahedral mesh
with a maximum element size 0.2 mm, whereas the AlN thin-film is modeled as a
swept, structured mesh with 2 mesh layers in the thickness direction and the
same element size along the surface as for the glass. We used quadratic
Lagrangian test functions, and used the mesh convergence test described in
Ref. Muller  _et al._ (2012), to verify that the numerical resolution was
satisfactory.

We used an adaptive step-size algorithm to compute the impedance spectra.
Close to a resonance peak, the step size went down to 15 Hz, but far from
resonance peaks, it increased to 1600 Hz. The resonance frequency was
determined as the position of the maximum of a Lorentzian fit containing the
11 closest points. A typical small 0.1-3-MHz spectrum contained 450-500
points, whereas a long 0.1-10-MHz spectrum contained 2500-2800 points. A short
and long spectrum took respectively around 0h:40min and 2h:50min to compute on
our HP-G4 workstation with a 3.7-GHz Intel Core i9 7960X-2.8 16-core processor
and with a memory of 128 GB RAM.

Figure 4:  (a) The resistive real part $\operatorname{Re}[Z]$ of the impedance
as a function of frequency $f$ measured for the 5 devices C2 - G2 of Fig. 3(d)
(gray curves $f^{\mathrm{exp}}$) and simulated numerically using the parameter
values from the literature, listed in Table 3, (red curve
$f^{\mathrm{sim,A}}$), and using the fitted parameter values of Table 4 (blue
curve $f^{\mathrm{sim,B}}$). (b) Similar plot of $\operatorname{Re}[Z]$ versus
$f$, but for the 7 devices A8 - G8 of Fig. 3(c).

##  VI Comparing simulation with experimental results

We use the numerical model of Section V to simulate electrical impedance
spectra $Z^{\mathrm{sim}}(f)$ through Eq. (34). In the following, we compare
these simulation results with experimental spectra $Z^{\mathrm{exp}}(f)$
described in Section IV.

Two examples of such a comparison in the frequency range from 0.1 to 3 MHz are
seen in Fig. 4(a) and (b) for device geometry 2 and 8, respectively, where we
are plotting the resistive part $\operatorname{Re}[Z]$ of the electrical
impedance as an alternative way to reveal the resonance peaks, instead of
using the phase spectra $\theta(f)$ as Fig. 3(c) and (d). The single red curve
in each plot is the simulation results obtained with the parameter values
listed in Table 3, whereas the family of gray curves is the experimental
result for all available devices with the geometry. It is seen that the
experimental impedance spectra are reproduced fairly well by the simulation
regarding the position of the resonance peaks and their relative heights. We
note, that similar to Fig. 3, the off-resonance level of the experimental
curves (gray) for the real part $\operatorname{Re}[Z]$ coincide for the
devices C2-G2 in Fig. 4(a) that have been discharged, while they are displaced
for the devices A8-G8 in Fig. 4(b) that have not been discharged. Moreover, we
note that in Fig. 4(a), we have fitted the off-resonance level in the
numerical model for $Z$ in Eq. (34) using $\Gamma_{\varepsilon}(f)$ given in
Table 3.

The agreement between experiment and simulation is quantified in Table 2,
where it is seen that for the 173 resonance peaks in our study, the relative
deviation $\Delta^{\mathrm{A}}_{k,n}$ of each simulated resonance frequency
$f^{\mathrm{sim,A}}_{k,n}$ deviates only a few percent from the experimental
value $f^{\mathrm{exp}}_{k,n}$. The mean value and the standard deviation of
the mean value of all relative deviations $\Delta^{\mathrm{A}}_{k,n}$ is
$(-4.2\pm 0.04)$%, which is a fairly good result. However, the simulation
results show a statistical significant systematic underestimate of the
resonance frequency. This is also clear from the insets in Fig. 4(a) showing a
zoom-in on four resonance peaks: In all cases, the measured (gray) peaks of
the seven nominally identical devices fall on top of each other, with
variations much smaller than the peak width, whereas in contrast, the
simulated (red) peaks are shifted downward in frequency by a few peak widths.

This systematic shift of the simulated resonance peaks leads us to
investigate, if a fitting of the two most significant but not so accurately
determined material parameters, Young’s modulus $E$ and Poisson’s ratio $s$ of
the SLS float-glass substrate, could lead to better numerical agreement with
the experiment. Indeed, by using a simple stepping procedure, we find that
such an improvement is obtained by the substitution

$E=72~{}\textrm{GPa}\xrightarrow[\text{fit}]{}75.5~{}\textrm{GPa},\qquad
s=0.23\xrightarrow[\text{fit}]{}0.189.$ (35)

This procedure extends the methods for obtaining the electromechanical
coefficients of piezoelectric transducers Lahmer  _et al._ (2008); Kiyono  _et
al._ (2016) to determining the elastic coefficients of any elastic solid
driven by vanishingly small thin-film transducers. The parameter values
obtained for the best fitting to the experimental data are summarized in Table
4, and the resonance frequencies obtained from the simulation based on these
values are denoted $f^{\mathrm{sim,B}}_{k,n}$.

The improved agreement between simulation and experiment is seen qualitatively
in Fig. 4, where the blue curves representing $f^{\mathrm{sim,B}}_{k,n}$ more
accurately reproduce the gray experimental impedance spectra. The improved
agreement is also evident from the insets in Fig. 4(a), where the blue curve
$f^{\mathrm{sim,B}}_{k,n}$ now falls on top of the gray experimental curves
within a fraction of the peak width. Quantitatively, the improved agreement
between experiment and simulation is seen in Table 2, where the relative
deviation $\Delta^{\mathrm{B}}_{k,n}$ of each simulated resonance frequency
$f^{\mathrm{sim,B}}_{k,n}$ deviates less than one percent from the
experimental value $f^{\mathrm{exp}}_{k,n}$. The mean value and the standard
deviation of the mean value of all relative deviations is
$\Delta^{\mathrm{B}}_{k,n}=(-0.03\pm 0.04)$%, or nearly 140 times better than
the first result $(-4.2\pm 0.04)$%.

Table 4:  Fitted values for Young’s modulus $E^{\prime}$, Poisson’s ratio $s^{\prime}$, and the damping coefficient $\Gamma_{\mathrm{sl}}$ of SLS float glass for improving the agreement between experiments and simulations. The values of $C^{\prime}_{11}$, $C^{\prime}_{12}$, $C^{\prime}_{44}$, $c_{\mathrm{tr}}^{\prime}$, and $c_{\mathrm{lo}}^{\prime}$ are functions of $E^{\prime}$ and $s^{\prime}$ Karlsen and Bruus (2015). Parameter | Value |  | Parameter | Value
---|---|---|---|---
$E^{\prime}$ | 75.5 GPa |  | $s^{\prime}$ | 0.189
$C^{\prime}_{11}$ | 82.8 GPa |  | $C^{\prime}_{44}$ | 31.8 GPa
$C^{\prime}_{12}$ | 19.2 GPa |  | $\Gamma_{\mathrm{sl}}(f)$ |  $f^{-1}$ 0.018 MHz
$c_{\mathrm{lo}}^{\prime}$ | 5731 $\textrm{m}\,\textrm{s}^{-1}$ |  | $c_{\mathrm{tr}}^{\prime}$ | 3549 $\textrm{m}\,\textrm{s}^{-1}$

The improved numerical agreement with experiments using the fitted values of
Young’s modulus and Poisson’s ratio for the glass substrate suggests that our
numerical fitting simulation algorithm may be used as an improved _in-situ_
determination of acoustic parameters in the 0.1-10 MHz range, and possibly for
higher frequencies. At the present stage we do not know, if the relative
deviation of roughly 4% between the literature values in Table 3 and the
fitted values in Table 4 is within the uncertainty generated by sample-to-
sample variations in SLS float glass, or if is a direct result of the
fabrication process involving the deposition of the AlN-thin-film transducers
at elevated temperatures. The improved simulation-to-experiment agreement also
suggest that the impedance measurements coupled with simulations is a
promising way to investigate the resonance modes, since there is little
ambiguity in which peaks are which. Also from the measurements shown in Fig.
4, it is clear that by combining simulation with more varied device geometries
it will be possible to predict particularly desired resonance frequencies with
good accuracy and consequently to design appropriate device geometries.

##  VII Conclusion

We have presented fabrication of 570-µm-thick, millimetric-sized SLS float-
glass blocks with a 1-µm-thick AlN-thin-film piezoelectric transducer
sandwiched between thin metallic electrodes and deposited on the top surface.
The electromechanical properties of the devices were characterized by
electrical impedance measurements in the frequency range from 0.1 to 10 MHz
with a peak-to-peak voltage of 0.5 V applied to the electrodes.

We have measured the electrical impedance spectra of 35 different device
geometries, all of width 2 mm, but with 9 different lengths ranging from 2 to
6 mm and with 2-7 copies of each individual geometry. Each impedance spectrum
exhibits many resonance peaks, of which we carefully measured the 5 most
prominent ones in each spectrum. Comparing the resulting 173 experimental
resonance frequencies with the result of a finite-element method simulation,
using material parameter values from the manufacturer, results in an average
relative deviation of the 173 simulated resonance frequencies from the
experimental ones of $(-4.2\pm 0.04)$%.

We have shown that by using Young’s modulus and Poisson’s ratio of the SLS
float glass as fitting parameters in the simulation, we could reduce the
relative deviation of the 173 simulated resonance frequencies from the
experimental ones by a factor of nearly 140 to $(-0.03\pm 0.04)$%. This result
suggests that our numerical fitting procedure could lead to an improved method
for _in-situ_ determination of the acoustic parameters at ultrasound
frequencies of any elastic solid onto which a thin-film transducer can be
attached.

## Acknowledgements

This work was supported by the BioWings project funded by the European Union’s
Horizon 2020 Future and Emerging Technologies (FET) programme, grant No.
801267.

## References

  * Ruby (2017) R. Ruby, The ’how & why’ a deceptively simple acoustic resonator became the basis of a multi-billion dollar industry, in  _2017 IEEE 30th International Conference on Micro Electro Mechanical Systems (MEMS)_ (2017) pp. 308–313.
  * Dubois and Muralt (1999) M. A. Dubois and P. Muralt, Properties of aluminum nitride thin films for piezoelectric transducers and microwave filter applications, Appl. Phys. Lett. 74, 3032 (1999).
  * Piazza  _et al._ (2006) G. Piazza, P. J. Stephanou, and A. P. Pisano, Piezoelectric aluminum nitride vibrating contour-mode MEMS resonators, J. Microelectromech. S. 15, 1406 (2006).
  * Zaghloul and Piazza (2014) U. Zaghloul and G. Piazza, Highly scalable nems relays with stress-tuned switching voltage using piezoelectric buckling actuators, IEEE T. Electron. Dev. 61, 3520 (2014).
  * Sinha  _et al._ (2009) N. Sinha, G. E. Wabiszewski, R. Mahameed, V. V. Felmetsger, S. M. Tanner, R. W. Carpick, and G. Piazza, Piezoelectric aluminum nitride nanoelectromechanical actuators, Appl. Phys. Lett. 95, 053106 (2009).
  * De Pastina  _et al._ (2018) A. De Pastina, D. Maillard, and L. Villanueva, Fabrication of suspended microchannel resonators with integrated piezoelectric transduction, Microelectron. Eng. 192, 83 (2018).
  * Olsson, III  _et al._ (2009) R. H. Olsson, III, K. E. Wojciechowski, M. S. Baker, M. R. Tuck, and J. G. Fleming, Post-CMOS-compatible aluminum nitride resonant MEMS accelerometers, J. Microelectromech. S. 18, 671 (2009).
  * Zhang  _et al._ (2006) H. Zhang, W. Pang, H. Yu, and E. S. Kim, High-tone bulk acoustic resonators on sapphire, crystal quartz, fused silica, and silicon substrates, J. Appl. Phys. 99, 124911 (2006).
  * Masson  _et al._ (2007) J. Masson, G. Martin, R. Boudot, Y. Gruson, S. Ballandras, A. Artieda, P. Muralt, B. Belgacem, and L. Chomeloux, On the dispersive behaviour of AlN/Si high overtone bulk acoustic resonators, in  _2007 IEEE International Frequency Control Symposium Joint with the 21st European Frequency and Time Forum_ (2007) pp. 741–744.
  * Fujikura  _et al._ (2008) T. Fujikura, O. Matsuda, D. M. Profunser, O. B. Wright, J. Masson, and S. Ballandras, Real-time imaging of acoustic waves on a bulk acoustic resonator, Appl. Phys. Lett. 93, 261101 (2008).
  * Muralt (2008) P. Muralt, Recent progress in materials issues for piezoelectric MEMS, J. Am. C. Soc. 91, 1385 (2008).
  * Howell  _et al._ (2019) K. M. Howell, W. Bashir, A. De Pastina, R. Matloub, P. Muralt, and L. G. Villanueva, Effect of AlN seed layer on crystallographic characterization of piezoelectric AlN, J. Vac. Sci. Technol. A 37, 021504 (2019).
  * Matloub  _et al._ (2011) R. Matloub, A. Artieda, C. Sandu, E. Milyutin, and P. Muralt, Electromechanical properties of al0.9sc0.1n thin films evaluated at 2.5 ghz film bulk acoustic resonators, Appl. Phys. Lett. 99, 092903 (2011).
  * Matloub  _et al._ (2013) R. Matloub, M. Hadad, A. Mazzalai, N. Chidambaram, G. Moulard, C. S. Sandu, T. Metzger, and P. Muralt, Piezoelectric Al1-xScxN thin films: A semiconductor compatible solution for mechanical energy harvesting and sensors, Appl. Phys. Lett. 102, 152903 (2013).
  * Fu  _et al._ (2017) Y. Q. Fu, J. K. Luo, N. T. Nguyen, A. J. Walton, A. J. Flewitt, X. T. Zu, Y. Li, G. McHale, A. Matthews, E. Iborra, H. Du, and W. I. Milne, Advances in piezoelectric thin films for acoustic biosensors, acoustofluidics and lab-on-chip applications, Prog. Mater. Sci. 89, 31 (2017).
  * Lee  _et al._ (2007) J.-B. Lee, D.-H. Cho, D.-Y. Kim, C.-K. Park, and J.-S. Park, Relationships between material properties of piezo-electric thin films and device characteristics of film bulk acoustic resonators, Thin Solid Films 516, 475 (2007).
  * Dubois and Muralt (2001) M. Dubois and P. Muralt, Stress and piezoelectric properties of aluminum nitride thin films deposited onto metal electrodes by pulsed direct current reactive sputtering, J. Appl. Phys 89, 6389 (2001).
  * Lee  _et al._ (2002) C. K. Lee, F. Placido, S. Cochran, and K. J. Kirk, Growth of sputtered aln thin film on glass in room temperature, in  _2002 IEEE Ultrasonics Symposium, 2002. Proceedings._, Vol. 2 (2002) pp. 1119–1122 vol.2.
  * (19) _Float glass substrates_ , [soda-lime-silicate glass: SiO2 (70.8%), Na2O (13.9%), CaO (8.4%), MgO (4.4%), Al2O3 (1.5%), K2O (0.4%), SO3 (0.3%), Fe2O3 (0.08%)], SCHOTT, SCHOTT Suisse SA, 2 Rue Galilée, Yverdon, Switzerland, https://www.epfl.ch/research/facilities/cmi/wp-content/uploads/2020/05/float_glass.pdf, accessed 16 November 2020.
  * Radovic  _et al._ (2004) M. Radovic, E. Lara-Curzio, and L. Riester, Comparison of different experimental techniques for determination of elastic properties of solids, Mat Sci Eng A-Struct 368, 56 (2004).
  * (21) A. Cachiaras, L. Gilde, J. Swab, P. Patel, and G. Quinn,  _Soda-Lime-Silicate Float Glass: A Property Comparison_ , US Army Research Laboratory, ARL, Aberdeen Proving Ground, MD 21005-5069, USA, https://www.researchgate.net/publication/320871872, accessed 16 November 2020.
  * Skov  _et al._ (2019a) N. R. Skov, J. S. Bach, B. G. Winckelmann, and H. Bruus, 3D modeling of acoustofluidics in a liquid-filled cavity including streaming, viscous boundary layers, surrounding solids, and a piezoelectric transducer, AIMS Mathematics 4, 99 (2019a).
  * Caro  _et al._ (2015) M. A. Caro, S. Zhang, T. Riekkinen, M. Ylilammi, M. A. Moram, O. Lopez-Acevedo, J. Molarius, and T. Laurila, Piezoelectric coefficients and spontaneous polarization of ScAlN, J. Phys.-Condens. Mat. 27, 245901 (2015).
  * Skov  _et al._ (2019b) N. R. Skov, P. Sehgal, B. J. Kirby, and H. Bruus, Three-dimensional numerical modeling of surface-acoustic-wave devices: Acoustophoresis of micro- and nanoparticles including streaming, Phys. Rev. Applied 12, 044028 (2019b).
  * (25) COMSOL Multiphysics 5.4, http://www.comsol.com (2018).
  * Muller and Bruus (2015) P. B. Muller and H. Bruus, Theoretical study of time-dependent, ultrasound-induced acoustic streaming in microchannels, Phys. Rev. E 92, 063018 (2015).
  * Muller  _et al._ (2012) P. B. Muller, R. Barnkob, M. J. H. Jensen, and H. Bruus, A numerical study of microparticle acoustophoresis driven by acoustic radiation forces and streaming-induced drag forces, Lab Chip 12, 4617 (2012).
  * Lahmer  _et al._ (2008) T. Lahmer, M. Kaltenbacher, B. Kaltenbacher, R. Lerch, and E. Leder, Fem-based determination of real and complex elastic, dielectric, and piezoelectric moduli in piezoceramic materials, IEEE Trans. Ultrason. Ferroelec. Freq. Contr. 55, 465 (2008).
  * Kiyono  _et al._ (2016) C. Y. Kiyono, N. Perez, and E. C. N. Silva, Determination of full piezoelectric complex parameters using gradient-based optimization algorithm, Smart Mater. Struct. 25, 025019 (2016).
  * Karlsen and Bruus (2015) J. T. Karlsen and H. Bruus, Forces acting on a small particle in an acoustical field in a thermoviscous fluid, Phys. Rev. E 92, 043010 (2015).


# Calculating entropy by Ma’s method for a system of $k$-mers on a square
lattice

Denis V Kokosinskii1 Mikhail V Ulyanov1,2 1 Computational Mathematics and
Cybernetics, M V Lomonosov Moscow State University, Moscow 119991, Russia  2 V
A Trapeznikov Institute of Control Sciences of RAS, Moscow 117997, Russia
deniskokossskappa<EMAIL_ADDRESS>
###### Abstract

Random walk of rectangular particles on a square lattice leads to a pattern
formation when the hard-core interaction between the particles is assumed. To
estimate changes in the entropy during this random walk, we propose a
modification of Ma’s method. The 2D sliding window technique was used to
divide a system under consideration into subsystems. We used Ma’s
“coincidence” method to estimate the total number of possible states for such
subsystems. In this study, the accuracy of Ma’s method is studied in a simple
combinatory model, both experimentally and theoretically. We determine which
definition of “coincidence” for this scheme leads to greater accuracy. Ma’s
estimate of the number of possible states for a system of $k$-mers correlates
well with the estimate obtained using a “naive” method.

##  1 Introduction

Boltzman’s entropy $S$ is defined as $S=k_{B}\ln(N),$ where $N$ is the number
of possible microstates, corresponding to the system’s macrostate, and $k_{B}$
is the Boltzmann’s constant. Boltzmann’s entropy is an important
characteristic of any evaluating system. Calculating Boltzmann’s entropy
directly as the logarithm of the total number of microstates for a current
macrostate is difficult for large systems. For complex, dense systems, direct
calculation of the Boltamzann’s entropy may require a great deal of
computation. The obvious solution to this problem is to obtain only an
estimate for the entropy, as this is faster to calculate. Therefore,
approximations of entropy, that are easier to compute, can help. Different
ways to estimate entropy were proposed [1, 2, 4, 5].

In our study, we propose a modification of Ma’s method [1], apply it to a
system of $k$-mers, and investigate its accuracy. In sec:modelandmethod, the
system under consideration is described, the modification of Ma’s method, and
the accuracy of Ma’s method is investigated in a simple combinatory model. In
sec:application, we apply Ma’s method to a system of $k$-mers. sec:concl
summarizes the main results.

##  2 Model and Methods

###  2.1 System under consideration

A square lattice of size $256$ by $256$ cells is used as a substrate. At
moment 0, elongated rectangular particles ($k$-mers) of size $1\times k$ or
$k\times 1$ lattice cells are placed onto the substrate using random
sequential adsorption (RSA) [6]. The RSA process [6] stops when no further
$k$-mer can be added to the lattice. A Monte Carlo simulation of the random
walk of $k$-mers is then begun. We shall, here give only a brief description
of this process, although a full description can be found at [7]. One Monte
Carlo step is $M$ attempts to shift a randomly selected $k$-mer in a random
direction by 1 lattice cell, where $M$ is the total number of $k$-mers on the
lattice. A $k$-mer can shift either along its’ orientation or perpendicular to
it, rotation being forbidden. The $k$-mers are “blind”, they are not aware of
the positions of other $k$-mers. If the position after the shift would cover a
position already occupied by other $k$-mers, the shift is not made and the
shift is attempted with a different $k$-mer. To avoid boundary effects, the
lattice is considered to have toroidal boundary conditions. Thus, a $k$-mer
moving beyond the right boundary of the lattice will appear on the left.

The evolution of the system of $k$-mers on the lattice is illustrated in
fig:system. Seemingly, this isolated system evolves from chaos to order.
Random walk of rectangular particles of two mutually perpendicular
orientations on a square lattice with periodic boundary conditions
demonstrated a pattern formation. The system tends to self-organize, as shown
in [8]. This pattern formation is supposed to be entropy-driven. Since the
complexity of counting the possible microstates grows exponentially with the
number of $k$-mers in the system, any direct computation of entropy changes is
hardly reasonable. An estimation of the entropy seems promising.

The complexity of the direct calculation of the Boltzmann’s entropy of this
system depends exponentially on the total number of $k$-mers. A “naive”
estimate of entropy at any given step can be calculated as

$S=\sum_{p}\ln(\text{DoF}(p)),$

where $p$ is any $k$-mer, and $\text{DoF}(p)$ is the number of degrees of
freedom of the particle $p$, in other words, the number of directions, in
which $p$ can shift. Boltzmann’s constant is omitted here and subsequently,
because we are interested in the relative values of entropy at the different
steps. This estimate does not take into consideration potential blocking of
$k$-mers shifts after some previous $k$-mers have moved.

Figure 1: Temporal evolution of the system. (a) $t_{MC}=0$, initial jammed
state, (b) $t_{MC}=100$, end of fluidization, (c) $t_{MC}=10^{4}$, labyrinth
patterns, (d) $t_{MC}=10^{5}$, labyrinth patterns are transforming into
stripes, (e) $t_{MC}=10^{6}$, final stage of stripe formation, (f)
$t_{MC}=10^{7}$, steady state in the form of diagonal stripes, (g) final
“yin–yang” pattern on a torus (the same pattern as in (f)); for clarity, the
major radius of the torus is significantly exaggerated. Illustration adopted
from [8].

###  2.2 Ma’s method in a simple combinatory model

Ma proposed estimating the total number of system states by counting the
number of “coincidences”, i.e. paired microstates. We aim to investigate this
method in a simple combinatory model. Let $N$ be a total number of boxes, $n$
be the number of balls thrown at these $N$ boxes with uniform distribution,
and $N_{p}$ be the number of resulting pairs (coincidences) in boxes. The
index of a box where the ball is thrown can be interpreted as this ball’s
microstate, so the total number of registered microstates would equal $n$. $N$
would then be the total number of possible microstates of the system. Suppose
$n$ and $N_{p}$ are given and we would like to determine $N$. According to Ma,
$N$ can be estimated as

$N=\frac{n(n-1)}{2N_{p}}.$ (1)

A question arises of how to correctly count the pairs. We can count pairs as
the number of boxes with exactly two balls ($N_{p1}$) or as the number of
boxes with two or more balls ($N_{p2}$). Another approach is to count all the
possible pairs in each box and then sum these values for all boxes ($N_{p3}$).
If $k$ balls are in a box, there are $\binom{k}{2}$ possible pairs of them. It
is important to notice, that in the case of $n\ll N$, all three approaches
produce close results. The probability of exactly $k$ balls in any box can be
calculated as

$p_{k}=\binom{n}{k}{\left(\frac{1}{N}\right)}^{k}{\left(1-\frac{1}{N}\right)}^{n-k},$
(2)

so

$\frac{p_{3}}{p_{2}}=\frac{n-2}{3N}\left(1-\frac{1}{N}\right).$

In the case $n\ll N$, $p_{3}\ll p_{2}$.

The results of experiments for estimating $N$ using different pair counting
strategies are presented in fig:paircounts. It is noticeable that method 3
produces good estimation of $N$ even when $n\sim N$. The errors of the other
methods appear to grow as $n$ increases. Let us examine this method in terms
of probability theory.

Figure 2: $N$ estimation using equation eq:ma with different approaches to
pair counting. Number of boxes $N=10000$. Estimate averaged over $1000$
experiments for each value of $n$.

Let $Z_{i}$ be a random variable, $Z_{i}=\binom{k}{2}$, when $k$ balls are in
a given box $i$. Using eq:pk, we can calculate the mathematical expectation of
$Z_{i}$:

$E(Z_{i})=\sum_{j=2}^{N}{\binom{j}{2}p_{j}}.$ (3)

We can rewrite $p_{k}$ as

$p_{k}=\binom{n}{k}{\left(\frac{1}{N}\right)}^{k}{\left(1-\frac{1}{N}\right)}^{n-k}=\left[\,\frac{\binom{n}{k}}{\binom{n}{2}}{\left(\frac{1}{N}\right)}^{k-2}{\left(1-\frac{1}{N}\right)}^{-k+2}\right]\,\cdot
p_{2}$

or

$p_{k}=\left[\,\frac{\binom{n}{k}}{\binom{n}{2}}{\left(\frac{1}{N-1}\right)}^{k-2}\right]\,p_{2}.$
(4)

Substituting $p_{k}$ in eq:EZi with eq:pkfromp2 we get

$E(Z_{i})=\sum_{j=2}^{N}{\left[\,\binom{j}{2}\frac{\binom{n}{j}}{\binom{n}{2}}{\left(\frac{1}{N-1}\right)}^{j-2}p_{2}\right]\,}.$

We then use the properties of binomial coefficients’ and move the common
factor out

$E(Z_{i})=p_{2}\sum_{j=2}^{N}{\left[\,{\binom{n-2}{j-2}}{(\frac{1}{N-1})}^{j-2}\right]\,}$

and finally apply the binomial formula

$E(Z_{i})=p_{2}{\left(1+\frac{1}{N-1}\right)}^{n-2}.$ (5)

Now, if we introduce a random variable $Z=\sum_{i=1}^{N}{Z_{i}}$, it takes the
values of $N_{p3}$. The mathematical expectation of $Z$ is

$E(Z)=\sum_{i=1}^{N}{E(Z_{i})}=Np_{2}{\left(\frac{N}{N-1}\right)}^{n-2}.$

Using representation eq:pk, we get:

$E(Z)=N\left[\,\binom{n}{2}{\left(\frac{1}{N}\right)}^{2}{\left(\frac{N-1}{N}\right)}^{n-2}\right]\,{\left(\frac{N}{N-1}\right)}^{n-2}=\binom{n}{2}\frac{1}{N}=\frac{n(n-1)}{2N}.$

We can now get an estimate of N:

$N=\frac{n(n-1)}{2E(Z)}$

We have proved that the formula $N=\frac{n(n-1)}{2N_{p3}}$ provides an
unbiased estimate of $N$ regardless of $n$.

##  3 Results: Application of Ma’s method to a system of $k$-mers

To apply Ma’s method, we first need to define a microstate for the system of
$k$-mers on the lattice. The lattice can be subdivided into small non-
overlapping squares. We can define a microstate as the system state inside one
of these squares. However, in this case, the results of computer experiments
show, that no pairs (equal states of the system in different squares) can be
found at any step of a Monte Carlo simulation of random walk. The existence of
at least one pair is required to apply Ma’s method.

We can use another definition of a microstate, based on [4]. The 2D sliding
window approach was used in this article to investigate the various structural
properties of the same system and process that we are studying. At any given
Monte Carlo step, a small square window of size $L$ by $L$ lattice cells moves
over all the positions on the system. There is a total of $256\times 256$
positions of the window, because the boundary conditions are toroidal. We can
count identical states of the window in different positions and use Ma’s
method to estimate the total number of possible window states. Here, and
subsequently, we count the number of pairs according to the best pair counting
approach from sec:modelandmethod. The two definitions of a microstate are
illustrated in fig:windowdef.

Figure 3: Two definitions of a microstate for a system of $k$-mers. (a) non-
overlapping squares; (b) 2D sliding window. A “paired” microstate is shown in
(b).

We need to choose the window size $L$. There are two important considerations
for $L$. First, there must be at least one pair of window states at any Monte
Carlo step. Second, the bigger the size of the window, the better it
represents the state of the whole system. The total number of pairs of window
states during the Monte Carlo process for various window sizes $L$ is
presented in fig:windowsize. A window size $L=7$ satisfies the requirements.

Figure 4: The total number of pairs for different window sizes $L$ during the
Monte Carlo simulation of random walk for the system of $k$-mers.

##  4 Conclusion

Our study demonstrates, that, in the system under consideration, the
rearrangement of particles and pattern-formation are entropy-driven, i.e., in
its initial state, a particle has less possibilities to change its location in
comparison with any succeeding state [9]. The Ma’s entropy estimate for a 2D
sliding window during a Monte Carlo simulation of diffusion of a system of
$k$-mers on a lattice is presented in fig:results. The “naive” method,
introduced in sec:intro, is also presented as the baseline. Both estimates
grow, despite the fact, that the system self-organizes. Worthy of note, that
Ma’s estimate is, in practice, several magnitudes smaller, than the “naive”
estimate. As the positions of the sliding window overlap, so it may not be
correct to assume, that its’ states are independent and uniformly distributed.
However, Ma’s method, if using the correct definition of a pair from
sec:modelandmethod, will give close estimates if these conditions are
satisfied. It might be worth testing to estimate the entropy in systems where
the microstates are independent and uniformly distributed.

Figure 5: Entropy estimates during the Monte Carlo simulation of the diffusion
of $k$-mers on a lattice. Entropy is normalized to the entropy at step 0.

We acknowledge funding from the Russian Foundation for Basic Research, Project
No. 18-07-00343. We are thankful to Yuri Tarasevich for his stimulating
discussions.

## References

## References

  * [1] Ma S k 1981 J. Stat. Phys. 26 221
  * [2] Romá F, Ramirez-Pastor A J, Riccardo J L 2001 J. Chem. Phys. 114 10932
  * [3] Vogel E E, Saravia G, Bachmann F, Fierro B, Fischer J 2009 Physica A 388 4075
  * [4] Ulyanov M V, Smetanin Y G, Shulga M M, Eserkepov A V and Tarasevich Y Y 2018 J. Phys. Conf. Ser. 1141 012137
  * [5] Vogel E E, Saravia G, Ramirez-Pastor A J, Pasinetti M 2020 Phys. Rev. E 101 022104
  * [6] Evans J W 1993 Rev. Mod. Phys. 65(4) 1281
  * [7] Lebovka N I, Tarasevich Y Y, Gigiberiya V A and Vygornitskii N V 2017 Phys. Rev. E 95(5) 052130
  * [8] Tarasevich Y Y, Laptev V V, Burmistrov A S and Lebovka N I 2017 J. Stat. Mech: Theory Exp. 2017 093203
  * [9] Frenkel D 1993 Physics World 6 24