# Mastering Complex Modes: A New Method for Real-Time Modal Identification of
Vibrating Systems

Satyam Panda
Department of Civil Engineering
Indian Institute of Technology Guwahati
<EMAIL_ADDRESS>
& Sanghamitra Das
Department of Mechanical Engineering
Indian Institute of Technology Guwahati
<EMAIL_ADDRESS>
& Basuraj Bhowmik
Department of Civil Engineering
Indian Institute of Technology (BHU), Varanasi
<EMAIL_ADDRESS>
& Budhaditya Hazra
Department of Civil Engineering
Indian Institute of Technology Guwahati
<EMAIL_ADDRESS>

###### Abstract

A novel algorithm for real-time modal identification in linear vibrating
systems with complex modes is introduced, utilizing a combination of first
order eigen-perturbation and second order separation techniques. In practical
settings, structures with complex modes are frequently encountered and their
presence often poses a challenge in accurately estimating the source signal in
real-time. The proposed methodology addresses this issue by incorporating the
right angle phase shift of the response in the sensor output and updating the
second order statistics of the complex response through first order eigen-
perturbation. Empirical evidence of the efficacy of the technique is
demonstrated through numerical case studies and validation using various
numerically modeled systems, as well as a standard ASCE-SHM benchmark problem
with complex modes, highlighting the capability of the proposed method to
achieve precise real-time modal property identification and online source
separation with a minimal number of initially required batch data.

_Keywords_ Complex modes $\cdot$ modal identification $\cdot$ blind source
separation $\cdot$ eigen perturbation $\cdot$ second order separation

##  1 Introduction

Vibration-based condition monitoring is a diagnostic methodology utilized to
detect anomalies in machinery [1, 2], aerodynamics models [3], structures,
such as buildings and bridges [4, 5], to name a few. The vibrational response
of these systems under external excitation can vary depending on their
structural and energy dissipation properties [6]. Classical viscous damping
models – such as Rayleigh damping – are often employed to represent the
damping forces of the system for dissipative forces. Although these
assumptions are accurate for simpler models, real systems often contain many
different materials and are combined together in different ways to form
complex spatial geometries, these assumptions are sometimes impractical and
constricting. Non-classical damped linear systems can exhibit unique complex
mode shapes [7] due to non-proportional damping matrices [8]. Identification
of these intricate modes is crucial for accurate characterization of the
system and can be achieved through modal identification techniques, a well-
researched area within vibration-based condition monitoring [9, 10, 11, 12].
However, traditional algorithms are not suited for dynamic objects, like
windmills and bridges, due to their offline nature. Recent studies have
introduced first order eigen perturbation (FOEP) techniques as a means to
continuously update the eigenspace data and provide real-time vibration modes
[13]. Although these methods perform well for real valued modes [14], they are
not applicable to complex valued modes. In light of this, the present study
proposes a novel method for real-time blind source modal identification that
combines eigen perturbation and second order source separation in a recursive
framework. This method enables accurate estimation of complex modes and modal
response at each moment in time.

The past decade has seen widespread application of eigen perturbation theory
[13] in the field of vibrating systems, particularly for the detection of
damage instances [15], modal identification [13], enhancement of system
response, tuning of tune mass dampers (single and multi mode), and filtering
of non-stationary signals [16]. These FOEP techniques are data-driven and
exhibit reduced dependence on baseline parameters, making them a valuable tool
for real-time condition monitoring. Despite their efficacy for the
identification of real valued mode shapes, FOEP techniques exhibit limitations
in the estimation of complex valued modes. In the context of infrastructure
asset management, these methods can aid in the effective monitoring and
maintenance of critical assets. Modal analysis applications continue to rely
on real valued mode shapes for tuning, updating, and damage identification.
However, real-time algorithms that update the eigenspace at each time instant
can lead to an accumulation of errors due to the lack of consideration of the
complex phasing of the system during updates. This ignorance of the complex
nature of vibrational modes can result in inaccurate results for in-situ
problems that require high precision. It is important to note that complex
modes can exist in structures without damping – such as those with closely
spaced modes [17], subjected to aerodynamic loads [18], or characterized by
asymmetry [19]. Ignoring these complex modes in modal analysis can negatively
impact the accuracy of results.

The traditional real-time FOEP algorithms update the eigenspace of the system
through numerical simulation or sensor-acquired data. The eigenvector is then
checked for permutation ambiguity to determine the modes of the system.
However, many systems possess inherently complex modes that are overlooked by
these algorithms. To overcome this limitation, the proposed algorithm involves
taking two simultaneous samples and performing a Hilbert transform to obtain a
90∘ phase shift, resulting in a complex response. The FOEP algorithm is then
used to update the complex eigenspace with the updated eigenspace serving as
the basis for estimating the whitening matrix and updating the delayed
covariance matrix. This novel approach offers improved accuracy in
characterizing complex modes, which has significant commercial and sectorial
benefits in the field of civil engineering, including improved infrastructure
asset management, enhanced vibration-based damage identification, and advanced
tuning and filtering techniques. The proposed real-time algorithm employs the
principles of eigen perturbation and second order source separation to
accurately estimate the complex-valued modes of ever-changing systems. This is
accomplished through the simultaneous acquisition of two samples, the
calculation of the Hilbert transform, and the updating of the eigenspace and
delayed covariance matrix using FOEP techniques. The resulting complex de-
mixing matrix separates the source signals, enabling the estimation of complex
modes in real-time. The development of this algorithm provides significant
benefits for the commercial and industrial sector, particularly in the field
of civil engineering, by enabling the accurate and efficient assessment of
infrastructure assets.

The proposed framework in this study constitutes a major advancement in the
field of real-time condition monitoring for civil and structural engineering
applications. Specifically, it offers the following key contributions:

  1. 1.

Development of a mathematically consistent algorithm for the real-time
estimation of complex modes in dynamic systems, which has been an understudied
area in the field of structural health monitoring and damage detection.

  2. 2.

Incorporation of second-order delayed statistics updates into traditional FOEP
techniques to account for the constantly changing complex phasing in systems
with dynamic modes.

  3. 3.

The proposed framework is applicable for both real and complex modal
identification in a variety of simulations and practical problems, providing a
more comprehensive solution for the blind source problems faced in real-life
structural health monitoring and damage detection applications.

Overall, this study represents a significant step towards the accurate
identification of complex modes in real-time, which is a crucial aspect in the
development of effective structural health monitoring and damage detection
systems.

The paper’s structure comprises of the following segments: Section 2
delineates a succinct exposition of the complex state theory within the
purview of structural dynamics and the conventional batch Blind Source
Separation (BSS) technique, namely the Second Order Blind Identification
(SOBI), utilized for the calculation of complex modes. Section 3 propounds a
mathematically rigorous formulation of the proposed methodology for
determining complex modes of general dynamic systems, incorporating an
algorithm and flowchart. Section 4 presents numerical simulations on selected
categories of linear systems to exemplify the application of the proposed
methodology in both simulation and real-world situations, validating its
efficacy against the Structural Health Monitoring (SHM) benchmark structure
outlined by the International Association for Structural Control and
Monitoring - American Society of Civil Engineers (IASC-ASCE). Lastly, Section
5 summarizes the work and provides insightful observations.

##  2 Background

This section explains the necessity of identification of the complex modes in
structural systems with a concise summary of blind modal identification
techniques.

###  2.1 Complex modes in structural dynamics

The concept of complex modes in the light of structural dynamics can be
understood by considering the generic dynamical system of the form: where,
mass, stiffness and damping matrices are represented by M, C and K with X as
the displacement vector and $\textbf{F}(t)$ as the input excitation. Since the
matrices M, C and K are diagonalized in the modal coordinate the modal
responses (${\bf{Q}}={[{q_{1}},{q_{2}},\ldots{q_{i}}\ldots{q_{n}}]^{T}}$) of
the Eqn. Eq. 2 can be written as:

${{\ddot{q}}_{i}}(t)+2{\zeta_{i}}{\omega_{n,i}}{{\dot{q}}_{i}}(t)+{\omega^{2}}_{n,i}{q_{i}}(t)={\frac{1}{m_{i}}}{f_{i}}(t)$
(1)

where, ${f_{i}}(t)={v_{i}}^{T}{\bf{F}}(t)$ is the modal force corresponding to
the ${v_{i}}$ mode. $\zeta_{i}$, ${\omega^{2}}_{n,i}$ and $m_{i}$ are the
modal damping, modal frequency and modal mass for $i^{th}$ mode. The solution
of the Eq. 1 is obtained by using Duhamel Integral:
${q_{i}}(t)=\int\limits_{0}^{\infty}{{f_{i}}(\tau){h_{i}}(t-\tau)}{\kern
1.0pt}d\tau$, where
${h_{i}}(t-\tau)=\frac{1}{{{m_{i}}{\omega_{d,i}}}}{e^{-{\zeta_{i}}{\omega_{n,i}}(t-\tau)}}\sin{\kern
1.0pt}{\kern 1.0pt}{\omega_{d,i}}(t-\tau)$ and ${\omega^{2}}_{d,i}$ is the
damped modal frequency. The solution matrix ${\bf{Q}}$ in the vectorial form
can be expressed as (using convolution property):

${\bf{M\ddot{X}}}(t)+{\bf{C\dot{X}}}(t)+{\bf{KX}}(t)={\bf{F}}(t)$ (2)
$\xi_{ij}(t)=P(x_{t}=i,x_{t+1}=j|y,v,w;\theta)={\frac{\alpha_{i}(t)a^{w_{t}}_{ij}\beta_{j}(t+1)b^{v_{t+1}}_{j}(y_{t+1})}{\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}(t)a^{w_{t}}_{ij}\beta_{j}(t+1)b^{v_{t+1}}_{j}(y_{t+1})}}$
(3)

${{\bf{Q}}_{i\times
1}}(t)=\int\limits_{0}^{\infty}{\underbrace{{{\left[{\begin{array}[]{*{20}{c}}{{h_{1}}(t-\tau)}&0&\cdots&0\\\
0&{{h_{2}}(t-\tau)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&{{h_{i}}(t-\tau)}\end{array}}\right]}_{i\times
i}}}_{{{\bf{H}}^{Q}}(t-\tau)}}\underbrace{{{\left\\{{\begin{array}[]{*{20}{c}}{{f_{1}}(\tau)}\\\
{{f_{2}}(\tau)}\\\ \vdots\\\ {{f_{i}}(\tau)}\end{array}}\right\\}}_{i\times
1}}}_{{{{F}}^{Q}}(\tau)}d\tau$ (4)

The solution of Eq. (2) can be written as ${\bf{X}}={\bf{VQ}}$, where X is the
measurement matrix of size $m\times N$ and Q is the corresponding modal
response of size $m\times N$ with $m$ as number of degrees of freedom and $N$
as the sampling size. V is a orthogonal transformation matrix of size $m\times
m$ yielding mode matrix such that the mode shape matrix are orthogonal to each
other with respect to the matrix M. In general, mode shape matrices of
structural dynamical systems, whether with constant system parameters or with
varying mass and stiffness properties, are often approximated as proportional
modes. However, a complex modal solution is necessary when addressing system
modifications with respect to factors such as damping or non-linearity. If
proportional mode approximations are used, errors can accumulate due to the
complex phasing resulting from the modified system, which the proportional
approximation overlooks.

Complex modes can arise in non-proportionally damped systems either for a
narrow range of system parameters or for systems with closely spaced modes.
Complex modes can also occur in aeroelastic systems under non-linear
aerodynamic loads. Furthermore, the presence of analysis and measurement
errors can also lead to complex modes [7]. Complex eigenvectors can also
result from an asymmetric velocity-dependent system matrix, even if there is
no damping and only a gyroscopic term present in the ${\bf{C}}$ matrix. This
complexity can also arise in axisymmetric structures, where modes are often
identical-eigenvalue pairs and hence, have a degree of degeneracy. It is
crucial for practitioners dealing with real-time modal identification problems
to utilize algorithms that can handle both proportional and complex modes to
avoid incorrect modal identification. This necessitates the effective use of
Blind Source Separation (BSS) algorithms with specific capabilities for real-
time mechanism and complex modal identification. To address these problems, a
succinct understanding of the BSS approach and its application to structural
dynamics is required, which is discussed in the following section.

###  2.2 Blind source separation (BSS)

Blind Source Separation (BSS) refers to the process of extracting source
signals from output sensor signals without any prior knowledge of the system.
In BSS, only the mixed output signals are available for separation in order to
retrieve the source signals. Numerous studies have been conducted to solve
modal identification problems using BSS, which can act as a useful inverse
methodology for retrieving the source signals of an unknown Multiple-
Input/Multiple-Output (MIMO) mixing system. The input signals are estimated
based on the output sensor signals. Several BSS techniques for estimating the
source signals include Independent Component Analysis (ICA) and its extensions
(Topographic ICA, Multidimensional ICA, Kernel ICA, Tree-dependent Component
Analysis, Subband Decomposition -ICA), Sparse Component Analysis (SCA), Sparse
PCA (SPCA), Non-negative Matrix Factorization (NMF), Smooth Component Analysis
(SmoCA), Parallel Factor Analysis (PARAFAC), Time-Frequency Component Analyzer
(TFCA), and Multichannel Blind Deconvolution (MBD) [20]. BSS finds wide
applications in many areas, such as array signal processing, seismic signal
processing, and blind equalization [21]. It assumes that the output data is a
linear combination of the modal coordinates, which are estimated from the
output signals only. This makes BSS a useful technique in the estimation of
modal responses from ambient data. By comparing the superposition of vibration
modes with the BSS model’s specification, it is possible to understand the
physical significance of the sources in system identification problems. The
mode shapes, natural frequencies, and damping factors of a vibrating structure
can be estimated using BSS, providing a numerical solution for the deflection
patterns of vibration when the system vibrates at one of its natural
frequencies. BSS has various applications in the field of signal processing
and engineering, including array signal processing, seismic signal processing,
and blind equalization. By assuming the output data is a linear combination of
modal coordinates and only using the output signals for estimation, linear
mixing models can be used to estimate modal responses from ambient data. These
models can then be used to determine the natural frequencies, damping factors,
and mode shapes of a vibrating structure, which provide insight into its
physical behavior during vibration. BSS plays an important role in system
identification problems by facilitating the comprehension of the physical
importance of sources.

###  2.3 Second Order Blind Identification (SOBI)

For system identification problems, the popular BSS technique known as Second
Order Blind Identification (SOBI) is frequently used [22].SOBI has a wide
variety of applications in modal identification of vibrating structures [17].
It has proven to be a viable substitute to produce accurate results for
output-only modal analysis of vibrating structures [23]. Independent component
analysis (ICA) is found to give better results for identification of mode
shapes as compared to Principal component analysis (PCA). The reason behind
this is the poor estimation performance and high computation time at higher
order statistics and damping ratios greater than 1 percent [22]. A suitable
alternative to estimate the mode shapes for second order statistics is the
second order blind identification (SOBI) which assumes uncorrelated sources
with variable spectral contents[24]. SOBI estimates real valued mode shapes
with real values mixing matrix. It uses the time information contained in the
signals in contrary to considering them as random variables in ICA. SOBI has
also shown to be a promising technique for time series dimension reduction.
SOBI is based on the premise of simultaneous diagonalization of two covariance
matrices ${\bf\hat{R}}_{\bf Y}\left(0\right)$ and ${\bf\hat{R}}_{\bf
Y}\left(p\right)$, defined by the relations as under:

$\left.\begin{array}[]{l}{\bf
R}_{Y}\left(p\right)=E\left\\{{{Y}\left(k\right){s}^{T}\left({k-p}\right)}\right\\}\\\
{\bf
R}_{Y}\left(0\right)=E\left\\{{{Y}\left(k\right){Y}^{T}\left(k\right)}\right\\}=A{\bf
R}_{Y}\left(0\right)A^{T}\\\ {\bf
R}_{Y}\left(p\right)=E\left\\{{{Y}\left(k\right){Y}^{T}\left({k-p}\right)}\right\\}\\\
\end{array}\right\\}$ (5)

for some non-zero time-lag $p$. The simultaneous diagonalization is performed
in three basic steps: whitening, orthogonalization, unitary transformation. A
framework for second order blind identification method can be summarized as
follows:

  1. 1.

Obtain output dataset, assign to ${{Y}_{0}}\,(t)$. Compute cross correlation
functions from the data, and assign to ${{Y}_{0}}\,(t)$.

${{R}_{Y\,}}(0)\,=\,\frac{1}{N-1}\,{{Y}^{T}}$ (6)

  2. 2.

Whiten the data to obtain whitened data, $\bar{Z}$. Retain the whitening
matrix and its inverse. Whitening is a linear transformation in which
${\bf\hat{R}}_{\bf Y}\left(0\right)=\left({{\raise
3.01385pt\hbox{$1$}\\!\mathord{\left/{\vphantom{1N}}\right.\kern-1.2pt}\\!\lower
3.01385pt\hbox{$N$}}}\right)\left({\sum\limits_{K=1}^{N}{\textbf{Y}\left(k\right)\textbf{Y}^{T}\left(k\right)}}\right)$
is first diagonalized using singular value decomposition that is accomplished
as ${\bf\hat{R}}_{\bf Y}\left(0\right)={\bf V}_{\bf Y}{\bf\Lambda}_{\bf Y}{\bf
V}_{\bf Y}^{\bf T}$. The whitened signal is given as:

${Y}\left(k\right)={\bf W}Y\left(k\right)={\bf\Lambda}_{Y}^{-\frac{1}{2}}{\bf
V}_{Y}^{\bf T}{Y}\left(k\right)$ (7)

where ${\bf W}$ is the whitening matrix.

  3. 3.

Apply joint approximate diagonilization (JAD) to the whitened data to obtain
the joint diagonalizer. Orthogonal transformation is applied to diagonalize
the matrix ${\bf\hat{R}}_{Z}\left(p\right)$ The Eigen value decomposition of
${\bf\hat{R}}_{Z}\left(p\right)$ has the form

${\bf\hat{R}}_{Z}\left(p\right)={\bf U}_{Z}{\bf\Sigma}_{Z}{\bf U}_{Z}^{T}.$
(8)

  4. 4.

Compute the de-mixing matrix, $\hat{\bf A}$, and the mixing matrix, ${\bf A}$.

${\bf\hat{A}}=\,{\bf U}_{Z}^{\dagger}\,{\bf W}$ (9)

If the diagonal matrix ${\bf\Sigma}_{Z}$ has distinct eigen values then the
mixing matrix can be estimated uniquely by the following equation:

${\bf A}=\textbf{W}^{+}{\bf U}_{Z}$ (10)

where, $(.)^{+}$ is the Moore-Penrose pseudoinverse of a matrix.

  5. 5.

Compute the modal response. $X\,=\,{{\bf A}^{T}}\,Y$

###  2.4 General eigen perturbation theory

Recently, advancements have been made in the field of real-time estimation of
vibrational modes through the implementation of eigen perturbation techniques.
The focus of these algorithms lies in the recursive update of the eigenspace,
providing an alternative method for updating the full rank data covariance
matrix. The proposed approach utilizes the framework of First-Order Eigen
Perturbation (FOEP) technique for the recurrent update of the complex
eigenspace, thereby updating the delayed whitening covariance matrix. As a
preface to the mathematical intricacies of the proposed algorithm, it is
imperative to have a comprehensive understanding of the generalized higher
order eigen perturbation theory and the derivation of FOEP from it. In the
context of a vibrating system, the symmetric eigen decomposition can be
mathematically represented as:

${\bf{R}}_{Y}{\bf{V}}={\bf{\Lambda}}{\bf{V}}$ (11)

where, ${{\bf{R}}_{X}}\in{\mathbb{R}^{n\times n}}$, and ${{\bf{\Lambda}}}$,
${{\bf{V}}}$ are eigenvalues and eigenvectors matrices, respectively. As the
eigenvector matrix is symmetric and eigenvalue matrix is diagonal, the
following orthogonality conditions will hold

$\begin{array}[]{l}{V}_{j}^{T}{{V}_{i}}={{\delta_{ij}}},\quad
V_{j}^{T}{\bf{R}}_{Y}{V_{i}}={\Lambda_{i}}{\delta_{ij}}\quad\forall\quad
i,j=1,2,3,\ldots,n\end{array}$ (12)

where ${\delta_{ij}}$ is the Kronecker delta. For a generalized $n^{th}$ order
perturbation, Eq. (11) can be written as a linear combination of $n$
perturbation terms added to the original matrix, given by the expression:

$\begin{array}[]{l}{{\bf{R}}_{Y}}=\left({\bf{V}}+\sum_{i=1}^{n}\delta^{i}{\bf{V}}\right)^{T}\left({{\lambda}+\sum_{i=1}^{n}\delta^{i}\lambda}\right)\left({{{\bf{V}}}+\sum_{i=1}^{n}\delta^{i}{\bf{V}}}\right)\end{array}$
(13)

This generalized eigen perturbation theory provides a method for updating the
data covariance matrix without requiring repeated eigen decomposition of the
data covariance matrix at each time step. However, the higher the perturbation
order involved in the expansion, the more accurate the algorithm becomes, but
also the more mathematically complex and computationally demanding. Hence,
real-time system monitoring algorithms generally adopt either First [13] or
Second Order perturbation terms [25]. Nevertheless, recent studies [14, 15,
16] have incorporated error correction mechanisms to mitigate the impact of
neglected higher-order perturbation terms and enhance the convergence and
damage detection capability of traditional FOEP techniques. However, in this
work, the focus is on the recursive whitening of signals using the eigenspace
update of complex data, therefore the use of higher-order perturbation or
error feedback mechanisms would increase computational complexity without
significantly improving the robustness of the algorithm. The proposed
framework combines the First Order Eigen Perturbation (FOEP) technique with
the concept of Second Order Blind Source Separation (SOBSS) in a recursive
manner, providing a mathematically consistent approach to the estimation of
complex modes, source signal extraction, and will be discussed in detail in
the following section.

##  3 Recursive modal identification in dynamical systems

The basic principle of the recursive modal identification of the systems with
complex mode shapes primarily premise on three primary segments operating
simultaneously in a recursive framework: Firstly whitening of the complex
response using RPCA to obtained the whitening matrix, its inverse and whitened
data for each sample, Secondly updating the lagged covariance matrix and
Finally estimation of unitary matrix from joint diagonalization of the updated
lagged covariance matrix and in turn estimating the mixing matrix. These
process require an initialization of some parameters which are obtained using
the batch algorithms for a few initial samples (ideally 100-500samples).

The incorporation of the phase shifted data obtained through the Hilbert
transformation enhances the accuracy of the real-time modal identification
algorithm in handling non-proportional damping, non-linearity, gyroscopic
effect, and closely placed modes, among others. The real-time feeding of the
sensor data and phase shifted data to the algorithm forms a complex response,
thereby allowing for more accurate identification of the modal response in an
online fashion. Thus the proposed algorithm supplements the recorded response
$Y(t)$ with the phase shifted data $Y_{90}(t)$ obtained through employing the
Hilbert transformation (HT) and taking its imaginary value as 90∘ phase
shifted response. The key point is that the response is fed to the algorithm
in real time to obtain the phase shifted data and is supplemented to the
sensor data for formation of the complex response as,

${\bar{Y}}(t)=Y(t)+\text{i}\ Y_{90}(t)$ (14)

This leads to the problem of finding the complex valued modal response through
recursive BSS techniques and can be written in the form of complex modal
response as,

$\left\\{Y_{k}+\text{i}\ Y_{90,k}\right\\}=\left[{\bf{A}}_{0,k}+\text{i}\
{\bf{A}}_{90,k}\right]\left\\{X_{k}+\text{i}\ X_{90,k}\right\\}$ (15)

After obtaining the complex response, the objective shifts towards the
estimation of the whitening matrix, where, the mathematically consistent
formulation of RPCA for the complex data is necessary. Towards understanding
the mathematical formulation of RPCA, consider the co-variance matrix of the
form: ${{\bf{C}}_{Y}}=\frac{1}{N}{\bf{Y}}{{\bf{Y}}^{{T}}}$. For any
multivariate data set the recursive estimation of the data co-variance matrix
(${{\bf{C}}_{k}}$) at time instant $k$ can be expressed in terms of the data
vector at $k^{th}$ instant i.e. ${Y_{k}}$ as:

${{\bf{C}}_{k}}=\frac{{k-1}}{k}{{\bf{C}}_{k-1}}+\frac{1}{k}\left\\{Y_{k}+\text{i}\
Y_{90,k}\right\\}\left\\{Y_{k}+\text{i}\ Y_{90,k}\right\\}^{T}$ (16)

where ${{\bf{C}}_{k-1}}$ is the co-variance estimate at $(k-1)^{th}$ instant.
The primary aim of RPCA is to estimate the update of co-variance matrix at
each time instant without actually performing EVD on the block co-variance
matrices thereby reducing the time and memory complexity. First order eigen
perturbation (FOEP) technique facilitates the estimation of eigenspace without
actually performing the EVD recursively. The individual data co-variance
estimates at $k^{th}$ instant can be expressed by its EVD as,
${{\bf{C}}_{k}}={{\bf{V}}_{k}}{\bf{\Sigma}}_{k}{{\bf{V}}_{k}}^{T}$ where
${\bf{V}}_{k}$ and ${\bf{\Sigma}}_{k}$ represents the eigenvector and
eigenvalue matrices at $k^{th}$ instant. For non-stationary processes the
drift in mean level can be accommodated in the estimate of co-variance update
by using mean shift in the data samples as follows:

$\displaystyle{{\bf{V}}_{k}}{\bf{\Sigma}}_{k}{{\bf{V}}_{k}}^{T}=\frac{{(k-1)}}{k}{{\bf{V}}_{k-1}}{\bf{\Sigma}}_{k-1}{{\bf{V}}_{k-1}}^{T}+\frac{1}{k}\left[{\left\\{Y_{k}+\text{i}\
Y_{90,k}\right\\}-{\mu_{k}}}\right]{\left[{\left\\{Y_{k}+\text{i}\
Y_{90,k}\right\\}-{\mu_{k}}}\right]^{T}}$ (17)

where, the recursive mean at $k^{th}$ instant is estimated as:
${\mu_{k}}=\frac{{k-1}}{k}{\mu_{k-1}}+\frac{1}{k}\left\\{Y_{k}+\text{i}\
Y_{90,k}\right\\}$. However, for structural systems the response data vectors
are generally evolve as a zero mean processes enabling the Eq. 17 to be
written without any mean shift. Further simplifying the Eq. 17 with the
relation $\left\\{Y_{k}+\text{i}\ Y_{90,k}\right\\}={{\bf{V}}_{k-1}}{P_{k}}$
one gets:

${{\bf{V}}_{k}}k{\bf{\Sigma}}_{k}{{\bf{V}}_{k}}^{T}={{\bf{V}}_{k-1}}\\{(k-1){\bf{\Sigma}}_{k-1}+{P_{k}}{P_{k}}^{T}\\}{{\bf{V}}_{k-1}}^{T}$
(18)

Under finitely large sample-size $k$ the term
$\\{(k-1){\bf{\Sigma}}_{k-1}+{P_{k}}{P_{k}}^{T}\\}$ exhibits a diagonally
dominant characteristics since ${\bf{\Sigma}}_{k-1}$ is diagonally dominant.
Towards obtaining a mapping between both the sides of Eq. 18, Gershgorin’s
theorem ensures the EVD of the term to be of the form
${\bf{\Psi}}_{k}{{{\bf{\Omega}}_{k}}}{\left({{\bf{\Psi}}_{k}}\right)^{T}}$,
where ${{\bf{\Psi}}_{k}}=\left({\bf I}+\delta{\bf V}\right)$ and
${{\bf{\Omega}}_{k}}=\left({{\bf{\Sigma}}_{k-1}}+\delta{\bf\Sigma}\right)$ are
the orthonormal eigenvector and eigenvalue perturbation matrices,
respectively, at $k^{th}$ instant. Upon substitution and simplification the
Eq. 18 can be rewritten as follows:

${{\bf{V}}_{k}}k{\bf{\Sigma}}_{k}{{\bf{V}}_{k}}^{T}={{\bf{V}}_{k-1}}{\bf{\Psi}}_{k}{{{\bf{\Omega}}_{k}}}{\left({{\bf{\Psi}}_{k}}\right)^{T}}{{\bf{V}}_{k-1}}^{T}$
(19)

On close observation of the Eq. 19, final updates for the eigenvector and
eigenvalues can be identified as:

${\bf{V}}_{k}={{\bf{V}}_{k-1}}{\bf{\Psi}}_{k}\quad\&\quad{\bf{\Sigma}}_{k}=\frac{1}{k}{{\bf{\Omega}}_{k}}$
(20)

A drawback for using FOEP technique to estimate the eigenspace at each instant
is associated with permutation ambiguity i.e. the recursive eigenvectors and
eigenvalues are not arranged in same sequence as the theoretical values.
However the shortcoming can be overcome by rearranging the eigenvalues in
decreasing order, and, correspondingly the eigenvectors. The complex whitening
data at a particular time instant can be extracted as follows:

$\left\\{Z_{k}+\text{i}\
Z_{90,k}\right\\}={{\bf{\Sigma}}_{k}}^{-1/2}{{\bf{V}}_{k}}^{T}\left\\{Y_{k}+\text{i}\
Y_{90,k}\right\\}={{\bf{W}}_{k}}^{T}\left\\{Y_{k}+\text{i}\ Y_{90,k}\right\\}$
(21)

Here, $Z_{k}$ is the whitening data and $Z_{90,k}$ is its 90∘ phase shift of
the whitening data, and, ${{\bf{W}}_{k}}^{T}$ is the whitening matrix. This
accomplishes the second objective of the proposed algorithm and the next step
is to update the lagged covariance matrix formed through the whitened data.
For a series of lags ($\tau=1,\ldots,\mathcal{T}$), the augmented lagged
whitening covariance matrix is formed as follows,

$\left[{\begin{array}[]{*{20}{c}}{{\bf{R}}_{k,k-1}}\\\ {{\bf{R}}_{k,k-2}}\\\
\vdots\\\
{{\bf{R}}_{k,k-\mathcal{T}}}\end{array}}\right]=\frac{{k-1}}{k}\left[{\begin{array}[]{*{20}{c}}{{\bf{R}}_{(k-1),(k-1)-1}}\\\
{{\bf{R}}_{(k-1),(k-1)-2}}\\\ \vdots\\\
{{\bf{R}}_{(k-1),(k-1)-\mathcal{T}}}\end{array}}\right]+\frac{1}{k}\left[{\begin{array}[]{*{20}{c}}\left\\{Z_{k}+\text{i}\
Z_{90,k}\right\\}\left\\{Z_{k-1}+\text{i}\ Z_{90,k-1}\right\\}^{T}\\\
\left\\{Z_{k}+\text{i}\ Z_{90,k}\right\\}\left\\{Z_{k-2}+\text{i}\
Z_{90,k-2}\right\\}^{T}\\\ \vdots\\\ \left\\{Z_{k}+\text{i}\
Z_{90,k}\right\\}\left\\{Z_{k-\mathcal{T}}+\text{i}\
Z_{90,k-\mathcal{T}}\right\\}^{T}\end{array}}\right]$ (22)

where, ${\bf{R}}_{k,k-\tau}$ is the lagged whitening covariance matrix with
lag of $\tau$. After obtaining augmented lagged whitening covariance matrix
the problem of modal identification reduces to finding a best suitable
candidate for unitarily diagonalizing $\left[{\bf{R}}_{k,k-1}\
{\bf{R}}_{k,k-2}\ \ldots\ {\bf{R}}_{k,k-\mathcal{T}}\right]$ at each instance
of time. Denoting
$\mathbf{D}^{\tau}=\mathbf{U}_{k}^{T}\left[{\bf{R}}_{k,k-\tau}\right]\mathbf{U}_{k}$,
the problem of finding the minimum of the performance index $J$ can be
mathematically written as,

$J({\bf U},\tau)\mathop{\rm=}\sum\limits_{\mathcal{T}}{\sum\limits_{1\leq
i\neq j\leq m}{\left|{D_{ij}^{\tau}}\right|^{2}}}$ (23)

Then, for minimum $J$ over fixed $h$ iterations the unitary matrix
$\mathbf{U}$ is said to be an approximate joint diagonalizing the augmented
lagged whitening covariance matrix. Having obtained the whitening and unitary
matrix, the mixing matrix is obtained using the following relation

$\left[{\bf{A}}+\text{i}\ {\bf{A}}_{90}\right]=\mathbf{U}^{\dagger}\mathbf{W}$
(24)

The mixing matrix is then operated on the obtained sensor response $Y_{k}$ to
estimate the complex modal response as,

$\left\\{\left(X\right)_{k}+\text{i}\
X_{90,k}(t)\right\\}=\left[{\bf{A}}_{0,k}+\text{i}\
{\bf{A}}_{90,k}\right]^{\dagger}Y_{k}$ (25)

Thereafter, both the complex identified mode shape matrix and modal response
are normalized to obtain the real valued modal mixing matrix and responses
which are independent to each other.

$\begin{array}[]{ll}\left({\bf{A}}_{ij}\right)_{k}&=\frac{1}{k}\|X_{k}+\text{i}\
X_{90,k}\|_{i}\
sgn\left(\left({\bf{A}}_{ij}\right)_{0,k}\right)\left\|\left({\bf{A}}_{ij}\right)_{0,k}+\text{i}\
\left({\bf{A}}_{ij}\right)_{90,k}\right\|\\\
X_{k}&=\frac{1}{k}\|X_{k}+\text{i}\ X_{90,k}\|\cdot X_{k}\end{array}$ (26)

Figure 1: A framework for the proposed technique Algorithm 1 Algorithm for the
proposed method

${\bf{X}}=[X_{1},X_{2},...,X_{L}]^{T}\quad\forall\quad X\in{\mathbb{R}^{n}}$
$\triangleright$ Response vector of system

Appending 90∘ phase shifted response $\triangleright$ Complex response, Eq.
(14)

Initialization: complex eigenspace, whiten data and its lagged covariance
matrix

for $k=1$ : $N_{iteration}$ do$\triangleright$ Iteration for sample length

define memory depth parameter $\lambda$

Define ${P_{k}}={{\bf{V}}_{k-1}^{T}\left\\{Y_{k}+\text{i}\
Y_{90,k}(t)\right\\}}$

Estimation of the perturbation matrices ${\bf{\Psi}}_{k}$ and
${\bf{\Omega}}_{k}$

Update of eigenvector and eigenvalue matrices $\triangleright$ Eigenspace
update, Eq. (20)

Estimation of whitened data $\left\\{Z_{k}+\text{i}\ Z_{90,k}(t)\right\\}$ and
whitening matrix ${\bf{W}}_{k}$ $\triangleright$ Eq. 21

Update of lagged whitened covariance matrix ${\bf{R}}_{k,k-\tau}$
$\triangleright$ Eq. (22)

Joint approx diagonalization on ${\bf{R}}_{k,k-\tau}$:
$\mathbf{D}^{\tau}=\mathbf{U}_{k}^{T}\left[{\bf{R}}_{k,k-\tau}\right]\mathbf{U}_{k}$

Estimation of complex mixing matrix using whitening matrix ${\bf{W}}_{k}$ and
unitary matrix ${\bf{U}}_{k}$ $\triangleright$ Eq. (24)

Estimation of complex modal response $\triangleright$ EVD of Covariance, Eq.
(25)

Estimation of real valued modal matrix and response $\triangleright$
Normalization, Eq. (25)

Output: Normalized real valued modes (${\bf A}_{k}$) and Modal response
(${X}_{k}$)

Figure 2: Representation of 3DOF linear lumped mass model

##  4 Case studies for systems with complex mode shapes

The proposed algorithm’s efficacy in identifying the modal response of
structural systems with complex modes is demonstrated through numerical case
studies on linear systems with closely spaced modes and non-proportional
damping. The algorithm’s performance is evaluated on three structures: (i) a
3-DOF linear structural system with closely spaced modes, (ii) a 3-DOF linear
system with non-proportional damping, and (iii) an IASC-ASCE SHM benchmark
structure with complex modes. The numerical simulation results and modal
identification outcomes for these structures are presented and discussed in
detail.

Figure 3: Frequency domain response of the (a) physical (b) modal displacement
of 3DOF linear system with closely spaced modes

###  4.1 3DOF linear structural system with closely spaced modes

A 3-DOF linear vibrating structure modelled as mass, spring and dash-pot
system as presented in Fig. 2 with its governing expression given as in Eq.
(2). The purpose of using WGN is to excite all the modes of a structure in
order to obtain accurate modal parameters and mode shapes. The Monte Carlo
method is a variance reduction technique that can be used to obtain an
ensemble average of the mode shape matrix, source signal, and modal parameters
by generating multiple realizations of the input excitation and system
response. This ensures that the results are more robust and accurate, as the
average of multiple simulations helps to reduce the impact of random
fluctuations and measurement noise. The state-space representation for the
system subjected to a force vector $F(t)$ can be written as,

$\begin{array}[]{l}\dot{U}={\bf{A}}U+{\bf{B}}F\\\ Y={\bf{C}}U\end{array}$ (27)

where, $U$ = vector of states and $Y$ = system response vector governed by the
C matrix. The state matrix, A, and the excitation matrix B are given by,

$\begin{array}[]{l}{\bf{A}}=\left[{\begin{array}[]{*{20}{c}}{{{\left[{\bf{0}}\right]}_{3\times
3}}}&{{{\left[{\bf{I}}\right]}_{3\times 3}}}\\\
{-{{\bf{M}}^{-1}}{\bf{K}}}&{-{{\bf{M}}^{-1}}{\bf{C}}}\end{array}}\right]\\\
{\bf{B}}={\left[{\begin{array}[]{*{20}{c}}0&0&0&diag\left({\bf{M}}^{-1}\right)^{T}\end{array}}\right]^{T}}\end{array}$
(28)

where, ${\bf{M}}$,${\bf{C}}$ and ${\bf{K}}$ are the general system matrix
given as,

${\bf M}=\left[{\begin{array}[]{*{20}c}{1.5}&{0}&{0}\\\ {0}&{2}&{0}\\\
{0}&{0}&{1.3}\end{array}}\right],\;{\bf
K}=\left[{\begin{array}[]{*{20}c}{7}&{-2}&{0}\\\ {-2}&{4}&{-2}\\\
{0}&{-2}&{5}\end{array}}\right],\quad{\bf C}=\alpha{\bf M}+\beta{\bf K}$ (29)

The natural frequencies corresponding to the undertaken system parameters are
$0.1649Hz$, $0.3248Hz$, $0.3655Hz$. The damping ratio for the system is kept
at $\zeta=2.0\%$ critical for each mode. The real valued mode shape can be
obtained from the eigenvalue decomposition of ${\bf M^{-1}K}$ matrix, however,
the actual nature of the mode shapes (real or complex) can be depicted by the
eigen decomposition of the state matrix ${\bf A}$. For the sufficiently close
mode shapes, the mode shape matrix will be complex in nature and for the
undertaken system the real valued and complex mode shapes without
normalization are evaluated as,

$\begin{array}[]{ll}{\bf{{\phi}^{act}}}&=\left[{\begin{array}[]{*{20}c}{-0.2299}&{-0.4061}&{-0.6700}\\\
{-0.6193}&{-0.1528}&{0.3051}\\\
{-0.3437}&{0.7369}&{-0.3288}\end{array}}\right],\\\
{\bf\phi_{c}^{act}}&=\left[{\begin{array}[]{*{20}c}{0.0043+0.2143i}&{0.0052+0.2089i}&{-0.0033-0.3318i}\\\
{0.0116+0.5774i}&{0.0020+0.0786i}&{0.0015+0.1511i}\\\
{0.0064+0.3204i}&{-0.0095-0.3790i}&{-0.0016-0.1628i}\end{array}}\right]\end{array}$
(30)

Figure 4: Actual and Identified real valued mode shapes of the 3DOF linear
system with closely spaced modes

The proposed framework takes the complex response (obtained through appending
the 90∘ phase shifted response) as input and updates the data eigenspace
recursively through FOEP to obtain the whitening matrix. Thereafter, the
whitened data is utilized for the updation of initial lagged covariance
matrix, which are then jointly diagonalized for the estimation of unitary
matrix. Thereafter both recursive whitening and unitary matrices are used to
obtain the mixing and demixing matrices which are then made use of to extract
the system’s source signal. Fig. 3 describes the physical and modal
displacement response in frequency domain. In Fig. 3(b) it can be observed
that the algorithm accurately separates the modal responses of all the modes.
The above observations can further be explored broadly through the MAC values,
which provides a correlation between the actual and identified mode shapes.
The MAC value for the identified modes are obtained as $0.9999$, $0.9910$ and
$0.9835$ for modes 1,2 and 3, respectively, which validates the accuracy of
the identified modes of the system. The plot of actual and identified real
valued modes in Fig. 4 demonstrates the usefulness of the proposed framework
in real-world applications where traditional modal identification methods may
not be effective due to the presence of complex modes.

###  4.2 3DOF linear structural system with non-proportional damping

As shown in the previous section, proportional damping is often used as a
simplified approach to model the effect of damping in linear vibrational
mechanical systems. However, there are cases in which a general viscous
damping is needed to simulate the dynamic of the system with sufficient
accuracy. To validate our technique for the complex modes arising due to the
non-proportional damping, consider the generic differential equation given as
in Eq. (2) with system matrices as follows:

${\bf M}=\left[{\begin{array}[]{*{20}c}{3}&{0}&{0}\\\ {0}&{2}&{0}\\\
{0}&{0}&{1}\end{array}}\right],\;{\bf
C}=\left[{\begin{array}[]{*{20}c}{0.1856}&{0.2290}&{-0.9702}\\\
{0.2290}&{0.0308}&{-0.0297}\\\
{-0.9702}&{-0.0297}&{0.1241}\end{array}}\right],\;{\bf
K}=\left[{\begin{array}[]{*{20}c}{4}&{-2}&{0}\\\ {-2}&{4}&{-2}\\\
{0}&{-2}&{10}\end{array}}\right]$ (31)

The input excitation is modelled as White Gaussian noise (WGN) with sufficient
intensity to excite all the modes properly, and then a variance reduction
technique viz, Monte Carlo is employed to obtain the ensemble average of the
mode shape matrix, source signal and modal parameters. The actual real valued
and complex mode shapes for the undertaken system without normalization is
provided in Eq. (32).

Figure 5: Frequency domain response of the (a) physical (b) modal displacement
of 3DOF linear system with non-proportional damping

$\begin{array}[]{ll}{\bf{{\phi}^{act}}}&=\left[{\begin{array}[]{*{20}c}0.4575&0.3521&0.0090\\\
0.4264&-0.5510&-0.1206\\\ 0.0918&-0.1446&0.9852\end{array}}\right],\\\
{\bf\phi_{c}^{act}}&=\left[{\begin{array}[]{*{20}c}-0.0370-0.0021i&-0.0115+0.2947i&-0.5458+0.0000i\\\
0.0050+0.0364i&-0.0059-0.4345i&-0.5156-0.0126i\\\
-0.0051-0.2988i&-0.0594-0.1167i&-0.1100-0.0521i\end{array}}\right]\end{array}$
(32)

Figure 6: Actual and Identified real valued mode shapes of the 3DOF linear
system with non-proportional damping Figure 7: MAC between actual and
identified modes with increase in samples representing accurate identification
of modes for 3DOF system

Comparison of the estimated modal frequency response in Fig. 5(b) with the
original frequency response in Fig. 5(a) provides a clear depiction of
accurate modal responses separation. This is achieved through the use of
proposed framework on complex responses (obtained through appending the 90∘
phase shifted response). The above observations are substantiated through the
MAC values, which for the identified modes are obtained as $0.9996$, $0.9986$
and $0.9878$ for modes 1,2 and 3, respectively. This validates the accuracy of
the identified modes of the system, and Fig. 6 shows the approximate overlap
between identified modes over true modes. From the results presented in the
study it is evident that the modal identification method was able to estimate
the complex modes. This establishes the applicability of the technique to
estimate complex modes for closely spaced modes and non-proportional damping.

###  4.3 3DOF linear structural system with sudden change in damping

Most real-world phenomena have to accommodate dynamic changes in system
matrices due to various operational and environmental conditions. To overcome
this limitation, the proposed algorithm uses a recursive method to update the
system matrices in real-time, which allows the identification of the changing
modes. This real-time monitoring and updating of the modes provides more
accurate information about the system behavior, which can be useful in many
applications, such as real-time SHM, damage detection, and control systems.
The proposed algorithm provides a robust solution for the identification of
complex modes in dynamic systems, which can handle non-proportional damping,
non-linearity, gyroscopic effects, and closely spaced modes. In order to
validate the applicability of the proposed algorithm, a general 3DOF linear
system with proportional damping is considered with the system matrices
defined as,

${\bf M}=\left[{\begin{array}[]{*{20}c}{3}&{0}&{0}\\\ {0}&{2}&{0}\\\
{0}&{0}&{2}\end{array}}\right],\;{\bf
K}=\left[{\begin{array}[]{*{20}c}{4}&{-2}&{0}\\\ {-2}&{4}&{-2}\\\
{0}&{-2}&{10}\end{array}}\right],\quad{\bf C}=\alpha{\bf M}+\beta{\bf K}$ (33)

Here, the damping coefficient is 2$\%$ and the system is excited with WGN for
a time period of 50s and sampling frequency is taken as 50Hz. This system is
considered to be in its prime state and has real valued modes. At 25s the
damping matrix suddenly becomes non-proportional and is given as,

${\bf C}=\left[{\begin{array}[]{*{20}c}{0.1856}&{0.2290}&{-0.9702}\\\
{0.2290}&{0.0308}&{-0.0297}\\\
{-0.9702}&{-0.0297}&{0.1241}\end{array}}\right]$ (34)

The perturbations in the damping matrix convert the real-valued modes into
complex-valued modes – presenting a more formidable challenge to conventional
blind source separation (BSS) algorithms – as the precision of the mode
identifications can not be substantiated in practical conditions, nor can the
alterations be monitored. Nonetheless, the proposed algorithm encompasses the
requisite tools to effectively monitor such alterations and accurately
determine both real-valued modes (i.e. proportional damping within the
interval of 0 to 25 seconds) and complex-valued modes (i.e. non-proportional
damping within the interval of 25 to 50 seconds). To further verify the
applicability of the proposed algorithm in practical scenarios, a benchmark
system is described in the subsequent section.

Figure 8: ASCE-SHM benchmark model adopted from [26]

###  4.4 IASC-ASCE SHM benchmark structure with complex modes

Having established the viability of the proposed algorithm through numerical
simulation of simple systems, this study progresses to the examination of its
application to the IASC-ASCE Structural Health Monitoring (SHM) Benchmark
structure. The latter constitutes a quarter-scale model of a four-story, two-
bay by two-bay steel frame structure located at the University of British
Columbia’s Earthquake Engineering Research Laboratory. It boasts a planar
dimension of 2.5 m x 2.5 m and a height of 3.6 m, with members fabricated from
hot rolled grade 300 W steel, possessing a nominal yield stress of 300 MPa.
The structural sectional properties can be consulted in [26]. The parameters
utilized to generate the response of the benchmark structure are: a damping
ratio of 0.01, a sampling time step size of 0.002s, a total time duration of
50s, a noise level of 1, and a force coefficient of 150.

Figure 9: Frequency domain response of the (a) physical (b) modal displacement
of modified ASCE SHM benchmark structure Figure 10: MAC between actual and
identified modes representing accurate identification of modes for ASCE-SHM
benchmark structure Figure 11: Original and identified real valued mode shapes
of the ASCE SHM benchmark system with non-proportional damping for different
samples

In the IASC-ASCE SHM Benchmark structure, the columns are aligned in the x
direction, which exhibits superior resistance to bending, while the floor
beams are aligned to withstand strong vertical bending. Each bay on each floor
features a single floor slab, with the first, second, third, and fourth levels
comprising four slabs weighing 800 kg, 600 kg, 600 kg, and 400 kg,
respectively, for the symmetrical building configuration. In contrast, the
fourth level in the unsymmetrical building configuration includes three 400 kg
slabs and one 550 kg slab. To generate simulated response data, a finite
element model is devised for various conditions. In the absence of damage, a
12-degree-of-freedom (DOF) shear-building model that restricts all motion
(with the exception of translations in the x and y directions and a single
rotation $\theta$ about the center column) by modeling the floor beams and
floor slabs as rigid bodies is created, resulting in three DOFs per floor. The
columns and floor beams are modeled as Euler-Bernoulli beams in the finite
element model. A diagram of the analytical model is shown in Fig. 8 (case 1)
with no damage in the structure (damage pattern 0) is taken for the validation
of the proposed algorithm. For the undamaged case, the mass (${\bf{M}}$) and
stiffness (${\bf{K}}$) matrices of the structure are given by the following
matrices:

$\begin{array}[]{l}{\bf{M}}=\left[{\begin{array}[]{*{20}{c}}{3452.4}&0&0&0&0&0&0&0&0&0&0&0\\\
0&{3452.4}&0&0&0&0&0&0&0&0&0&0\\\ 0&0&{3819.4}&0&0&0&0&0&0&0&0&0\\\
0&0&0&{2652.4}&0&0&0&0&0&0&0&0\\\ 0&0&0&0&{2986.1}&0&0&0&0&0&0&0\\\
0&0&0&0&0&{2652.4}&0&0&0&0&0&0\\\ 0&0&0&0&0&0&{2652.4}&0&0&0&0&0\\\
0&0&0&0&0&0&0&{2986.1}&0&0&0&0\\\ 0&0&0&0&0&0&0&0&{1809.9}&0&0&0\\\
0&0&0&0&0&0&0&0&0&{1809.9}&0&0\\\ 0&0&0&0&0&0&0&0&0&0&{1809.9}&0\\\
0&0&0&0&0&0&0&0&0&0&0&{2056.9}\end{array}}\right]\text{kg}\\\ \end{array}$

(35)

$\begin{array}[]{l}{\bf{K}}=\left[{\begin{array}[]{*{20}{c}}{213.20}&0&0&{-106.60}&0&0&0&0&0&0&0&0\\\
0&{135.81}&0&0&{-67.90}&0&0&0&0&0&0&0\\\
0&0&{464.04}&0&0&{-232.02}&0&0&0&0&0&0\\\
{-106.60}&0&0&{213.20}&0&0&{-106.60}&0&0&0&0&0\\\
0&{-67.90}&0&0&{135.81}&0&0&{-67.90}&0&0&0&0\\\
0&0&{-232.02}&0&0&{464.04}&0&0&{-232.02}&0&0&0\\\
0&0&0&{-106.60}&0&0&{213.20}&0&0&{-106.60}&0&0\\\
0&0&0&0&{-67.90}&0&0&{135.81}&0&0&{-67.90}&0\\\
0&0&0&0&0&{-232.02}&0&0&{464.04}&0&0&{-232.02}\\\
0&0&0&0&0&0&{-106.60}&0&0&{106.60}&0&0\\\
0&0&0&0&0&0&0&{-67.90}&0&0&{67.90}&0\\\
0&0&0&0&0&0&0&0&{-232.02}&0&0&{232.02}\end{array}}\right]\text{MN/m}\end{array}$

(36)

The IASC-ASCE SHM Benchmark structure is symmetric, except in cases where
damage results in loss of symmetry, and the applied forces are limited to the
y direction. In such scenarios, the output measurements indicate no response
in the x or $\theta$ direction, and only the y direction measurements are
utilized in the identification process to reduce any numerical instability.
Given that the system in question has sparsely separated modes and mass- and
stiffness-proportional damping, the actual modes are real and can be obtained
as follows:

$\begin{array}[]{ll}{\bf{{\phi}^{act}}}&=\left[{\begin{array}[]{*{20}c}0.2422&-0.6226&0.5029&0.2129\\\
0.4414&-0.4297&-0.5032&-0.5591\\\ 0.5803&0.1956&-0.3494&0.6559\\\
0.6399&0.6241&0.6098&-0.4603\end{array}}\right]\end{array}$ (37)

To introduce complex modes into the benchmark structure for the purposes of
our study, an additional mass with a stiffness and damper is added to the top
story, affecting only the y direction. The new dynamic system has a mass of
$172.6$ kg, a stiffness of $6.8$ MN/m, and a 2% damping ratio. With the
addition of this new dynamic system, the modes of the system in question
become complex in the y direction and can be characterized as follows:

$\begin{array}[]{ll}{\bf\phi_{c}^{act}}&=\left[{\begin{array}[]{*{20}c}0.0980-0.6579i&0.2123-2.0306i&0.1644+3.8407i&-0.2735-3.7886i\\\
-0.0310+1.8428i&0.1897+2.1281i&-0.2647+2.7033i&-0.2480-7.2326i\\\
-0.0123-2.1900i&-0.0904+1.4109i&-0.2072-1.1655i&-0.2244-9.6458i\\\
0.0105+1.5428i&-0.0524-2.5213i&-0.0388-3.8303i&-0.2230-10.6845i\end{array}}\right]\end{array}$

(38)

The proposed algorithm was also applied to a benchmark structure subjected to
stochastic excitation. The comparison between the physical and estimated modal
responses in the frequency domain is presented in Figures 9(a) and 9(b),
respectively, which shows a clear separation of the modal responses,
indicating that the modes are indeed complex. The accuracy of the identified
modes was further validated through the Maximum Absolute Coefficient (MAC)
values, which are plotted against the number of samples in Figure 10. The
comparison between the actual and identified real-valued modes, as shown in
Figure 11, demonstrates a close match between the identified modes and the
true modes, as evidenced by the high MAC values. These results demonstrate the
applicability of the proposed algorithm for accurately estimating complex
modes in a wide range of numerical simulations and real-life structural
models.

##  5 Conclusions

The findings of this research have firmly established the capability of the
proposed algorithm in identifying complex modes in a wide range of dynamic
systems in real-time. The algorithm is built upon the principle of first-order
eigen perturbation (FOEP) and is designed to synchronize with the complex
response, updating the second-order statistics recursively. Detection of
features of interest may not necessarily be related to sudden or gradual
damage, but can also be related to parameters representing their performance
with respect to control, repair or change in other operational conditions. In
many practical vibrating systems, the modes can either inherently be or become
complex due to dynamic modifications of the system. This study therefore
provides compelling evidence for the real-time estimation of complex modes and
the accurate identification of modal properties, even in the presence of
dynamic changes to the system’s characteristics. The proposed algorithm
effectively overcomes the limitations of traditional batch algorithms and
real-time techniques, providing accurate real-time identification of complex
modes even in the presence of closely spaced modes, non-proportional damping,
dynamically changing damping matrices, and the addition of new dynamical
systems. The results of this study have been demonstrated through extensive
numerical simulations and benchmark case studies, making it a highly reliable
and robust approach for real-time evaluation of complex modes.

## Acknowledgement

Satyam Panda acknowledges the financial support received from Prime Minister
Research Fellowship. Budhaditya Hazra gratefully acknowledge the financial
support received from Science and Engineering Research Board (SERB),
Department of Science and Technology (DST), Government of India (Project no.
IMP/2019/000276).

## References

  * [1] WF Xue, J Chen, JQ Li, and XF Liu.  Acoustical feature extraction of rotating machinery with combined wave superposition and blind source separation.  Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 220(9):1423–1431, 2006.
  * [2] B Hazra, A Sadhu, and S Narasimhan.  Fault detection of gearboxes using synchro-squeezing transform.  Journal of Vibration and Control, 23(19):3108–3127, 2017.
  * [3] Shashank Chauhan, Dmitri Tcherniak, and Morten H Hansen.  Dynamic characterization of operational wind turbines using operational modal analysis.  Proceedings of China Wind Power, pages 13–15, 2010.
  * [4] A Sadhu, S Narasimhan, and J Antoni.  A review of output-only structural mode identification literature employing blind source separation methods.  Mechanical Systems and Signal Processing, 94:415–431, 2017.
  * [5] YANJIE ZHU, ANDRE JESUS, IRWANDA LAORY, and YUAN REN.  Thermal effect identification and bridge damage disclosure by using blind source separation method.  Structural Health Monitoring, pages 527–534, 2017.
  * [6] Maria Giovanna Masciotta and Daniele Pellegrini.  Tracking the variation of complex mode shapes for damage quantification and localization in structural systems.  Mechanical Systems and Signal Processing, 169:108731, 2022.
  * [7] M Imregun and DJ Ewins.  Complex modes-origins and limits.  In Proceedings-SPIE The International society for optical engineering, pages 496–496. SPIE INTERNATIONAL SOCIETY FOR OPTICAL, 1995.
  * [8] Ulrich Fuellekrug.  Computation of real normal modes from complex eigenvectors.  Mechanical systems and signal processing, 22(1):57–65, 2008.
  * [9] Jérôme Antoni and S Chauhan.  A study and extension of second-order blind source separation to operational modal analysis.  Journal of Sound and Vibration, 332(4):1079–1106, 2013.
  * [10] Budhadtiya Hazra, Ayan Sadhu, Aaron J Roffel, and Sriram Narasimhan.  Hybrid time-frequency blind source separation towards ambient system identification of structures.  Computer-Aided Civil and Infrastructure Engineering, 27(5):314–332, 2012.
  * [11] Rune Brincker and Carlos Ventura.  Introduction to operational modal analysis.  John Wiley & Sons, 2015.
  * [12] Wei Liu, Wei-cheng Gao, and Yi Sun.  Application of modal identification methods to spatial structure using field measurement data.  Journal of vibration and acoustics, 131(3), 2009.
  * [13] Basuraj Bhowmik, Tapas Tripura, Budhaditya Hazra, and Vikram Pakrashi.  First-order eigen-perturbation techniques for real-time damage detection of vibrating systems: Theory and applications.  Applied Mechanics Reviews, 71(6), 2019.
  * [14] Satyam Panda, Tapas Tripura, and Budhaditya Hazra.  First-order error-adapted eigen perturbation for real-time modal identification of vibrating structures.  J Vib Acoust,, 143(5), 2021.
  * [15] Satyam Panda, Tapas Tripura, and Budhaditya Hazra.  Online damage detection of earthquake-excited structure based on near real-time envelope extraction.  Structural Health Monitoring, 21(2):298–319, 2022.
  * [16] Basuraj Bhowmik, Satyam Panda, Budhaditya Hazra, and Vikram Pakrashi.  Feedback-driven error-corrected single-sensor analytics for real-time condition monitoring.  International Journal of Mechanical Sciences, 214:106898, 2022.
  * [17] SI McNeill and DC Zimmerman.  A framework for blind modal identification using joint approximate diagonalization.  Mechanical Systems and Signal Processing, 22(7):1526–1548, 2008\.
  * [18] Kiran D’Souza, Chulwoo Jung, and Bogdan I Epureanu.  Analyzing mistuned multi-stage turbomachinery rotors with aerodynamic effects.  Journal of Fluids and Structures, 42:388–400, 2013.
  * [19] Rakesh K Goel.  Seismic behaviour of asymmetric buildings with supplemental damping.  Earthquake engineering & structural dynamics, 29(4):461–480, 2000\.
  * [20] Seungjin Choi, Andrzej Cichocki, Hyung-Min Park, and Soo-Young Lee.  Blind source separation and independent component analysis: A review.  Neural Information Processing-Letters and Reviews, 6(1):1–57, 2005\.
  * [21] Xi-Ren Cao and Ruey-wen Liu.  General approach to blind source separation.  IEEE Transactions on signal Processing, 44(3):562–571, 1996.
  * [22] Adel Belouchrani, Karim Abed-Meraim, J-F Cardoso, and Eric Moulines.  A blind source separation technique using second-order statistics.  IEEE Transactions on signal processing, 45(2):434–444, 1997.
  * [23] Carlo Rainieri.  Perspectives of second-order blind identification for operational modal analysis of civil structures.  Shock and Vibration, 2014, 2014.
  * [24] Gaëtan Kerschen, Fabien Poncelet, and J-C Golinval.  Physical interpretation of independent component analysis in structural dynamics.  Mechanical Systems and Signal Processing, 21(4):1561–1575, 2007\.
  * [25] Paul Mucchielli, Basuraj Bhowmik, Budhaditya Hazra, and Vikram Pakrashi.  Higher-order stabilized perturbation for recursive eigen-decomposition estimation.  Journal of Vibration and Acoustics, 142(6), 2020.
  * [26] Erik A Johnson, Heung-Fai Lam, Lambros S Katafygiotis, and James L Beck.  Phase i iasc-asce structural health monitoring benchmark problem using simulated data.  Journal of engineering mechanics, 130(1):3–15, 2004.


# Simulations of hybrid charge-sensing single-electron-transistors and CMOS
circuits

Tetsufumi Tanamoto  Department of Information and Electronic Engineering,
Teikyo University, Toyosatodai, Utsunomiya 320-8551, Japan
<EMAIL_ADDRESS>Keiji Ono  Advanced device laboratory, RIKEN,
Wako-shi, Saitama 351-0198, Japan

###### Abstract

Single-electron transistors (SETs) have been extensively used as charge
sensors in many areas such as quantum computations. In general, the signals of
SETs are smaller than those of complementary metal-oxide semiconductor (CMOS)
devices, and many amplifying circuits are required to enlarge these signals.
Instead of amplifying a single small output, we theoretically consider the
amplification of pairs of SETs, such that one of the SETs is used as a
reference. We simulate the two-stage amplification process of SETs and CMOS
devices using a conventional SPICE (Simulation Program with Integrated Circuit
Emphasis) circuit simulator. Implementing the pairs of SETs into CMOS circuits
makes the integration of SETs more feasible because of direct signal transfer
from the SET to the CMOS circuits.

Single-electron transistors (SETs) have been intensively investigated owing to
their advantages of low-power operation, which is desirable for application to
logic and memory elements Grabert ; Likharev . Since Tucker’s proposal Tucker
, many approaches have been developed to replace the elements in complementary
metal-oxide semiconductor (CMOS) circuits Chen ; Tiwari ; Yano ; Mahapatra ;
Inokawa ; Uchida . SETs were also investigated after they were directly
embedded in CMOS circuits Yano ; Uchida2 ; Inokawa . Currently, conventional
Si transistors are much smaller than SETs. SETs have attracted attention as
charge sensors Field ; Berman ; Ionescu , which have been used for the
readouts of silicon qubits Morello ; Gonzalez-Zalba ; Shaji , or as current
standards Fujiwara ; Nishiguchi ; Fujiwara2 .

The SET consists of a small metallic island or a degenerated semiconductor
island surrounded by a source and drain via tunneling junctions with low
capacitances. The SET yields periodic outputs referred to as Coulomb
oscillations, which vary as a function of the gate voltage, and the
oscillation period corresponds to the change in the number of electrons in the
island. The maximum signal change in an SET is the difference between the peak
and the trough of the Coulomb oscillation. The SET current is sensitive to
electric potential variations on the SET island Field . The standard readout
process of a spin qubit (spin in localized states, such as a single impurity
or quantum dot (QD)) is a spin-to-charge conversion, and the change in charge
of the localized state is detected by the SET current Schoelkopf ; Lu ;
Morello ; Gonzalez-Zalba .

There are many methods for the detection of charges, such as the use of
radiofrequency SETs (rf-SETs) Schoelkopf ; Lu ; Gonzalez-Zalba  and the direct
measurement of the SET by a single GaAs amplifier Petersson ; Visscher .
Inokawa et al showed that a Coulomb oscillation is effectively outputted by
directly connecting the SET to an MOS field-effect transistor. They
experimentally demonstrated multiple-valued logic by SETs operating at 27 K.
The effectiveness of the series coupling of the SET with MOS transistors was
also experimentally investigated by Uchida et al. Uchida2 . Scalable SET
sensing systems require scalable circuits. However, previous SET sensors
Morello ; Gonzalez-Zalba ; Shaji ; Fujiwara ; Nishiguchi ; Fujiwara2  did not
explicitly consider the array of SETs as part of CMOS circuits.

Figure 1:  Concept of our two-stage amplification circuit of the charge-
sensing single-electron-transistor (SET). Instead of amplifying one SET
sensor, we consider amplifying pairs of SETs that have different states with
each other. Depending on the existence of the extra electron in the quantum
dot (QD) outside the SETs, the electric potential of the island of the SETs
changes. Throughout this paper, we model the effects of the electric potential
of the QDs and $V_{\rm SET}$ by the gate voltages $V_{\rm G1}$ and $V_{\rm
G2}$ of the SETs. The switches between the 1st and 2nd stages are typically
realized by wordline nMOS transistors, using which we can select pairs of SETs
among the SET array.

In this study, we theoretically consider a scalable detection circuit based on
the implementation of many SETs in CMOS circuits. Our basic concept of the
amplification of the SET signals consists of two stages, as shown in Fig. 1.
The first amplification stage of the process is conducted by direct connection
of the SET to pMOS transistors. In the second stage, we introduce a reference
and a target SET and amplify the difference between two SETs using standard
amplifier circuits, such as the differential amplifier (DA) circuits, and the
static random-access memory (SRAM) cells.

For smooth connection to the digital circuits, the output of the sensing
circuit should be a digital signal (either a ”0” or a ”1”). For this purpose,
it is better to compare the relative output of the target SET with that of the
reference SET. In our application of the SRAM cell, the relative voltage
difference between the target SET and the reference SET is found to be quickly
latched to ”0” or ”1.” This is in contrast to the measurements in Refs.
Morello ; Gonzalez-Zalba ; Shaji , wherein the results were obtained after the
analysis of a series of time-dependent SET currents. Our CMOS circuits are
assumed to be close to the SETs. In addition, the target and the reference
SETs are chosen by switching on the wordline transistors between the pairs of
SETs and amplifiers. Then, the circuit using the pairs of SETs becomes more
compact than the amplifying circuit of a single SET. Accordingly, our proposal
is suitable for an array of sensing SETs.

In this study, we implement the current characteristics of the orthodox theory
of the SET Amman  into the SPICE (Simulation Program with Integrated Circuit
Emphasis) circuit simulator based on the BSIM4 (Berkeley Short-channel
transistor Model, level=54) by using the standard modeling language of the
Verilog-A. We considered two types of SETs one of which is operated at low
(4.2 K), and the other is operated at high temperatures (-30∘C). At low
temperatures, such as 4.2 K, the threshold voltage of the MOS transistors
becomes higher because of the incomplete ionizations. Many studies regarding
cryo-CMOS Green ; Gaensslen ; Beckers ; Dijk  have been conducted to determine
the model parameters at low temperatures. However, the general compact model
is not available for such low temperature. In addition, the basic CMOS
operations are basically the same as that at room temperature (RT). Thus, we
applied the CMOS parameters that are available in the conventional SPICE
models to the CMOS parts for both types of the SETs. In the following, circuit
calculations are mainly conducted at $T=-30^{\circ}$C, which is in the range
of the conventional models.

We consider CMOSs with gate lengths of $L=$90 nm and $L=$65 nm, whose drain
voltages $V_{\rm D}$ are 1.2 V and 1.0 V, respectively. We also consider the
effects of small variations in the SETs and CMOS transistors. The purpose of
the comparison of pairs of SETs is to detect the changes in the Coulomb
oscillations of the target SET. Thus, given that the voltage difference
between the peak and trough currents of the Coulomb oscillations can be
distinguished, we will be able to detect the changes of the target SET even if
there are variations in the devices.

Table 1:  The four SETs we use here.  | Capacitance(aF) | Resistance($\Omega$) | Temp(K) |  $E_{c}$(meV)
---|---|---|---|---
SET | $C_{\rm up}$ | $C_{\rm dn}$ | $C_{g}$ | $R_{\rm up}$ | $R_{\rm dn}$ | $T$ |
LT-SET$r$ | 1 | 10 | 2 | 100k | 1M | 4.2 | 6.16
LT-SET$t$ | 1.1 | 0.9 | 2.2 | 110k | 0.9M | 4.2 | 6.51
RT-SET$r$ | 0.1 | 0.5 | 0.5 | 100k | 500k | 243 | 72.8
RT-SET$t$ | 0.11 | 0.45 | 0.45 | 110k | 450k | 243 | 79.3
Figure 2:  First-stage amplification circuits, where the LT-SET$t$ (see table
I) and pMOS are directly connected. The output voltage $V_{\rm out}$ is
amplified depending on the gate voltage $V_{\rm Gp}$ of the pMOS. Depending on
the width of the pMOS ($W_{p}=$ 0.5 $\mu$m and $W_{p}=$ 1 $\mu$m), the optimal
points change. The Coulomb oscillation is amplified up to approximately 10 mV
($V_{\rm D}=1.2$ V). Hereafter we use $W_{p}=$ 0.5 $\mu$m pMOS at the 1st
stage.

Herein, we consider four SETs, as listed in Table I. The LT-SET$t$ (target)
and LT-SET$r$ (reference) are operated at 4.2 K, and RT-SET$t$ and RT-SET$r$
can operate at -30∘C. The parameters of the LT-SET$t$ (RT-SET$t$) exhibit 10%
variations compared with those of the LT-SET$t$ (RT-SET$r$). The calculated
current–voltage characteristics ($I_{\rm D}$-$V_{\rm D}$) are listed in the
appendix. The charging energy is given by $E_{c}\equiv e^{2}/[2(C_{\rm
up}+C_{\rm dn}+C_{\rm g})]$. The magnitude between the peak and trough
currents is on the order of pA, and its voltage change estimated by pA$\times$
25.9 k$\Omega=$ 25.9 $\mu$V is very small (25.9 k$\Omega$ is a quantum
resistance Grabert ). On the other hand, the variations in the threshold
voltage $V_{\rm th}$ of conventional CMOSs are generally on the order of
millivolts. Therefore, we need to amplify the SET outputs before connecting
the SETs to conventional CMOS circuits.

Figure 3:  (a) The differential amplifier (DA) is applied at the 2nd stage
amplification of Fig.1. (b) Numerical results of $V_{\rm out1}$-$V_{\rm out2}$
(enhanced Coulomb oscillations) for $L=$90 nm at $V_{\rm G1}=0$. $V_{\rm
Gp1}=V_{\rm Gp2}=V_{\rm Gp}$. The pMOS and nMOS widths of the DA are given by
$W_{\rm pa}=1\mu$m, and $W_{\rm na}=$10 $\mu$m, respectively. The widths of
the wordline nMOS and 1st stage pMOS are given by $W_{\rm n}=5$ $\mu$m and
$W_{\rm p}=$ 0.5 $\mu$m, respectively. $V_{\rm D}=1.2$ V. We can see that the
amplitude of $V_{\rm out2}-V_{\rm out1}$ is approximately 40 mV, which is
larger than the amplitude of the input $V_{\rm in2}-V_{\rm in1}$. (c)
Histogram of the Monte Carlo simulation of the output difference when the
$V_{\rm th}$ of all CMOS transistors varies by 10% over 300 simulations.

We start from the first amplification stage in which the SET is directly
connected to the CMOS transistors, as shown in Fig. 2(a). By a series
connection, the low SET current value increases with CMOS Inokawa . Figure
2(b) shows that the amplification of the SET signal becomes prominent at
$V_{\rm Gp}\gtrsim 0.9$ V. The distortions of the waveforms compared with the
original Coulomb oscillations originate from the nonlinearity of the $I_{\rm
D}$-$V_{\rm D}$ characteristics of the pMOS transistors. We analyze the
amplification mechanism based on the standard long-channel model Taur . The
$I_{\rm D}$-$V_{\rm D}$ of CMOS transistors depends on the triode region
$|V_{\rm DS}|<|V_{\rm GS}-V_{\rm th}|$ and the saturation region $|V_{\rm
DS}|>|V_{\rm GS}-V_{\rm th}|$. At present, amplification is observed in
$|V_{\rm DS}|=|V_{\rm out}-V_{\rm D}|\approx$1.05 V$>|V_{\rm GS}|=|V_{\rm
Gp}-V_{\rm D}|=0.15$ V (saturation region). We did not observe a sufficiently
large amplification in the triode region. The SET current $I_{\rm SET}$ under
large source-drain voltage is approximately described by $I_{\rm SET}\propto
V_{\rm SET}$ (see appendix). More explicitly we assume that $I_{\rm
SET}\approx V_{\rm SET}/R_{\rm D}$. Thus, we can write the $I_{\rm D}$-$V_{\rm
D}$ characteristics in the saturation region; further, the SET given by

$\displaystyle I_{\rm D}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\beta_{p}(V_{\rm Gp}-V_{\rm D}-V_{\rm
thp})^{2}(1+\lambda(V_{\rm D}-V_{\rm out})),$ (1) $\displaystyle V_{\rm out}$
$\displaystyle\approx$ $\displaystyle R_{\rm D}I_{\rm D}.$ (2)

The solution of these equations is given by

$V_{\rm out}=\frac{1+\lambda V_{\rm D}}{1+\lambda I_{D0}R_{\rm D}}I_{D0}R_{\rm
D},$ (3)

where $I_{D0}\equiv\frac{1}{2}\beta_{p}(V_{\rm Gp}-V_{\rm D}-V_{\rm
thp})^{2}$. Because $I_{D0}R_{\rm D}<V_{\rm D}$, the SET output is enhanced by
$\frac{1+\lambda V_{\rm D}}{1+\lambda I_{D0}R_{\rm D}}>1$. Note that $\lambda$
is conspicuous in both the $L=$90 nm and $L=$65 nm pMOS transistors (see
appendix), and we can see that the output voltage oscillates around $\sim$10
mV in Fig. 2(b).

Figure 4: (a) Six transistor static random-access memory (SRAM) cell applied
in the 2nd-stage amplification of Fig. 1. (b) Time dependent voltage behaviors
of the SRAM setup of $L$=90 nm for widths $W_{n}=0.5$ $\mu$m and
$W_{p}=1.2W_{n}$. The width of the $WL$ transistor is 0.4 $\mu$m. $V_{\rm
Gp}=0.55$ V and $V_{\rm G1}=0.0$ V using the LT-SET$t$s. The change in $V_{\rm
out1}$ and $V_{\rm out2}$ is in the range of $V_{\rm G2}=\\{$0, 0.05, 0.1,
0.15, 0.2, 0.25 0.3, 0.35, 0.4, 0.45, 0.5 $\\}$ V. (c) Replotting of (b) as a
function of $V_{\rm G2}$ of (b) and the result for $L=$65 nm.

Figure 5:  (a) Dynamic RAM(DRAM)-like detection circuits are applied at the
2nd stage amplification of Fig. 1. The stray capacitance of 1 pF is added to
the nodes $V_{\rm out1}$, $V_{\rm out2}$, inputs of the SETs, and equalizer
circuits. (b)(c) Time-dependent characteristics of the different SET pairs
(the LT-SET pair in (b) and RT-SET pair in (c)) in Fig. (a) with 10% threshold
voltage variations in the MOS transistors. The pulse sequence is constituted
following the standard DRAM sequence of Ref. DRAM1 ; DRAM2 . $V_{\rm G1}=0$ is
fixed and $V_{\rm G2}$ changes. Depending on whether $V_{\rm G1}>V_{\rm G2}$
or $V_{\rm G1}<V_{\rm G2}$, the outputs $V_{\rm out1}$ and $V_{\rm out2}$
change from $V_{\rm out1}>V_{\rm out2}$ to $V_{\rm out1}<V_{\rm out2}$. The
widths of the nMOS and pMOS are 0.5 $\mu$m and 0.6 $\mu$m. The width of the
equalizer nMOS is 20 $\mu$m. The width of the $WL$ transistor is 0.6 $\mu$m.
($V_{\rm Gp1}=V_{\rm Gp2}=0.87$ V, $L=$90 nm).

The second-stage amplification is conducted using standard amplifier circuits.
Figure 3 shows the results obtained from a basic DA circuit. Note that the
gate voltages $V_{\rm G1}$ and $V_{\rm G2}$ represent the shifted electrical
potential $V_{\rm SET}$ by the additional sensing QD in Fig. 1. Conventional
DAs amplify two inputs with opposite phases. In the case of SETs, two
different phases of Coulomb oscillations are input. We consider that the two
SET signals possessing different current peaks mimic the two input signals of
the conventional DA with different phases. In Figs. 3(b), we show the results
for the output voltage difference of $V_{\rm out2}-V_{\rm out1}$ for a fixed
$V_{\rm G1}=0$. The voltage difference $V_{\rm out2}-V_{\rm out1}$ increases
by approximately 40 mV for $L=90$ nm (the results for $L=65$nm are shown in
the appendix), and can be detected by conventional CMOS amplifier circuits
(the following circuits after $V_{\rm out2}$ and $V_{\rm out1}$ are not
shown). This enhancement of the Coulomb oscillations is the result of the two-
stage amplification of SET signals.

The magnitude of the enhancement of the Coulomb oscillation in Fig 3(b)
changes depending on the threshold voltage variations of the eight MOS
transistors. In Fig.3(c), we provide the distribution of the difference in the
peak and trough of $V_{\rm out2}-V_{\rm out1}$ obtained through over 300 Monte
Carlo simulations of the threshold variations with LT-SET$t$ and LT-SET$r$ at
$V_{\rm G1}=0$ and $V_{\rm G2}=0.2$ V. The number of small amplitudes of
$|V_{\rm out2}-V_{\rm out1}|$ ($<$ 5 meV) is nine out of 300 samples. Small
amplitudes of $V_{\rm out2}-V_{\rm out1}$ can be avoided by applying different
voltages such as $V_{\rm Gpi}$ and $WL_{i}$ on each part of the DA.

We now consider the application of the standard SRAM cell containing six MOS
transistors Seevinck  to detect a pair of SETs, as shown in Fig. 4(a). Figure
4(b) shows the simulation results of the time-dependent SRAM cell outputs
$V_{\rm out1}$ and $V_{\rm out2}$ of the two LT-SET$t$s. We can see that
$V_{\rm out1}<V_{\rm out2}$ for $V_{\rm G2}=0.1$ V, but $V_{\rm out1}>V_{\rm
out2}$ for $0.15$ V$\leq V_{\rm G2}\leq 0.4$V for $L$= 90 nm devices. This
implies that the shift in the electric potential of the sensing QD from
$V_{\rm G2}=0.1$ V to $V_{\rm G2}=0.15$ V changes the electric potential of
the target SET island, resulting in a change in the relative magnitude between
$V_{\rm out1}$ and $V_{\rm out2}$. Figure 4(c) replots $V_{\rm out2}-V_{\rm
out1}$ in Fig. 4(b) as a function of $V_{\rm G2}$ with the results for $L=$ 65
nm at time $10^{-8}$ s. Thus, we can detect the change in the target SET
represented by $V_{\rm G2}$ by measuring the relative outputs of the SRAM
cell. Herein, the initial voltage at the SRAM cell input is set to $V_{\rm
D}/2$, and stray capacitances of 0.2 pF are included at the input nodes of the
SRAM cell (Figures not shown). As the stray capacitance increases, the time to
split increases.

In general, SRAM cells undergo initial threshold voltage $V_{\rm th}$
variations PUF , and we have to consider these variations in the MOS
transistors as well as the two SETs. Here, we extend the SRAM cell circuit to
a dynamic random-access memory (DRAM)-like structure DRAM1 ; DRAM2  in Fig.
5(a), where it is considered that the equalizer circuit mitigates the voltage
difference of the wordlines between the two SETs. Figures 5(b)(c) show two
types of readouts for the two types of SETs. The equalizer voltage is switched
on at 1 ns and stopped after 2 ns. The SAN and SAP voltages are switched on
after the equalizer at $t=6$ ns and $t=8.5$ ns. The wordlines are switched on
at 4 ns. We can see that the change in the gate voltage $V_{\rm G2}$, which
corresponds to the existence of the charge sensor QD, causes the outputs
$V_{\rm out1}$ and $V_{\rm out2}$ to change from $V_{\rm out1}>V_{\rm out2}$
to $V_{\rm out1}<V_{\rm out2}$.

In a realistic situation, it is possible that the temperature changes.
Therefore, we also calculated the temperature dependence of the amplifier
response. It is desirable that the relative magnitude of $V_{\rm out1}$ to
$V_{\rm out2}$ does not change even when the operating temperature changes.
The temperature dependencies of the characteristics of the DA (Fig. 3) and
DRAM (Fig. 5) are robust against temperature changes as long as the
temperature change is sufficiently small on the order of several tens of
degrees (see appendix).

Considering that there are variations in SETs and MOS transistors, we may have
to check and record the basic characteristics of each device at the first
calibration stage of the chip. The $I$–$V$ characteristics of each SET should
be clarified before using the charge-sensing SET. This information is stored
in the extra SRAM or other memories, and becomes the overhead of the system.
An effective method for determining the optimal biases automatically is a
future issue. It is possible that the applied voltages destroy the SET charge
states. However, herein, we neglected the effects of the backaction from the
CMOS circuits. Hence, the assessment of the backaction remains a future issue.

In conclusion, we proposed two-stage amplification circuits for SETs. Based on
serial connections with MOS transistors and a comparison with the reference
SET, we numerically show that the readout of the charge-sensing SET is greatly
enhanced. We also considered the effects of variations in the MOS transistors
and SETs, and show that as long as the variations are small, the two SETs can
be compared effectively.

## DATA AVAILABILITY

The data that supports the findings of this study are available within the
article.

###### Acknowledgements.

We are grateful to T. Mori and H. Fuketa for the fruitful discussions. We are
also grateful to Y. Yamamoto for his technical support in using SmartSpice.
This work was partly supported by MEXT Quantum Leap Flagship Program (MEXT
Q-LEAP) Grant Number JPMXS0118069228, Japan.

##  Appendix A SET characteristics determined by the orthodox theory

Figure 6 shows the current-gate voltage characteristics of the four SETs used
in the main text, calculated based on the orthodox theory Amman . The Coulomb
staircase becomes clearer for a large asymmetry between the two tunnel
junctions; however, the amplitude of the oscillation becomes smaller.
Therefore, we consider the case of a smaller asymmetry. The wide-range view of
the current–voltage ($I_{\rm D}-V_{\rm D}$) characteristics o f the SET shown
in Figs. 6(b)(d)(f)(h) enables the approximation of $V_{\rm D}\propto I_{\rm
D}$ of Eq.(2) in the main text.

Figure 6:  Current-voltage characteristics of SETs calculated based on the
orthodox theory. (a)(c)(e)(g) Temperature dependence of the $I_{\rm
SET}-V_{\rm G}$ characteristics at $V_{\rm D}=0.01$ V. (b)(d) Three-
dimensional (3D) view of $I_{\rm SET}$ for $T$ = 0.1 K. (f)(g) 3D view of
$I_{\rm SET}$ for $T=-30^{\circ}$C. (a)(b)$C_{\rm up}=1$ aF, $C_{\rm dn}=10$
aF, $C_{g}=2$ aF, $R_{\rm up}=100$ k$\Omega$, $R_{\rm dn}=1$ M$\Omega$ (LT-
SET$t$ in Table I). (c)(d)$C_{\rm up}=1.1$ aF, $C_{\rm dn}=9$ aF, $C_{g}=2.2$
aF, $R_{\rm up}=110$ k$\Omega$, $R_{\rm dn}=0.9$ M$\Omega$.(LT-SET$r$ in Table
I). As the asymmetry of the two junctions becomes larger, the Coulomb
staircase becomes clearer. (e)(f)$C_{\rm up}=0.1$ aF, $C_{\rm dn}=0.5$ aF,
$C_{g}=0.5$ aF, $R_{\rm up}=100$ k$\Omega$, $R_{\rm dn}=500$ k$\Omega$ (RT-
SET$t$ in Table I). (g)(h)$C_{\rm up}=0.11$ aF, $C_{\rm dn}=0.45$ aF,
$C_{g}=0.45$ aF, $R_{\rm up}=90$ k$\Omega$, $R_{\rm dn}=450$ k$\Omega$ (RT-
SET$r$ in Table I).  Figure 7:  $I_{\rm D}-V_{\rm D}$ characteristics of the
pMOS transistors at -30∘ C. (a) $L=$90 nm and (b) $L$=65 nm.

##  Appendix B Additional results for the transistors and amplification
circuits

Figure 7 shows the basic $I_{\rm D}$-$V_{\rm D}$ characteristics of the pMOS
transistors used herein (Berkeley short-channel transistor model BSIM-4), and
Fig. 8 shows the characteristics of the DA circuits of Fig. 3 (b) for $L$ = 65
nm. Although amplification is not observed in the triode region, herein, we
show that amplification can be expected even in the triode region. In the
triode region, Eq. (1) is replaced by

$\displaystyle I_{\rm D}\\!$ $\displaystyle=$
$\displaystyle\\!\beta_{p}(V_{\rm Gp}-V_{\rm D}-V_{\rm thp}-[V_{\rm D}-V_{\rm
out}])[V_{\rm D}-V_{\rm out}],\ \ \ \ $ (4)

where $\beta_{p}\equiv\mu_{p}C_{\rm ox}\frac{W}{L}$ and $\lambda(<1)$ are the
channel length modulation coefficients ($L$, $W$, $C_{\rm ox}$, and $\mu_{p}$
are the length, width, gate capacitance and mobility of the pMOS,
respectively). The solution to the equations is given by

$V_{\rm out}\approx V_{\rm D}+2v_{\rm gp}-\frac{1}{R_{\rm D}\beta_{p}},$ (5)

where $v_{\rm gp}=V_{\rm Gp}-V_{\rm D}-V_{\rm thp}$. Eq. (5) shows that the
SET output changes the rate of $R_{\rm D}\partial V_{\rm out}/\partial R_{\rm
D}\approx 1/[\beta_{p}R_{\rm D}]$. When we take $\mu_{p}=0.025$ m2/(Vs),
$C_{\rm ox}=\epsilon_{0}\epsilon/t_{\rm ox}$ with $\epsilon=3.9$ and $t_{\rm
ox}=4$ nm for $L=$90 nm and $W_{p}=0.5~{}\mu$m, we have $\beta_{p}=0.0002397$
AV-2 and $\beta_{p}R_{\rm D}=2.397\times 10^{2}$ V-1 for $R_{\rm
D}=1~{}$M$\Omega$. Then the change rate $R_{\rm D}\partial V_{\rm
out}/\partial R_{\rm D}\approx 4.17$ mV. Whether the triode region or the
saturation region is better for the amplification depends on the current
characteristics of the MOS transistors.

Figure 8:  Numerical results of the two outputs $V_{\rm out1}$ and $V_{\rm
out2}$ of the differential amplifier circuit for $L=$65 nm for a fixed $V_{\rm
G2}=0$. The transistor sizes are the same as those in Fig. 3 in the main text.

##  Appendix C Temperature dependence

We calculated the temperature dependence of the amplifier response. It is
desirable that the relative magnitude of $V_{\rm out1}$ to $V_{\rm out2}$ does
not change as long as the operating temperature changes in the range of
approximately -50${}^{\circ}\sim 0^{\circ}$C. Figures 9 and 10 show the
temperature dependence of the characteristics of the DA (Fig. 3), and DRAM
(Fig. 5). For example, depending on the value of $V_{\rm Gp}$, we can observe
the relative change in the magnitude between $V_{\rm out1}$ and $V_{\rm
out2}$. We also observe that even when there are variations, the two SETs can
be differentiated by changing the gate voltage $V_{\rm G1}$ or $V_{\rm G2}$.
The strength to the change of temperature depends on each operation region.

Figure 9:  Temperature dependence of $V_{\rm out1}-V_{\rm out2}$ in Fig. 3(a)
of $L$=90 nm MOS transistors for different SETs as a function of $V_{\rm G1}$.
$V_{\rm G2}$=0.1 V. The left SET is the RT-SET$t$ and the right SET is the RT-
SET$r$.  Figure 10:  Temperature-dependent characteristics for $V_{\rm
G2}=0.2~{}$V($L$=90 nm). The widths of the nMOS and pMOS in the circuits are
0.5 $\mu$m and 0.6 $\mu$m. The width of the equalizer nMOS is 20 $\mu$m.
$V_{\rm Gp1}=V_{\rm Gp2}=0.87$ V.

Figure 11:  (a) Cross-coupled type readout of two pairs of SETs. (b)(c) Time-
dependent output characteristics of the cross-coupled type assuming 10 %
variations in the threshold voltage of the MOS transistors (the LT-SET pair in
(b) and RT-SET pair in (c)). The two wordline voltages are adjusted to change
the relation between $V_{\rm G1}$ and $V_{\rm G2}$ when $V_{\rm G2}$ changes.
The wordline voltages are switched on at 10 ps with a rising time of 10 ps.
$W_{n}=5~{}\mu$m, $W_{p}=0.1~{}\mu$m and $V_{\rm G1}$=0 V.

##  Appendix D Detection by cross-coupled circuit

Figure 11(a) shows our other proposal based on the cross-coupled type. In this
case, two SETs are series-connected to cross-coupled MOS transistors, where
the changes in the SET signals are directly transferred to the MOS
transistors. Herein, the 10% variations between the SETs and among the MOS
transistors are considered. To observe the change in $V_{\rm G2}$ to the fixed
$V_{\rm G1}$, we have to adjust the wordline voltages. In Figs. 11(b) and
11(c), the wordline voltage of the left SET is reduced to $WL_{1}\approx 0.84$
V, and $WL_{1}\approx 1.08$ V, respectively. In Fig. 11(b)((c)), $V_{\rm
out2}$ becomes higher than $V_{\rm out1}$ between $V_{\rm G2}=-0.55$ V and
$V_{\rm G2}=-0.5$ V ($V_{\rm G2}=-0.3$ V and $V_{\rm G2}=-0.2$ V). When there
are no variations, the left and right hand circuits are exactly the same, and
$V_{\rm out1}$ and $V_{\rm out2}$ are clearly split depending on the gate
voltages of the two SETs. Temperature changes are not observed, similar to the
case with the SRAMs mentioned above (Figures not shown).

## References

  * (1) H. Grabert and M. H. Devoret, ”Single Charge Tunneling: Coulomb Blockade Phenomena In Nanostructures” (Nato Science Series B:)(Springer, New York, 1992).
  * (2) K. K. Likharev, Proc. IEEE, 87, 606 (1999).
  * (3) J. R. Tucker, J. Appl. Phys. 72, 4399 (1992).
  * (4) R. H. Chen, A. N. Korotkov, and K. K. Likharev, Appl. Phys. Lett. 68, 1954 (1996).
  * (5) S. Tiwari et al., Appl. Phys. Lett. 68, 1377 (1996).
  * (6) K. Uchida et al., Jpn. J. Appl. Phys. 38, 4027 (1999).
  * (7) S. Mahapatra and A. M. Ionescu, 4th IEEE Conference on Nanotechnology p.287 (2004).
  * (8) H. Inokawa, A. Fujiwara and Y. Takahashi, IEEE Trans. Electron Devices, 50, 462 (2003).
  * (9) K. Yano et al., 1999 IEEE International Electron Devices Meeting 87, 633 (1999).
  * (10) K. Uchida et al., Digest of ISSCC 2002, 206 (2002).
  * (11) M. Field, C. G. Smith, M. Pepper, D. A. Ritchie, J. E. F. Frost, G. A. C. Jones, and D. G. Hasko, Phys. Rev. Lett. 70 1311 (1993).
  * (12) D. Berman J. Vac. Sci & Tech. B 15, 2844 (1997)
  * (13) A. M. Ionescu, M. J. Declercq, S. Mahapatra, K. Banerjee and J. Gautier, Proceedings 2002 Design Automation Conference pp. 88-93 (2002).
  * (14) A. Morello,et al., Nature 467, 687 (2010).
  * (15) Gonzalez-Zalba et al., Nat Commun 6, 6084 (2015).
  * (16) N. Shaji et al., Nature Phys 4, 540 (2008).
  * (17) A. Fujiwara, et al., Appl. Phys. Lett. 88, 053121 (2006).
  * (18) K. Nishiguchi and A. Fujiwara, 2007 IEEE International Electron Devices Meeting, 791 (2007).
  * (19) S. P. Giblin et al., Metrologia 57, 025013 (2020).
  * (20) R. J. Schoelkopf et al., Science 280, 1238 (1998).
  * (21) W. Lu, Z. Ji, L. Pfeiffer, K.W. West, and A.J. Rimberg, Nature 423, 422 (2003).
  * (22) J. Petersson et al., Phys. Rev. B 53, 13272 (1996).
  * (23) E. H. Visscher et al., Appl. Phys. Lett. 68, 2014 (1996).
  * (24) M. Amman, R. Wilkins, E. Ben-Jacob, P. D. Maker, and R. C. Jaklevic, Phys. Rev. B 43, 1146 (1991).
  * (25) R. R. Green, Rev. Sci. Instrum., 39, 1495 (1968).
  * (26) F. H. Gaensslen, V. L. Rideout, E. J. Walker and J. J. Walker, IEEE Trans. Electron Devices, 24, 218 (1977).
  * (27) A. Beckers, F. Jazaeri, A. Ruffino, C. Bruschini, A. Baschirotto and C. Enz, 2017 47th European Solid-State Device Research Conference (ESSDERC), p. 62, (2017).
  * (28) J. van Dijk et al., 2020 IEEE Custom Integrated Circuits Conference (CICC), p. 1, (2020).
  * (29) Y. Taur and T. H. Ning, ”Fundamentals of Modern VLSI Devices” (Cambridge University Press, 1998).
  * (30) E. Seevinck, F. J. List and J. Lohstroh, IEEE J. Solid-State Circuits, 22, 748 (1987).
  * (31) D. E. Holcomb et al.,IEEE Trans. Comput. 58, 1198 (2009).
  * (32) B. Keeth and R. J. Baker, ”DRAM Circuit Design: A Tutorial” (John Wiley & Sons, 2001).
  * (33) B. Jacob, D. Wang, and S. Ng, ”Memory Systems: Cache, DRAM, Disk” (Morgan Kaufmann; 1st edition, Amserdam 2007).