# Applications of the dynamical generator coordinate method to quadrupole
excitations

N. Hizawa  K. Hagino  K. Yoshida  Department of Physics, Kyoto University,
Kyoto 606-8502, Japan

###### Abstract

We apply the dynamical generator coordinate method (DGCM) with a conjugate
momentum to a nuclear collective excitation. To this end, we first discuss how
to construct a numerically workable scheme of the DGCM for a general one-body
operator. We then apply the DGCM to the quadrupole vibration of 16O using the
Gogny D1S interaction. We show that both the ground state energy and the
excitation energies are lowered as compared to the conventional GCM with the
same number of basis functions. We also compute the sum rule values for the
quadrupole and monopole operators, and show that the DGCM yields more
consistent results than the conventional GCM to the values from the double
commutator. These results imply that the conjugate momentum is an important
and relevant degree of freedom in collective motions.

##  I Introduction

In strongly correlated quantum many-body systems such as nuclei, collective
behaviors often emerge due to cooperation of many degrees of freedom. These
phenomena are called collective motions [1], for which a translational motion,
a rotational motion, a vibrational motion, and nuclear fission are typical
examples. To describe them based on microscopic degrees of freedom is one of
the major goals of many-body theory.

One of the basic microscopic theories describing quantum many-body systems is
the mean-field theory. In this theory, interacting quantum many-body systems
are approximated by those of independent particles in a common mean-field. As
particles move independently with each other, the many-body wave function in
the mean-field approximation cannot include explicitly the quantum correlation
among particles, except for the trivial Pauli correlation. Therefore, the pure
mean-field theory is not designed to describe collective motions.

A way to resolve this problem is to mix several mean-field configurations. An
important question in this method is what states should be mixed in the wave
function. The random phase approximation (RPA) can be regarded as a
configuration mixing method, in which the configurations are restricted to
one-particle one-hole (1p1h) excitations from a given mean field solution [1].
Although the RPA is powerful in describing a small oscillation around the mean
field solution, it is difficult to apply it to large amplitude motions, such
as nuclear fission. In order to describe such large amplitude collective
motions, the generator coordinate method (GCM) [1, 2] has often been used in
nuclear theory. In this method, many-body states which are generated along
some collective coordinate(s) are linearly superposed. Many numerical
calculations have been carried out in recent years [3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18].

A serious theoretical ambiguity inherent in the GCM is that it is difficult to
identify the appropriate collective coordinate for a motion of interest.
Notice that most of recent GCM calculations have been formulated based on a
collective coordinate chosen empirically, but it is not obvious that this
approach can appropriately describe the collective motion. We mention that a
collective coordinate can be unambiguously extracted in principle with the
self-consistent collective coordinate method (SCC) [19, 20] by projecting the
collective motion onto a specific function space. This is a self-contained
theory, and does not require any empirical manipulation. Furthermore, the
collective coordinates and momenta are treated on the same footing. However,
the SCC is formulated in such a way that a mixing of configurations does not
appear explicitly, and its application to the GCM is nontrivial. In this
sense, the GCM is not yet a rigorous framework unlike the density functional
theory (DFT) [21, 22] or the time-dependent DFT (TDDFT) [23, 24, 25], and one
should carefully examine the validity of a chosen generated coordinate. This
problem has been recognized for a long time.

In the case of the translational motion, Peierls and Thouless pointed out that
the naive GCM breaks certain symmetry, and they developed the double
projection method [26]. The double projection has also been formulated for the
rotational motion. Numerical calculations based on this method for the
rotational motion are becoming possible with the development of computer
technology [27, 28, 29, 30, 31, 30, 32], even though it is still
computationally rather expensive.

Since the symmetries are fundamental to this approach, the double projection
has not been extended to more general collective motions. In this connection,
Goeke and Reinhard improved the GCM by introducing the conjugate momentum to
the GCM. This method is referred to as the dynamical GCM (DGCM) [33]. However,
a practical application of this method has not been performed, except for the
one with the Gaussian overlap approximation [34]. In our previous work, we
proposed the implementation of the DGCM based on the constrained variational
method [35]. Based on this idea, we applied the DGCM to the particle number,
which is numerically the easiest case. This method, which could be regarded as
the double projection on the number of particles, is also easy to implement
due to the fact that the method of particle number projection is well known.
The purpose of this paper is to apply the DGCM actually to more general
collective motions and to investigate the effect of the conjugate momentum on
the collective motions. For this purpose, we consider a quadrupole
oscillation, which is often described with the conventional GCM.

The paper is organized as follows. In Sec. II, we give a brief review of the
GCM and the DGCM. We also discuss how to perform the DGCM numerically based on
the mean field theory for general collective motions. In Sec. III, we carry
out the DGCM calculation for 16O with a quadrupole operator and compare the
results with those obtained with the conventional GCM. We then summarize the
paper and discuss future perspectives in Sec. IV.

##  II Formulation

###  II.1 GCM

In this section, we give a brief review of the GCM. The GCM is a class of
configuration mixing, which is performed along some generator coordinate $q$.
Therefore, the GCM starts from creating many-body wave functions
$\\{|q\rangle\\}_{q}$ specified by the generator coordinate $q$. For a
description of collective motions, $q$ is referred to as a collective
coordinate. Of course, $q$ can be a complex number, and its dimension does not
need to be specified. However, for simplicity of the following discussion, we
consider in this paper a one-dimensional real number for $q$. Then, the
configuration mixing is expressed in the integral form as,

$|\Psi\rangle_{\mathrm{GCM}}=\int dq\,f(q)|q\rangle,$ (1)

where $f(q)$ is called a weight function. The subspace spanned by such a
linear combination of $|q\rangle$ is called a collective subspace. In the GCM,
a Hamiltonian of a given system $\hat{H}$ is diagonalized by restricting to
the collective subspace. Therefore, the fundamental equation of the GCM is
obtained by projecting the time-independent Schrödinger equation onto the
collective subspace,

$\langle
q^{\prime}|\left\\{(\hat{H}-E)|\Psi\rangle_{\mathrm{GCM}}\right\\}=0,\quad\mathrm{for}\;\;{}^{\forall}|q^{\prime}\rangle.$
(2)

We then obtain a generalized eigenvalue problem,

$\int dq\left(\langle q^{\prime}|\hat{H}|q\rangle-E\langle
q^{\prime}|q\rangle\right)f(q)=0,$ (3)

where $\mathcal{H}(q^{\prime},q):=\langle q^{\prime}|\hat{H}|q\rangle$ and
$\mathcal{I}(q^{\prime},q):=\langle q^{\prime}|q\rangle$ are referred to as a
Hamiltonian kernel and a norm kernel, respectively. This equation is called
the Hill-Wheeler equation [2].

The most serious problem in the GCM is how to prepare $|q\rangle$. First, we
have to decide what kind of function space to prepare the states $|q\rangle$.
In many cases, Slater determinant wave functions [3, 36] are used, but one may
also superpose more complicated wave functions. Next, we have to decide what
method we use to draw the path $|q\rangle$ on these function spaces. Above all
things, a collective motion is not mathematically well defined in general, and
thus the criteria of justification for a chosen $|q\rangle$ is unclear. The
dynamical GCM, which we review in the next subsection, is theoretically more
robust than the GCM and has certain good properties even when the collective
operator $\hat{Q}$ is chosen empirically.

###  II.2 Dynamical GCM

The dynamical GCM introduced by Goeke and Reinhard can be regarded as the GCM
that simultaneously considers the conjugate momentum $p$ in addition to a
collective coordinate $q$ [33]. Therefore, as in the case of the GCM, the
configuration mixing is performed as

$|\Psi\rangle_{\mathrm{DGCM}}=\iint dqdp\,f(q,p)|q,p\rangle,$ (4)

where this two-parameter path $|q,p\rangle$ is called a dynamical path. The
Hill-Wheeler equation is also obtained in the same way,

$\iint dqdp\left(\langle q^{\prime},p^{\prime}|\hat{H}|q,p\rangle-E\langle
q^{\prime},p^{\prime}|q,p\rangle\right)f(q,p)=0.$ (5)

Notice that Eqs. (4) and (5) alone do not guarantee that $p$ is the conjugate
momentum of $q$, and one additional equation is required. That is, the
conjugate momentum of $p$ has to be determined to satisfy the conjugation
condition,

$\langle
q,p|\overleftarrow{\partial}_{q}\overrightarrow{\partial}_{p}-\overleftarrow{\partial}_{p}\overrightarrow{\partial}_{q}|q,p\rangle=i,$
(6)

where $\overleftarrow{\partial}_{q}$ means the partial derivative of $q$
acting on the left side and similar for $\overleftarrow{\partial}_{p}$. This
condition can be seen as a kind of the canonical commutation relation.

Of course, this conjugate condition has to be supplemented by some boundary
condition to determine the dynamical path. For example, we can construct the
dynamical path $|q,p\rangle$ that satisfies the condition [35],

$|q,p\rangle=e^{ip\hat{Q}}|q\rangle,$ (7)

with the constrained condition, $\langle q|\hat{Q}|q\rangle=q$, where
$\hat{Q}$ is some Hermitian operator related to a collective motion. Since
there are innumerable states satisfying the condition $\langle
q|\hat{Q}|q\rangle=q$, the state $|q\rangle$ must be uniquely determined in
some way. When calculating $|q\rangle$ in the mean-field theory, one often
uses a constrained variational method. Therefore, we shall use this method in
this paper as well. Considering the DGCM constructed in this way, one can show
that under certain conditions the collective subspace defined with the DGCM
has good properties [37]. As a result, an unexpected quantum entanglement is
less likely to occur in the states obtained as solutions of the Hill-Wheeler
equation. See Ref. [37] for details.

It should be noted that the condition (7) assumes the existence of the
operator $\hat{Q}$. Identifying $\hat{Q}$ for a given collective motion is a
difficult problem in general and this remains so even with the DGCM.
Therefore, we choose $\hat{Q}$ empirically in the rest of this paper.

###  II.3 Dynamical path

In performing numerical calculations with the DGCM, we need to consider
specifically how to construct $e^{ip\hat{Q}}|q\rangle$. We first make an
assumption that $\hat{Q}$ is a one-body operator. We then expand $\hat{Q}$
using some representation of creation and annihilation operators
$\\{\hat{c}^{\dagger}_{i},\hat{c}_{i}\\}_{i=1}^{\infty}$ as

$\hat{Q}=\sum_{i,j=1}^{\infty}Q^{c}_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}.$ (8)

Here, $Q^{c}_{ij}$ must be a Hermitian matrix because of Hermicity of
$\hat{Q}$. Using the Campbell-Baker-Hausdorff formula, we obtain

$\displaystyle e^{ip\hat{Q}}\hat{c}^{\dagger}_{k}e^{-ip\hat{Q}}$
$\displaystyle=\sum_{l=1}^{\infty}\left(e^{ipQ^{c}}\right)_{lk}\hat{c}^{\dagger}_{l},$
(9) $\displaystyle e^{ip\hat{Q}}\hat{c}_{k}e^{-ip\hat{Q}}$
$\displaystyle=\sum_{l=1}^{\infty}\left(e^{-ipQ^{c}}\right)_{lk}\hat{c}_{l}.$
(10)

Let us consider the case where $|q\rangle$ is given in the Hartree-Fock
approximation as,

$|q\rangle=\prod_{i=1}^{A}\hat{c}^{\dagger}_{i}|0\rangle,$ (11)

where $A$ is the number of nucleons and $|0\rangle$ is the vacuum state. Note
that we choose the representation of
$\\{\hat{c}^{\dagger}_{i},\hat{c}_{i}\\}_{i=1}^{\infty}$ so that $|q\rangle$
can be written as Eq. (11). In this case, the dynamical path is calculated as

$e^{ip\hat{Q}}|q\rangle=\prod_{i=1}^{A}\left(\sum_{j_{i}}\left(e^{ipQ^{c}}\right)_{j_{i}i}\hat{c}^{\dagger}_{j_{i}}\right)|0\rangle.$
(12)

Thus, the only thing which we have to do for preparing the dynamical path is
to evaluate a unitary matrix,

$X^{\dagger}_{lk}:=\left(e^{ipQ^{c}}\right)_{lk}.$ (13)

The largest advantage of using the Hartree-Fock approximation is that Eq. (12)
remains a Slater determinant since $X$ is a unitary matrix. This implies that
the conventional numerical technique of the GCM, i.e., the generalized Wick
theorem, can be used also for the DGCM calculation. The same statement can be
applied also to the Hartree-Fock-Bogoliubov approximation.

An extention to a more general case of $|q\rangle$ is trivial given that any
state belonging to the Fock space is written by

$|\Psi\rangle=f(\\{\hat{c}^{\dagger}_{i},\hat{c}_{i}\\}_{i=1}^{\infty})|0\rangle,$
(14)

where $f$ is any polynomial function. We can calculate the dynamical path as
in the case of the Hartree-Fock approximation,

$e^{ip\hat{Q}}|\Psi\rangle=f\left(\left\\{\sum_{j}\hat{c}^{\dagger}_{j}X^{\dagger}_{ji},\,\,\sum_{j}X_{ij}\hat{c}_{j}\right\\}_{i=1}^{\infty}\right)|0\rangle.$
(15)

In general, it is important to note that even when the original $|\Psi\rangle$
has a simple form it may not be so after the operator $e^{ip\hat{Q}}$ is
multiplied to $|\Psi\rangle$. Therefore, even if we can evaluate
$e^{ip\hat{Q}}|\Psi\rangle$, it is another issue whether the DGCM can be
executed numerically. In this regard, the scheme (11) and (12) based on the
mean-field theory has the great advantage to implement the DGCM numerically.

In numerical cases, the full Fock space cannot be considered in practice. Thus
we must limit the Fock space to some model space. We define the model space
using a set of $M\in\mathbb{N}$ creation and annihilation operators,
$\\{\hat{b}^{\dagger}_{i},\hat{b}_{i}\\}_{i=1}^{M}$ and expand the operator
$\hat{c}_{i}$ as

$\hat{c}_{i}=\sum_{j=1}^{M}C_{ij}\hat{b}_{j},$ (16)

where $C$ is a $M\times M$ unitary matrix. Then, we can determine the
dynamical path using the operators

$\hat{d}_{i}=\sum_{j,k=1}^{M}X_{ij}C_{jk}\hat{b}_{k}=:\sum_{j,k=1}^{M}D_{ik}\hat{b}_{k}.$
(17)

When we take the limit $M\to\infty$, the matrix $D$ becomes a unitary matrix.
However, if $M$ is finite, $D$ is not a unitary matrix in general. In such
cases, we can not apply the generalized Wick theorem to the dynamical path,
and the advantage of using the mean field theory is lost.

There are several approximate ways to resolve this problem. From a numerical
point of view, the best way is to project the dynamical path onto the Slater
determinant type wave functions. Such projected state is obtained by
maximizing the overlap,

$\max_{|\psi\rangle}|\langle q,p|\psi\rangle|,$ (18)

where $|\psi\rangle$ spans the entire Slater determinant type wave functions
inside the Fock space constructed by
$\\{\hat{b}^{\dagger}_{i},\hat{b}_{i}\\}_{i=1}^{M}$. However, this method is
computationally expensive and it is difficult to evaluate $X$ in many cases.
It would be easier to introduce a cutoff to Eq. (8) for the one body operator
$\hat{Q}$ as

$\displaystyle\hat{Q}^{\mathrm{cut}}$
$\displaystyle:=\sum_{i,j=1}^{M}Q^{c}_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}=\sum_{i,j,k,l=1}^{M}C^{\dagger}_{ki}Q^{c}_{ij}C_{jl}\hat{b}^{\dagger}_{k}\hat{b}_{l}.$
(19)
$\displaystyle=:\sum_{i,j=1}^{M}Q^{b}_{ij}\hat{b}^{\dagger}_{i}\hat{b}_{j}$
(20)

Then, the dynamical path is calculated as

$X_{ij}^{\mathrm{cut}\,\dagger}=\left(e^{ipQ^{c}}\right)^{\mathrm{cut}}_{ij},$
(21)

where $(A)^{\mathrm{cut}}$ means that the matrix $A$ is evaluated as a
$M\times M$ matrix. In this way, $X_{ij}^{\mathrm{cut}}$ is a $M\times M$
unitary matrix, and we can treat the dynamical path as a Slater determinant
type wave function.

Using this simple cut-off method, we can calculate the dynamical path in the
following steps. At first, we diagonalize matrix $Q^{b}$ in Eq. (20),

$Q^{b}=U\Lambda^{Q}U^{\dagger},$ (22)

where $\Lambda^{Q}$ is a diagonal matrix with the eigenvalues of $Q^{b}$. Note
that all the matrices in this equation depend only on the basis representation
$\\{\hat{b}^{\dagger}_{i},\hat{b}_{i}\\}_{i=1}^{M}$ and the constraint
operator $\hat{Q}$. Next, we calculate the matrix $X^{\rm{cut}}$ as

$\left(X^{\mathrm{cut}}\right)^{\dagger}=CUe^{ip\Lambda^{Q}}U^{\dagger}C^{\dagger}.$
(23)

Finally, we get

$D=CUe^{-ip\Lambda^{Q}}U^{\dagger}.$ (24)

In this way, the dynamical path is entirely determined by the matrix $D$ and
one can only store $D$ in numerical calculations.

##  III Numerical applications of the DGCM to the quadrupole excitation

###  III.1 Numerical details

In this paper, we apply the DGCM to the dynamical path specified by the
quadrupole operator, which is often used in conventional GCM calculations.
This operator is defined as

$\hat{Q}_{20}=\int
d^{3}r\,r^{2}Y_{20}(\theta,\phi)\sum_{\sigma,\tau}\hat{\psi}^{\dagger}(\bm{r},\sigma,\tau)\hat{\psi}(\bm{r},\sigma,\tau),$
(25)

where $\hat{\psi}(\bm{r},\sigma,\tau)$ is a nuclear field operator with spin
$\sigma$ and isospin $\tau$. Since this operator has the dimension of the
square of the length, we introduce a dimensionless variable in a conventional
way,

$\hat{\beta}=\frac{\sqrt{20\pi}}{3r_{0}^{2}A^{5/3}}\hat{Q}_{20},$ (26)

where $r_{0}=1.2\,\mathrm{fm}$ and $A$ is the number of nucleons. We use
$\hat{\beta}$ for the collective operator $\hat{Q}$ in Eq. (8) in Sec. II. In
the following, we use the notation $(\beta,p_{\beta})$ for $(q,p)$.

We prepare the many-body state $|\beta\rangle$ using the constrained Hartree-
Fock method, i.e.

$\delta\left\\{\langle\hat{H}\rangle-\lambda(\langle\hat{\beta}\rangle-\beta)\right\\}=0,$
(27)

where the expectation value is defined as

${\langle\;\cdot\;\rangle}:=\frac{\langle\Psi|\;\cdot\;|\Psi\rangle}{\langle\Psi|\Psi\rangle},$
(28)

with the Slater determinant type wave function $|\Psi\rangle$. Here, $\lambda$
is a Lagrange multiplier. Solving this equation, we obtain the normalized DGCM
(GCM) basis function $|\beta\rangle$ satisfying
$\langle\beta|\hat{\beta}|\beta\rangle=\beta$. For the Hamiltonian $\hat{H}$,
we use the Gogny interaction with the D1S parameter set [38, 39]111Note that
the sign of the spin-orbit force parameter is wrong in the original D1S paper
[39]. The Coulomb interaction is treated without any approximation, and the
center-of-mass correction is evaluated by subtracting the energy of the
center-of-mass motion from $\hat{H}$ beforehand. Therefore, the Hamiltonian
$\hat{H}$ used in this paper reads,

$\hat{H}=\sum_{i=1}^{A}\frac{\hat{\bm{p}}_{i}^{2}}{2m}+\frac{1}{2}\sum_{i\neq
j}\hat{V}(i,j)-\frac{1}{2mA}\left(\sum_{i=1}^{A}\hat{\bm{p}}_{i}\right)^{2},$
(29)

where $\hat{\bm{p}}_{i}$ is the momentum operator of the $i$-th particle, and
$\hat{V}(i,j)$ means a sum of the Gogny and the Coulomb interactions between
the $i$-th and the $j$-th particles. Since both the Gogny and the Coulomb
interactions are finite range interactions, the basis expansion method is
effective. To this end, we use the 3-dimensional spherical harmonic oscillator
basis,

$\displaystyle\phi^{\mathrm{HO}}_{n_{x},n_{y},n_{z}}(\bm{r})=$ $\displaystyle
H_{n_{x}}\left(\frac{x}{b}\right)H_{n_{y}}\left(\frac{y}{b}\right)H_{n_{z}}\left(\frac{z}{b}\right)$
$\displaystyle\times\exp\left(-\frac{r^{2}}{2b^{2}}\right)\times\mathrm{normalization},$
(30)

where $b$ is the harmonic oscillator length and $H_{n}$ is a Hermitian
polynomial with the quantum number $n\in\mathbb{N}\cup\\{0\\}$. The basis
functions are truncated at $N=n_{x}+n_{y}+n_{z}\leq N_{0}$, which leads to
$(N_{0}+1)(N_{0}+2)(N_{0}+3)/6$ bases. The spin wave function is introduced in
the form of a direct product with the harmonic oscillator basis, while the
isospin mixing is not considered. Taking into account the spin degree of
freedom, the total number of bases is $(N_{0}+1)(N_{0}+2)(N_{0}+3)/3$ for
protons and neutrons, respectively. Note that we use the Hartree-Fock
approximation in this paper, ignoring the pairing effect.

The dynamical path
$|\beta,p_{\beta}\rangle=e^{ip_{\beta}\hat{\beta}}|\beta\rangle$ is prepared
according to Sec. II. Notice that $|\beta,0\rangle$ so constructed is equal to
$|\beta\rangle$. In this case, the cut-off dimension $M$ in Sec. II C is the
same as $(N_{0}+1)(N_{0}+2)(N_{0}+3)/3$ for neutrons and protons,
respectively. We evaluate the Hamiltonian kernel in the DGCM (GCM) using the
generalized Wick theorem [40]. Notice that the Gogny D1S interaction is the
effective one and has the density dependent term proportional to the 1/3 power
of the nucleon number density. Since the generalized Wick theorem cannot be
treated rigorously in this class of effective interaction, we use an ansatz of
substituting the local transition density in the density-dependent part [41].

Finally, we mention the issue of overcompleteness. In general,
$\\{|\beta\rangle\\}_{\beta}$ or
$\\{|\beta,p_{\beta}\rangle\\}_{\beta,p_{\beta}}$ are not linearly independent
to each other. This linear dependence implies that the Hill-Wheeler equation
has solutions with indefinite energy. This is referred to as the
overcompleteness problem. We can resolve it by artificially removing the zero
eigenvalue modes of the norm kernel. However, in numerical calculations, the
problem of overcompleteness becomes more serious. This is because small
eigenvalues of the norm kernel may cause numerical problems due to numerical
errors, even when the eigenvalues are not exactly zero. Therefore, we must
remove the small eigenvalue modes by hand to achieve stable numerical
calculations. In this paper, we remove the those modes with eigenvalues less
than $10^{-5}$.

###  III.2 Results

We now apply the GCM and the DGCM to the quadrupole excitations of 16O. To
this end, we include the harmononic oscillator basis up to $N_{0}=10$ with the
oscillator length of $b=1.457433\;\mathrm{fm}$, which is the optimized value
in the mean-field approximation with HFBTHO [42]. We first generate the
coordinate $\beta$ in increments of 0.02 in the range of $-$0.24 to 0.24. The
range of $\beta$ is chosen so that a level crossing does not occur, as we do
not include the pairing correlation. We then generate the conjugate momentum
$p_{\beta}/2\pi$ from $-$0.6 to 0.6 in increments of 0.05 for each $\beta$.

Figure 1 shows the energy surface defined as
$E(\beta,p_{\beta}):=\langle\beta,p_{\beta}|\hat{H}|\beta,p_{\beta}\rangle$.

Figure 1: The energy surface of the dynamical path $E(\beta,p_{\beta})$ for
the quadrupole motion of 16O, in units of MeV. The Gogny D1S interaction is
used to draw this figure. The minimum of the energy surface is
$E(0,0)=-129.569\,\mathrm{MeV}$, which is the same as the Hartree-Fock result
without the constraint condition. Because of the time reversal symmetry of
$|\beta,0\rangle$, the energy surface is symmetric about the $p_{\beta}-$axis.

As we can see, we determine the range of $p_{\beta}$ from the condition that
$E(\pm 0.24,0)$ and $E(0,\pm 0.6\times 2\pi)$ are similar to each other. Since
$|\beta\rangle=|\beta,0\rangle$ is an eigenstate of the time reversal operator
in our method, $|\beta,p_{\beta}\rangle$ and $|\beta,-p_{\beta}\rangle$ are
paired for the time reversal operation. This is leads to the inverted symmetry
of the energy surface with respect to the vertical direction in Fig. 1.

Figure 2 shows the root-mean-square (rms) radius
$r_{\mathrm{ms}}(\beta,p_{\beta})$ for each basis, which is defined as

Figure 2: The root-mean-square (rms) radius for the dynamical path
$r_{\mathrm{rms}}(\beta,p_{\beta})$, in units of fm. The rms radius increases
as the deformation $\beta$ increases, while it is flat in the direction of
$p_{\beta}$.

$r_{\mathrm{rms}}(\beta,p_{\beta})=\sqrt{\langle\beta,p_{\beta}|\int
d^{3}r\,r^{2}\hat{\rho}(\bm{r})|\beta,p_{\beta}\rangle/A},$ (31)

where $\hat{\rho}(\bm{r})$ is the particle number density operator of
nucleons. Notice that $r_{\mathrm{rms}}(\beta,p_{\beta})$ is almost
independent of $p_{\beta}$. This means that the operator
$e^{ip_{\beta}\hat{\beta}}$ does not significantly change the rms radius. This
result is not surprising, given that $\hat{\beta}$ and $\int
d^{3}r\,r^{2}\hat{\rho}(\bm{r})$ are commutative in the full space. However,
rigorously speaking, they are not commutative, since we now restrict the model
space to a finite space. Fig. 2 implies that our model space is large enough
and this noncommutativity is practically negligible. Notice that, as in the
rms radius, the operator $e^{ip_{\beta}\hat{\beta}}$ in
$|\beta,p_{\beta}\rangle$ keeps the expectation value of any local deformation
operator,

$\hat{Q}=\sum_{\sigma,\tau,\sigma^{\prime},\tau^{\prime}}\int
d^{3}r\,Q(\bm{r},\sigma,\tau,\sigma^{\prime},\tau^{\prime})\hat{\psi}^{\dagger}(\bm{r},\sigma,\tau)\hat{\psi}(\bm{r},\sigma^{\prime},\tau^{\prime}),$
(32)

invariant in the full space. This implies that dealing with the dynamical path
in the DGCM allows us to introduce the effect of internal excitations for each
$\langle\hat{Q}\rangle$. That is, $|\beta,p_{\beta}\rangle$ all have the same
expectation value of $\hat{Q}$ for each $\beta$, while the energy expectation
value is different. This effect cannot be captured by the conventional GCM
with $|\beta\rangle$.

In performing the DGCM, it is computationally expensive to include all of the
625 points in the $(\beta,p_{\beta})$ plane. We thus select 25 points,
$S^{\mathrm{DGCM}}_{25}=\\{(\beta,p_{\beta})\,|\,\beta\in\\{-0.2,-0.1,0,0.1,0.2\\},\quad
p_{\beta}/2\pi\in\\{-0.6,-0.3,0,0.3,0.6\\}\\}$, to carry out the DGCM
calculations. On the other hand, for the GCM calculations, we use 25 points of
$S^{\mathrm{GCM}}_{25}=\\{(\beta,0)\,|\,\beta\in\\{-0.24,-0.22,\cdots,0.24\\}\\}$,
to compare the DGCM to the GCM under the same conditions as much as possible.

Table 1: The GCM and the DGCM energy for the quadrupole excitations of 16O with the point sets $S^{\mathrm{GCM}}_{25}$ and $S^{\mathrm{DGCM}}_{25}$, respectively. The cut-off for the norm kernel is taken as $10^{-5}$.  state | GCM (MeV) | DGCM (MeV)
---|---|---
1 |  $-$129.682 |  $-$129.765
2 |  $-$107.993 |  $-$108.140
3 |  $-$92.260 |  $-$104.475
4 |  $-$77.911 |  $-$87.019
5 |  $-$64.097 |  $-$83.059

TABLE 1 shows the results of the GCM and the DGCM energies with these point
sets, obtained with the cut-off for the norm kernel of $10^{-5}$. We can see
that the energies are lower for the DGCM for all the states shown in the
table, even though we use the same number of basis functions. This tendency
does not change even when the cutoff is decreased to $10^{-10}$.

The better performance of the DGCM can be understood by investigating the
eigenvalue distributions of the norm kernel. In the DGCM, the smallest
eigenvalue is O($10^{-8}$). On the other hand, in the GCM, there are 10
eigenvalues which are less than 10-13. This implies that the DGCM contains
many more relevant states that are orthogonal to each other.

Figures 3 and 4 show the square of the GCM collective wave function
$|g(\beta)|^{2}$ and that of the DGCM $|g(\beta,p_{\beta})|^{2}$. The left and
the right panels in each figure show the collective wave functions for the
ground and the first excited states, respectively. In the GCM case, the
collective wave function is defined as

$g(\beta)=\int
d\beta^{\prime}\,\mathcal{I}^{1/2}(\beta,\beta^{\prime})f(\beta^{\prime}),$
(33)

where $\mathcal{I}^{1/2}$ is the square root of the norm kernel. In this case,
the collective wave functions for different GCM eigenstates are orthogonal to
each other. The corresponding collective wave function for the DGCM,
$g(\beta,p_{\beta})$, is defined in a similar way.

Figure 3: The square of the collective wave functions for the ground (the left
panel) and the first excited (the right panel) states. The black solid and the
red dashed lines represent the collective wave functions for the GCM and the
DGCM, respectively. For the DGCM wave function, $p_{\beta}$ is set to be 0.
Figure 4: The square of the collective wave functions for the ground (the left
panel) and the first excited (the right panel) states of the DGCM in the two-
dimensional $(\beta,p_{\beta})$ plane.

We can see that the collective wave functions are harmonic oscillator-like in
both the GCM and the DGCM cases. As is evidenced in the collective wave
functions of the DGCM, which are spread over both $\beta$ and $p_{\beta}$, the
conjugate momentum is an important degree of freedom.

Finally, we discuss the results of the sum rule. The sum rule in the GCM is
already known and can be formulated as in the usual case [43]. See Eqs. (44)
and (45) in Appendix A for the explicit forms of the sum rules for the
monopole operator, $\hat{Q}_{0}=\sum_{i=1}^{A}\hat{\bm{r}}_{i}^{2}$ and the
quadrupole operator
$\hat{Q}_{20}=\sum_{i=1}^{A}(2\hat{z}_{i}^{2}-\hat{x}_{i}^{2}-\hat{y}_{i}^{2})$,
respectively. Note that the validity of the usual sum rule depends on the
structure of the subspace (see Appendix A). This fact is known numerically
[44], and it can be seen also with the DGCM and the GCM. TABLE 2 shows the sum
rule results for the quadrupole and the monopole operators based on the GCM
and the DGCM.

Table 2:  Results of the sum rule calculations for the quadrupole and the monopole operators. The values in this table show the ratios of the left-hand side to the right-hand side of Eqs. (44) and (45).  state | quadrupole | monopole
---|---|---
$i$ | GCM | DGCM | GCM | DGCM
1 | 1.0005 | 1.0047 | 0.3855 | 0.9176
2 | 0.6287 | 1.0162 | 0.6533 | 0.9926
3 | 0.7309 | 1.0512 | 0.6184 |  $-$0.0066
4 | 0.6370 | 0.9414 | 0.5818 | 0.6449
5 |  $-$0.5878 | 0.6352 |  $-$0.8390 |  $-$0.0981

The values in this table are the ratios of the left-hand side to the right-
hand side of Eqs. (44) and (45), which should ideally be 1 if the sum rules
are perfectly satisfied. The label $i$ denotes the DGCM (or the GCM) states
corresponding to those in Table 1. For the quadrupole operator, both the GCM
and the DGCM are close to unity for the ground state ($i=1$). In the DGCM,
this is the case also for the excited states up to the third excited state
($i=4$). In contrast, the sum rule is satisfied much less satisfactory in the
GCM. It can thus be concluded that the DGCM is more likely to include relevant
states than the GCM for the quadrupole oscillation. This tendency does not
change even when we use a smaller cut-off for the eigenvalues of the norm
kernel. For example, if we take the cut-off of $10^{-10}$, the sum rule for
the excited states remain unsatisfactory in the GCM, whereas that in the DGCM
is good up to the 7-th excited state. This is related to the fact that under
certain ansatzs the collective subspace in the DGCM can be represented as a
tensor product of the collective and the non-collective degrees of freedom
[37].

Similarly, the results for the monopole operator indicate that the DGCM is
more superior than the GCM. Therefore, the DGCM can prepare different states
more efficiently than the GCM. However, we find that smaller cutoffs do not
improve the results much, both for the GCM and the DGCM. This fact indicates
that the collective subspace in this calculation does not contain enough
degrees of freedom relevant to the monopole operator. This is a natural
consequence of the fact that the collective subspace is prepared in this
calculations primarily for the quadrupole motion, without caring about the
monopole motion.

##  IV Summary and future perspectives

In this paper, we have discussed how to perform actual numerical calculations
for the DGCM, which is an extension of the GCM by including the conjugate
momentum. We have provided a concrete scheme to construct the dynamical path
approximately for general one-body operators. Using this scheme, we have
applied the DGCM to the quadrupole oscillation of ${}^{16}\mathrm{O}$ and have
compared the results with those of the conventional GCM. We have found that
the DGCM leads to more consistent results than the GCM both for the energies
and the sum rules. This suggests that, at least for the quadrupole
oscillations, it is essential to take into account the conjugate momentum.

In this paper, we have neglected the pairing correlation and used the Hartree-
Fock approximation. An obvious future extention of our work is to use the
Hartree-Fock-Bogoliubov (HFB) approximation, with which the pairing effect can
be fully considered. Such extension will be important to discuss large
amplitude collective motions, especially a nuclear fission. Since the HFB
breaks the U(1) symmetry with respect to the particle number, we should
perform a particle number projection to restore the symmetry. The methods in
this paper can be applied immediately to this case. It will be interesting to
see how the conjugate momentum alters the results of the conventional
calculation of the nuclear fission with the WKB approximation.

Theoretical aspects of the DGCM themselves are also interesting topics to
explore. Many conventional GCM calculations are based on microscopic degrees
of freedom, but have a strong phenomenological factor. The main reason for
this is that one does not have good information on a collective subspace. It
is expected that the DGCM is relatively easy to grasp the structure of the
space. This is one of the largest advantages of the DGCM, which allows us to
discuss new aspects for the collective motion, such as entanglement entropy.

## Acknowledgments

The authors thank Y. Hashimoto for numerical advice and N. Shimizu for useful
discussions. This work was supported by JSPS KAKENHI (Grant Nos. 21J22348,
JP19K03824, JP19K03861, and JP19K03872).

##  Appendix A Sum Rules in a truncated space

###  A.1 A general consideration

Let us consider a Hilbert space $\mathcal{H}_{\mathrm{tot}}$ and introduce
linear operators $\hat{A}$ and $\hat{B}$ which map from
$\mathcal{H}_{\mathrm{tot}}$ to itself, $\mathcal{H}_{\mathrm{tot}}$, that is,
$\mathcal{H}_{\mathrm{tot}}\to\mathcal{H}_{\mathrm{tot}}$. We assume that
$\hat{A}$ and $\hat{B}$ are Hermitian operators. We then consider an
orthonormal subset
$\\{|i\rangle\\}_{i\in\mathcal{X}}\subset\mathcal{H}_{\mathrm{tot}}$ that
satisfies

$\langle i|\hat{A}|j\rangle=a_{i}\langle i|j\rangle=a_{i}\delta_{ij},\quad
i,j\in\mathcal{X}.$ (34)

Because $\hat{A}$ is a Hermite operator, $a_{i}$ is a real number. Note that
this equation does not necessarily mean that $|i\rangle$ is an eigenstate of
$\hat{A}$ in the full space. That is, because of a truncation, the subspace
$\mathcal{H}_{\mathrm{sub}}$ spanned by $\\{|i\rangle\\}_{i\in\mathcal{X}}$ is
smaller than $\mathcal{H}_{\mathrm{tot}}$,
$\mathcal{H}_{\mathrm{sub}}\subseteq\mathcal{H}_{\mathrm{tot}}$. Equation (34)
is simply interpreted as a diagonalization of $\hat{A}$ within a restricted
space. The projection operator onto $\mathcal{H}_{\mathrm{sub}}$ can then be
constructed as

$\hat{\pi}_{\mathrm{sub}}=\sum_{i\in\mathcal{X}}|i\rangle\langle i|.$ (35)

This is an identity operator of $\mathcal{H}_{\mathrm{sub}}$.

We consider the following variable,

$\displaystyle S_{AB}(i):=\sum_{j\in\mathcal{X}}(a_{j}-a_{i})|\langle
j|\hat{B}|i\rangle|^{2}$ (36)

When $\hat{A}$ is the Hamiltonian $\hat{H}$, this is nothing but the energy
weighted sum rule. If we assume $[\hat{B},\hat{\pi}_{\mathrm{sub}}]=0$, we
find

$\displaystyle S_{AB}(i)$ $\displaystyle=\sum_{j\in\mathcal{X}}(\langle
j|a_{j}\hat{B}|i\rangle-\langle j|\hat{B}a_{i}|i\rangle)\langle
i|\hat{B}|j\rangle$ (37) $\displaystyle=\langle
i|\hat{B}[\hat{A},\hat{B}]|i\rangle.$ (38)

Since $S_{AB}(i)$ is real, this can also be written as

$\displaystyle S_{AB}(i)$ $\displaystyle=\frac{1}{2}(S_{AB}(i)+S_{AB}(i)^{*})$
(39) $\displaystyle=\frac{1}{2}\langle
i|[\hat{B},[\hat{A},\hat{B}]]|i\rangle.$ (40)

Note that Eq. (40) is trivially satisfied under the assumption
$[\hat{A},\hat{\pi}_{\mathrm{sub}}]=0$. Of course, in this case,
$[\hat{B},\hat{\pi}_{\mathrm{sub}}]$ does not have to be zero.

###  A.2 Explicit form for the quadrupole and the monopole fields

We here consider the case in which the operator $\hat{A}$ is a Hamiltonian
$\hat{H}$ and the operator $\hat{B}$ is a local one-body operatyor
$\sum_{i\in\mathcal{Y}}f(\hat{\bm{r}}_{i})$ with some
$\mathcal{Y}\subseteq\\{1,2,\cdots,A\\}$. Since the Hamiltonian $\hat{H}$
usually contains only up to the second order of the single particle momentum
operator $\hat{\bm{p}}_{\alpha}$, we can evaluate Eq. (40) explicitly. Notice
that the canonical commutation relation leads to a useful formula,

$[\hat{\bm{p}}_{\alpha},f(\hat{\bm{r}}_{i})]=-i\hbar\nabla
f(\hat{\bm{r}}_{i})\delta_{i,\alpha}.$ (41)

We then obtain

$\displaystyle[f(\hat{\bm{r}}_{i}),[\hat{\bm{p}}_{\alpha}\cdot\hat{\bm{p}}_{\beta},f(\hat{\bm{r}}_{k})]]$
$\displaystyle=\hbar^{2}\nabla f(\hat{\bm{r}}_{i})\cdot\nabla
f(\hat{\bm{r}}_{k})\delta_{i,\alpha}\delta_{k,\beta}$
$\displaystyle+\hbar^{2}\nabla f(\hat{\bm{r}}_{k})\cdot\nabla
f(\hat{\bm{r}}_{i})\delta_{i,\beta}\delta_{k,\alpha}.$ (42)

Using these equations, one finally obtains,

$\displaystyle[\hat{B},[\hat{H},\hat{B}]]$
$\displaystyle=\frac{\hbar^{2}}{m}\left(1-\frac{1}{A}\right)\sum_{i\in\mathcal{Y}}(\nabla
f(\hat{\bm{r}}_{i}))^{2}$
$\displaystyle-\frac{\hbar^{2}}{mA}\sum_{\begin{subarray}{c}i,k\in\mathcal{Y}\\\
i\neq k\end{subarray}}\nabla f(\hat{\bm{r}}_{k})\cdot\nabla
f(\hat{\bm{r}}_{i}).$ (43)

The first term in this equation is due to the contributions of the kinetic
energy term and the one-body part of the center-of-mas correction, while the
two-body part of the center-of-mass correction yields the second term.

In this paper, we consider the monopole operator
$\hat{Q}_{0}=\sum_{i=1}^{A}\hat{\bm{r}}_{i}^{2}$ and the quadrupole operator
$\hat{Q}_{20}=\sum_{i=1}^{A}(2\hat{z}_{i}^{2}-\hat{x}_{i}^{2}-\hat{y}_{i}^{2})$
for concrete examples of $\hat{B}$. Note that their normalization factors need
not be taken into account. Therefore, we obtain the sum rules

$\displaystyle
S_{H,Q_{0}}(i)=\frac{2\hbar^{2}}{m}\left(1-\frac{1}{A}\right)\langle
i|\hat{Q}_{0}|i\rangle-\frac{2\hbar^{2}}{mA}\langle
i|\sum_{\begin{subarray}{c}i,k=1\\\ i\neq
k\end{subarray}}^{A}\hat{\bm{r}}_{i}\cdot\hat{\bm{r}}_{k}|i\rangle,$ (44)
$\displaystyle
S_{H,Q_{20}}(i)=\frac{2\hbar^{2}}{m}\left(1-\frac{1}{A}\right)\langle
i|\hat{Q}_{a}|i\rangle-\frac{2\hbar^{2}}{mA}\langle i|\hat{Q}_{b}|i\rangle,$
(45)

where the operators $\hat{Q}_{a}$ and $\hat{Q}_{b}$ are defined as

$\hat{Q}_{a}=\sum_{i=1}^{A}(4\hat{z}_{i}^{2}+\hat{x}_{i}^{2}+\hat{y}_{i}^{2}),$
(46)

and

$\hat{Q}_{b}=\sum_{\begin{subarray}{c}i,k=1\\\ i\neq
k\end{subarray}}^{A}(4\hat{z}_{i}\hat{z}_{k}+\hat{x}_{i}\hat{x}_{k}+\hat{y}_{i}\hat{y}_{k}),$
(47)

respectively.

## References

  * Ring and Schuck [1980] P. Ring and P. Schuck,  _The nuclear many-body problem_ (Springer-Verlag, New York, 1980).
  * Hill and Wheeler [1953] D. L. Hill and J. A. Wheeler, Phys. Rev. 89, 1102 (1953).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Nikšić  _et al._ [2011] T. Nikšić, D. Vretenar, and P. Ring, Prog. Part. Nucl. Phys. 66, 519 (2011).
  * Egido [2016] J. L. Egido, Phys. Scr. 91, 073003 (2016).
  * Robledo  _et al._ [2018] L. M. Robledo, T. R. Rodríguez, and R. R. Rodríguez-Guzmán, J. Phys. G: Nucl. Part. 46, 013001 (2018).
  * Shinohara  _et al._ [2006] S. Shinohara, H. Ohta, T. Nakatsukasa, and K. Yabana, Phys. Rev. C 74, 054315 (2006).
  * Bender and Heenen [2008] M. Bender and P.-H. Heenen, Phys. Rev. C 78, 024309 (2008).
  * Yao  _et al._ [2010] J. M. Yao, J. Meng, P. Ring, and D. Vretenar, Phys. Rev. C 81, 044311 (2010).
  * Rodríguez and Egido [2011a] T. R. Rodríguez and J. L. Egido, Phys. Rev. C 84, 051307 (2011a).
  * Rodríguez and Egido [2011] T. R. Rodríguez and J. L. Egido, Phys. Lett. B 705, 255 (2011).
  * Rodríguez and Egido [2011b] T. R. Rodríguez and J. L. Egido, Phys. Rev. C 84, 051307 (2011b).
  * Yao  _et al._ [2013] J. M. Yao, M. Bender, and P.-H. Heenen, Phys. Rev. C 87, 034322 (2013).
  * Fukuoka  _et al._ [2013] Y. Fukuoka, S. Shinohara, Y. Funaki, T. Nakatsukasa, and K. Yabana, Phys. Rev. C 88, 014321 (2013).
  * Bally  _et al._ [2014] B. Bally, B. Avez, M. Bender, and P.-H. Heenen, Phys. Rev. Lett. 113, 162501 (2014).
  * Yao  _et al._ [2014] J. M. Yao, K. Hagino, Z. P. Li, J. Meng, and P. Ring, Phys. Rev. C 89, 054306 (2014).
  * Yao  _et al._ [2015] J. M. Yao, E. F. Zhou, and Z. P. Li, Phys. Rev. C 92, 041304 (2015).
  * Egido and Jungclaus [2020] J. L. Egido and A. Jungclaus, Phys. Rev. Lett. 125, 192504 (2020).
  * Marumori  _et al._ [1980] T. Marumori, T. Maskawa, F. Sakata, and A. Kuriyama, Prog. Theor. Phys. 64, 1294 (1980).
  * Matsuo [1986] M. Matsuo, Prog. Theor. Phys. 76, 372 (1986).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Runge and Gross [1984] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Nakatsukasa  _et al._ [2016] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Rev. Mod. Phys. 88, 045004 (2016).
  * Casida and Huix-Rotllant [2012] M. Casida and M. Huix-Rotllant, Ann. Rev. Phys. Chem. 63, 287 (2012), pMID: 22242728.
  * Peierls and Thouless [1962] P. Peierls and D. Thouless, Nucl. Phys. 38, 154 (1962).
  * Borrajo  _et al._ [2015] M. Borrajo, T. R. Rodríguez, and J. Luis Egido, Phys. Lett. B 746, 341 (2015).
  * Egido  _et al._ [2016] J. L. Egido, M. Borrajo, and T. R. Rodríguez, Phys. Rev. Lett. 116, 052502 (2016).
  * Rodríguez [2016] T. R. Rodríguez, Eur. Phys. J. A 52, 190 (2016).
  * Shimada  _et al._ [2016] M. Shimada, S. Tagami, and Y. R. Shimizu, Phys. Rev. C 93, 044317 (2016).
  * Shimada  _et al._ [2015] M. Shimada, S. Tagami, and Y. R. Shimizu, Prog. Theor. Exp. Phys 2015, 063D02 (2015).
  * Ushitani  _et al._ [2019] M. Ushitani, S. Tagami, and Y. R. Shimizu, Phys. Rev. C 99, 064328 (2019).
  * Goeke and Reinhard [1980] K. Goeke and P.-G. Reinhard, Ann. Phys. (N. Y.) 124, 249 (1980).
  * Klüpfel  _et al._ [2008] P. Klüpfel, J. Erler, P. G. Reinhard, and J. A. Maruhn, Eur. Phys. J. A 37, 343 (2008).
  * Hizawa  _et al._ [2021] N. Hizawa, K. Hagino, and K. Yoshida, Phys. Rev. C 103, 034313 (2021).
  * Kanada-En’yo  _et al._ [2012] Y. Kanada-En’yo, M. Kimura, and A. Ono, Prog. Theor. Exp. Phys 2012 (2012).
  * [37] N. Hizawa, to be published.
  * Dechargé and Gogny [1980] J. Dechargé and D. Gogny, Phys. Rev. C 21, 1568 (1980).
  * Berger  _et al._ [1991] J. Berger, M. Girod, and D. Gogny, Comput. Phys. Commun. 63, 365 (1991).
  * Balian and Brezin [1969] R. Balian and E. Brezin, Il Nuovo Cimento B (1965-1970) 64, 37 (1969).
  * Bonche  _et al._ [1990] P. Bonche, J. Dobaczewski, H. Flocard, P.-H. Heenen, and J. Meyer, Nucl. Phys. A 510, 466 (1990).
  * Perez  _et al._ [2017] R. N. Perez, N. Schunck, R.-D. Lasseri, C. Zhang, and J. Sarich, Comput. Phys. Commun. 220, 363 (2017).
  * Flocard and Vautherin [1976] H. Flocard and D. Vautherin, Nucl. Phys. A 264, 197 (1976).
  * Johnson  _et al._ [2020] C. W. Johnson, K. A. Luu, and Y. Lu, J. Phys. G: Nucl. Part. 47, 105107 (2020).


# Real-time detection of uncalibrated sensors using Neural Networks

Luis J. Muñoz-Molina  Ignacio Cazorla-Piñar  Juan P. Dominguez-Morales
Fernando Perez-Peña  Altran Innovation Center for Advance Manufacturing
Robotics and Technology of Computers Lab., Universidad de Sevilla, Seville,
41012, Spain.  School of Engineering, Universidad de Cadiz, Spain.

###### Abstract

Nowadays, sensors play a major role in several contexts like science, industry
and daily life which benefit of their use. However, the retrieved information
must be reliable. Anomalies in the behavior of sensors can give rise to
critical consequences such as ruining a scientific project or jeopardizing the
quality of the production in industrial production lines. One of the more
subtle kind of anomalies are uncalibrations. An uncalibration is said to take
place when the sensor is not adjusted or standardized by calibration according
to a ground truth value. In this work, an online machine-learning based
uncalibration detector for temperature, humidity and pressure sensors was
developed. This solution integrates an Artificial Neural Network as main
component which learns from the behavior of the sensors under calibrated
conditions. Then, after trained and deployed, it detects uncalibrations once
they take place. The obtained results show that the proposed solution is able
to detect uncalibrations for deviation values of 0.25º, 1% RH and 1.5 Pa,
respectively. This solution can be adapted to different contexts by means of
transfer learning, whose application allows for the addition of new sensors,
the deployment into new environments and the retraining of the model with
minimum amounts of data.

###### keywords:

Neural Networks , sensors , uncalibrations , sensor anomalies , transfer
learning

††journal: Expert Systems with Applications

##  1 Introduction

Nowadays, sensors play a major role in the connected world [1, 2, 3, 4].
Science, industry or even day-to-day devices integrate sensors to collect the
information coming from the surrounding environment. Ensuring the reliability
and consistency of the information collected is essential in order to
guarantee the proper use of the information. The reliability of the acquired
information depends on the capabilities of the sensor, but also on its
calibration status. An uncalibration is said to take place when the sensor is
not adjusted or standardized by calibration according to a ground truth value.
In general, anomalies can appear in several shapes [5], however the
uncalibration event usually appear in the form of long-term drifts with
different kind of responses such as linear, exponential, logarithmic or,
simply, an irregular drift. Furthermore, there are cases where the magnitude
to control is not directly controlled, that is, the set point is unknown.
Thus, in the latter case, the detection of potential uncalibrations can be
really challenging. As it happens with some other kind of anomalies in sensor
behaviour [6], the detection of uncalibrations is a key issue in order to
ensure fundamental aspects such as reliability, safety or cost-effectiveness
in crucial assets in different contexts, for instance, pharma, energy,
aeronautics, etc. All these properties, make the proposed solution a great
candidate to be deployed in a wide variety of context of both scientific and
industrial and, more precisely, in the context of industry 4.0.

The uncalibration detection problem can be tackled from different perspectives
or paradigms. The available information is a key factor in order to determine
what perspective should be used. On the one hand, the approaches based on the
use of classical techniques can be considered. This is the case of similarity-
based modeling and multivariate analysis [7] or the time series analysis [8].
A survey of some other classical approaches for anomaly derection can be seen
in [9]. It is important to note that the techniques mentioned in the previous
works were used for the detection of specific kinds of anomalies, not
uncalibrations.

On the other hand, due to the huge development of Machine Learning (ML) and
Artificial Intelligence (AI) experienced during the last years, the use of
techniques and approaches based on them has increased [10, 9]. An example of
this fact is the case of [11], where anomalies, defined as an unusually
different behaviour at a specific time, were detected by means of the use of
Hierarchical Temporal Memory. Another interesting example can be seen in [12],
where a solution based on Convolutional Neural Networks (CNN) was developed in
order to detect sensor failures. An example of the application of Neural
networks for anomaly detection in power plants can be found in [13].
Finally,an overview of Neural Network applications to fault diagnosis and
detection in engineering-related systems can be seen in [14]. As in the
classical approaches, it is important to note that all those developments were
devoted to specific anomalies detection, not uncalibration detection in the
sensor system.

In this work, an online uncalibration detector for temperature, humidity and
pressure sensors under unknown conditions (set points unknown) has been
developed. The major novelty of the paper is the use of an Artificial Neural
Network (ANN) as an expert on the sensors behaviour, so that the estimated
behaviour by the ANN can be compared to the one measured and thus, potential
uncalibrations can be detected. Different statistical techniques, such as
goodness of fit and confidence intervals, are used in order to check that the
estimation of the neural network can be considered as equal or different from
the obtained measures.

The developed architecture detects the ongoing uncalibrations and it exhibits
a great flexibility and scalability. Indeed, the obtained results suggest that
the solution can be trained with a generic set of sensors and then be
specialized for more specific contexts with a much lower amount of data. Thus,
for instance, the solution could be trained with a standard set of sensors and
then be specialized for a more specific set devoted to a more concrete task.
Also, including new sensors into the solution is very easy, and fairly low
amount of data is required. Finally, retraining the solution requires data
coming from about one week of data (at a sampling rate of one sample per
minute and per sensor).

The proposed detector has been tested in a real context in the facilities of a
company of the pharmaceutical sector. Thus, it is important to note that the
presented work has been developed under the production conditions of the
aforementioned company.

The rest of the paper is structured as follows: section 2 introduces the
materials used in this work, including the dataset (2.2). The implemented
methods are described in section 3, together with the Neural Network (NN)
model. Then, the results obtained are presented in section 4, which are
divided in two different experiments: first, the performance of a proposed NN
model is evaluated (4.1), and then its scalability and flexibility is
evaluated on (4.2). Section 5 presents the conclusions and the discussion of
this work, respectively.

##  2 Materials

This section describes the sensors used and the dataset obtained from them.

###  2.1 Sensors

Three different variables were measured, namely, temperature, humidity and
pressure. Temperature and humidity sensors are integrated within a single
device, the Novasina nSens-HT-ENS. The humidity measurement range of this
device goes from $0\%$RH to $100\%$ RH, with a measurement accuracy of $0.5\%$
RH. Regarding the temperature measurement range, it goes from -20°to +80°,
with a measurement accuracy of 0.1°. The sensor used to measure the pressure
is the Novasina Pascal-ST-ZB. It is a differential pressure sensor based on
static pressure measurement. Regarding its measurement range, it goes from -50
Pa to 50 Pa, with a measurement accuracy of $<0.5\%$ fs (full scale).

###  2.2 Dataset

The dataset used to obtain the results presented in this work consists of
temperature, humidity and pressure values coming from the sensors installed
across different clean rooms devoted to drug production. Rooms are located in
two different floors of a building: the second floor and the basement. More
precisely, the data coming from floor 2 is described in Table 1. This dataset
contains data collected during a whole year in different rooms of the second
floor. Each room is endowed with only one sensor per category (one
temperature/humidity sensor and one pressure sensor). Thus, sensor redundancy
is not available. The sensors used are described in 2.1. The data was
collected at a sample per minute rate during the aforementioned period of
time.

A smaller dataset containing data from rooms located in the basement was also
collected. This data is shown in table 1. Conversely to the second floor, the
basement dataset is smaller due to the difficulties to access the data
generated within this floor. This dataset contains, approximately, one hundred
samples per sensor, coming from a total of eleven sensors per category.
Besides, samples are not equally distributed. The data points (also called
samples) were collected during different periods of the year.

Table 1: Dataset summary. Note that -1 floor refers to the basement. Floor | Sensor Type | Number of sensors | Number of samples
---|---|---|---
2 | Temperature | 17 | 8935200
2 | Humidity | 17 | 8935200
2 | Pressure | 24 | 12614400
-1 | Temperature | 11 | 1100000
-1 | Humidity | 11 | 1100000
-1 | Pressure | 11 | 1100000

Since the dataset obtained from the second floor is more complete, it was used
as the main dataset. Thus, experiments in subsection 4.1 and 4.2 were
conducted taking this dataset for training purposes. However, in subsection
4.2, the dataset coming from the basement was used for testing the adaptation
of the solution proposed to a different environment.

Regarding the dataset slicing, two approaches were defined depending on the
kind of experiments that were carried out. On the one hand, in subsection 4.1,
60% of the samples of the data coming from the second floor were used as
training set for each type of sensor, 10% were used as validation set and the
remaining 30% were used as testing set. On the other hand, in subsection 4.2,
70% of the samples were used for training the main Neural Network, 10% as
testing set and 20% for the application of the Transfer Learning. In this last
case not only data coming from the second floor was used, but also data coming
from the basement was also included in order to test the transfer learning
application from one floor to another.

The use of the different slices is more detailed in 3.

##  3 Methodology

In this section, the architecture of the implemented solution, which is
capable for detecting sensor uncalibration, is presented.

The approach is based on ML techniques and it follows the block diagram shown
in Figure 1. The architecture is based on an ANN which, during the training
phase, learns how the sensors behave. Then, in the inference step, the ANN is
able to predict what the measurements should be based on the input stimuli.
Eventually, these predictions are compared with the real measurements. If
there exists any discrepancy between the prediction and the real measurements,
an uncalibration is taking place.

Figure 1: Block diagram of the approach presented. Five main stages can be
observed: the measurement stage, the estimation of the set point, the
evaluation of the estimated set point by the NN and the generation of the
prediction and, finally, the computation of the error with respect to the real
value.

###  3.1 Preprocessing

The data used to both train and test the network have to be pre-processed.
This preprocessing, using filtering techniques, is a standard procedure to
detect anomalous values. These values could yield undesired effects on the
performance of the algorithm if they are not properly detected and treated.
Thus, this preprocessing phase can be divided into three stages, namely, a
threshold-based filtering, the Mahalanobis distance [15] and the Savitzky-
Golay filter [16].

The threshold-based filtering is based on the thresholds defined for the
warning and alarm systems in the rooms. This alarm system is designed to
detect extreme values in the incoming measurements. Thus, potential sudden
failures in the sensor systems can be detected. It has been observed that
these extreme values add noise to the input signal and, therefore, worsen the
NN learning. These thresholds were used during this research as a way to
discard all of those extreme values that were already detected by the system
and which could have a deep impact in the NN learning process. Thus, a
measured value is dropped whenever one of the aforementioned thresholds is
exceeded. In this way, only calibrated and non-extreme uncalibrated values
remain in the signal.

Next, a filter based on the Mahalanobis distance [17, 15, 18] is applied.
Mahalanobis distance is an outlier detection method specifically designed for
dealing with multidimensional distributions. More precisely, this stage
focuses on determining whether an observation including information from the
whole set of sensors (i.e a vector with as many entries as sensors) can be
considered as an outlier. Thus, in this stage, non-extreme values whose
relative occurrence is low would be classified as outliers.

Finally, the third stage is focused on reducing the noise in the incoming
signal. Different kind of filters such as Butterworth or wavelet filters could
have been applied in this stage. However, the incoming signal presents a wide
variety of different behaviors along the data collection period. Thus, a more
flexible noise filtering process is required. The Savitzky-Golay’s filter
[16], which computes local polynomial regressions providing new values for
each point, is used. This filtering process preserves essential properties of
the signal, such as maximum and minimum values or trends, while noise is
efficiently dropped out.

After the data obtained from the sensors is preprocessed, the average value of
the measurements, at a time, of all the sensors involved is used as input to a
ANN.

###  3.2 Neural network

The proposed ANN model is a Multi Layer Perceptron (MLP) that consists of an
input layer with one neuron, five hidden layers with three hundred neurons
each and an output layer with as many output neurons as sensors in used. After
training the network with the average value of the preprocessed data from the
sensors, the network is able to decompose the input stimulus into the
predicted value for each sensor, which is the expected behavior that each
sensor should have based on that input stimuli.

###  3.3 Error computation

This is the final step of the proposed approach. In this stage, the difference
between the real measured value obtained from the sensors and the ones
predicted is computed, evaluated and classified as normal or anomalous. This
evaluation is based on the work of [11], where Hierarchical Temporal Memory
[19] was used to predict the next measurement. That prediction was then
compared to the real measurement and residuals were computed.

In our approach, goodness of fit tests were applied on the errors obtained
during the training stage. Thus, the underlying distribution of the error
could be used to compute specific confidence intervals with the desired
signification for each sensor. On the one hand, the absolute mean error of the
sensors measurements properly fits a normal distribution. On the other hand,
the mean squared error fits an exponential distribution.

Therefore, the confidence intervals were generated by using the normal
distribution associated to the error of each sensor. This choice was based on
the fact that the absolute mean error was approximately equal to the
uncalibration that was taking place.

Finally, a rejection takes place whenever the value exceeds the bounds of the
confidence interval. The density of rejections is computed as the rolling
ratio between the number of rejections for a window of a fixed time length
(with 1440 minutes as default value, which corresponds to the number of
minutes in a day) and the length of the window.

The rejection density is the variable that triggered the uncalibration
warning. A uncalibration is said to take place whenever a threshold is
overcome during a specified amount of time.

It should be considered that the proposed architecture does not really know
what being calibrated means. The expected behaviour is that the NN detects
small differences between the current measurement and the one from which it
learnt. Thus, in the eventuality of a maintenance task in any sensor, the
architecture proposed might detect this new event as a potential
uncalibration. This uncalibration event will remain active until this new
condition of the sensor is re-learnt as ‘calibrated’ by the model.

The standard procedure to overcome this issue is to re-train the model under
these new conditions. Thus, this process would require a very high temporal
cost, since collecting data during almost a whole year would be necessary.
Because of this reason, transfer learning arises as a feasible solution. Using
transfer learning, a small amount of data is needed to teach the new condition
to the NN.

###  3.4 Transfer learning

Transfer learning [20] is a ML technique aimed at specializing ML models with
a minimum amount of data. This technique consists of two main steps. The first
step is training a model for a more general task conceptually related with the
target or specialized task. For instance, the target task could be identifying
a specific face in a picture and, in that case, the general task would be
training a model for general face identification. The second step, would be to
re-train a subset of layers of the model with a training dataset made out of
elements corresponding to the specialized task (a specific face in the example
presented).

In the case of this work, multiple scenarios can benefit from transfer
learning. These scenarios can be summarized as it follows: changing the ground
truth value for all sensors or for a subset of them, including new sensors and
adapting to some other environments with, possibly, lack of information.

As it was presented in subsection 2.2, a different slicing of the dataset was
used for transfer learning. Transfer learning is applied by dividing the
original dataset into three smaller subsets, namely, a training subset for the
original model, which will be called model A, a re-training subset for the
application of transfer learning (model A will be renamed to model B after
this process) and, finally, a testing subset for both models. Regarding the
training subset, it contains 70% of the datapoints. This stage corresponds to
the training for the general task. Regarding the subset devoted to be used for
re-training the model and, thus, for the application of transfer learning, it
contains the 20% of the datapoints. The amount of required datapoints can be
reduced to up to 2.5% of the total (10000 datapoints per sensor) in the case
of temperature and humidity. Finally, the testing subset includes 10% of the
datapoints. This subset was used to test the performance of the neural
network, both before and after the application of the transfer learning. A
constant offset was added in the last part of the test set in order to
simulate the recalibration due to a maintenance task. It is important to note
that this offset is not applied to the whole test set in order to properly see
the change of conditions and the associated evaluation by both, model A and
model B.

As mentioned, one of the main purposes of the application of transfer learning
is the reduced amount of necessary samples for training a model. This amount
is directly associated to the number of parameters to retrain, which, at the
same time, is intimately related to the number of layers to retrain. The
nature of the target task, i.e., how abstract it is, and where (in which
layers) the concepts learnt are located within the NN, are crucial factors
which will determine the amount of required data. During the current research,
transfer learning has been applied on the very last hidden layer of the NN,
i.e., only the weighs connecting the last hidden layer with the output layer
has been retrained. The obtained results suggest that training this subset of
weights was enough for the existing requirements.

##  4 Results

This section shows the results obtained using the methods presented in section
3. Two different experiments were conducted: (1) to evaluate the capability
for the detection of uncalibrations of the current approach through, which has
been called, performance evaluation experiments; (2) to test the flexibility
and scalability of the solution, which have been called transfer learning
experiments.

###  4.1 Performance evaluation experiments

These experiments were conducted to test the detection of uncalibrated sensors
using architecture presented. The dataset used is the one shown in 2.2.
Uncalibrations were introduced in the testing set as drifts of different
nature (e.g linear, exponential or logarithmic) across time and in different
temporal frames (for instance, at the beginning, in the middle and at the end
of the year).

An uncalibration is considered to be detected once the rejection density value
has overcome a predefined threshold, which depends on the sensor, and whose
average value is 0.8. Besides, this threshold has to be overcome during a
predefined period of time, which also depends on the sensor, and can be
defined as two weeks, since it is enough time for all sensors and it is a
reasonably low period of time for uncalibration detection. An example of the
application of this technique can be seen in figure 2 where rejection density
for both the upper and lower bound of the confidence interval are shown.

Figure 2: Uncalibration detection through rejection density. In the upper
subplot, the blue trace shows the error corresponding to a linear
uncalibration and the red dashed line stands for the error zero value. In the
middle subplot, the value of the rejection density for error values exceeding
the lower bounds of the confidence interval is shown. As it can be seen, the
maximum value for this rejection is, approximately, 0.3. In the lower subplot,
the value of the rejection density for error values exceeding the upper bounds
of the confidence interval is shown. As it can be seen, this density reaches
the maximum possible value from a critical point and during the whole
uncalibration phenomenon.

The second and third plots show the rejection density values. It can be seen
that these values keep low whenever a uncalibration is not taking place.
Conversely, when a uncalibration takes place, the rejection density reaches
values close to one.

Regarding the performance of the approach, uncalibrations were detected for
all different types of sensors. However, two different scenarios can be
distinguished, namely, a first optimistic scenario, where the uncalibration is
detected within the tolerance ranges defined by the quality requirements of
the pharma industry (i.e. $\pm 0.5$°,$\pm 3\%$ and $\pm 0.5$P in temperature,
humidity and pressure respectively), and a non-optimistic scenario, where the
uncalibration is detected over the tolerance range. For all those sensors that
are in the optimistic scenario, an estimation of the remaining time until the
tolerance range gets exceeded can be provided.

Thus, for the case of the temperature sensors, the proposed approach was
tested on the available data coming from temperature sensors over 17 rooms.
All uncalibrations were detected for all sensors. In the best scenario, the
architecture detected uncalibrations associated to deviations of 0.25°C
(tolerance 0.5°C). This accuracy was reached for 16 rooms out of 17 of the
second floor. In the worst scenario, the uncalibration was detected for
deviations of, approximately, 0.5°C.

For the humidity sensors case, the proposed method was tested on the available
data coming from humidity sensors over 17 clean rooms. All uncalibrations were
detected for all sensors. In the best scenario, the algorithm detects
uncalibrations associated to deviations of 2% (tolerance 3%). This accuracy
was reached for 14 rooms out of 17 of the second floor. In the worst scenario,
the uncalibration was detected for deviations close to 3%.

Finally, for the pressure sensors case, the proposed method was tested on the
available data coming from over 24 pressure sensors. All uncalibrations were
detected for all sensors. In the best scenario, the algorithm detected
uncalibrations associated to deviations under 0.5 Pascals (tolerance 0.5
Pascals). This accuracy was reached for 6 out of 20 rooms of the second floor.
In the worst scenario, the uncalibration was detected for deviations close to
1.5 Pascals.

Some other ANN architectures were tested during the development of the present
research. More precisely, the same experiments were conducted using a Wide &
Deep Neural Network for temperature and humidity sensors and Recursive Neural
Networks (RNN) for the pressure sensors. The architecture of the Wide & Deep
Neural Network was based on the MLP used on the implemented approach, but also
includes an input layer per sensor, so that the precise information coming
separately from each sensor can be included. This can be seen in figure 3.

Figure 3: Wide and Deep model diagram. In this figure, the tested architecture
of the wide and deep model is shown. It can be seen $n+1$ different input
layers; an input layer with one single neuron, which is connected to a hidden
block of five hidden layers with three hundred neurons each. This input layer
takes the estimated global set point as input. Next, $n$ input layers with one
single neuron each. Each input layer takes the measurements coming from the
different associated sensors as input. Finally, these input layers are
concatenated to the output of the hidden block and connected to the output
layer, which has a neuron per sensor.

Wide and Deep architecture yielded better results for most of the rooms.
However, for two of them, the uncalibrations were not detected. It has been
observed that, for these two rooms, how the sensor is supposed to behave under
a specific condition is not learnt by the neural network, and the identity
function is approximated instead. Thus, when an anomaly of any kind was
introduced on the input signal, this same value was retrieved by the NN and,
therefore, the uncalibration was not detected. It has been observed that the
worst results were obtained for these two rooms independently of the
architecture. This fact suggests that the local conditions inside these rooms
had a deep impact on the learning capabilities of the different models.

Regarding the RNN, Long-Short Term Memory (LSTM) neurons were used. In this
case, uncalibrations were not detected either. This was probably caused due to
the long-term nature of uncalibrations. Uncalibrations are long term
phenomenons, thus, the deviation associated to an uncalibration will be
observed during long periods of time. Hence, the deviation will be included as
relevant by LSTM neurons during the information inclusion stage, that is, when
the input gate evaluates what information is relevant and what information is
not. In this case, the learning of this long-term phenomenons by the NN
provoke the worst results.

###  4.2 Transfer learning experiments

These experiments were conducted to test the capability of the architecture
for learning a new condition with the less temporal cost. Thus, these
experiments were split in three categories, namely, change in global or
specific conditions, adding new sensors and adapting to some other
environments with lack of information. To evaluate transfer learning, the
error obtained for the main model (model A) and the retrained model (model B)
under the specific new condition were compared.

####  4.2.1 Change of conditions

These experiments were conducted to test the capability of the architecture
proposed for learning a new condition as the current calibration status for
one or more sensors.

In the case of the individual offset, in figure 4 two different phases of the
error obtained from model A can be observed. The first one corresponds to the
model evaluation under normal conditions. The second one corresponds to the
model evaluation under the generation of a specific offset. It can be seen how
the error is displaced, approximately, three units which is the same amount of
the generated offset.

Figure 4: In this figure, the error associated to a scheduled recalibration
can be seen. The blue trace shows the error and the red dashed line stands for
the error zero value. The recalibration takes place at the mid point of the
figure. Since transfer learning is not applied, once the recalibration takes
place, the error immediately increases.

However, regarding model B, it can be seen how the error returns to the same
range of values as previously. This behaviour suggests that the expected
results have been obtained. This experiment was reproduced for different rooms
and the different types of sensors, obtaining similar results. It can be
observed in figure 5.

Figure 5: In this figure, the blue trace shows the error associated to a
scheduled recalibration and the red dashed line stands for the error zero
value. In this case, transfer learning is applied, thus, once the
recalibration takes place, the error immediately decreases. It is important to
note that this model is not applied in the previous moments to the
recalibration task.

It should be noted that the minimum number of datapoints to properly use
transfer learning is around ten thousand samples. This number of samples
corresponds, approximately, to one week of data collection.

Regarding the generation of an offset for the whole set of sensors, the model
did not detect any change in its behavior. It suggests that, if the behaviour
of the whole set of sensors is equally changed, then, no effect on the joint
behaviour can be detected. Thus, no uncalibration is taking place from the
point of view of the model. This fact proposes an interesting question
regarding how the joint behaviour could be learnt.

####  4.2.2 Adding new sensors to the model

These experiments were conducted to determine the time that the architecture
takes to learn and include the behaviour of new sensors in the loop. These
tests were made using the humidity and temperature sensors.

A first dataset with data coming from 13 sensors out of seventeen was taken as
main set. For this set of sensors, $o_{1}=322677$ datapoints were available.
In this case, model A was trained with this dataset. Next, model B was re-
trained with data coming from both the initial set of thirteen sensors and a
set of four new sensors. The obtained results show that uncalibrations on the
new set of sensors can be properly detected by model B. Figure 6 shows how
uncalibrations are detected in the new sensors that were included for the
temperature case.

Figure 6: Uncalibration detection for new included sensors. In this figure, we
can see the results of the application of transfer learning in order to
include new sensors in the solution. The detection of an uncalibration for one
of those sensors is shown. In the upper subplot,the error corresponding to a
linear uncalibration can be seen. The red dashed line stands for the error
zero value. In the middle subplot, the value of the rejection density for
error values exceeding the lower bounds of the confidence interval is
presented. In this case, the maximum value for this rejection is approximately
zero. In the lower subplot, the value of the rejection density for error
values exceeding the upper bounds of the confidence interval is presented. As
it can be seen, this density reaches the maximum possible value from a
critical point and during the whole uncalibration phenomenon.

####  4.2.3 Adapting to some other environment with lack of information

These experiments were conducted to test if the architecture could be trained
with data coming from one specific location and then be implemented on a
different location after a transfer learning process.

In this case, model A was trained with data coming from the second floor with
17 sensors of both types, temperature and humidity. Then, model B was re-
trained with a small amount of data coming from the basement. Approximately,
80000 datapoints were used. The obtained results show that uncalibrations are
indeed detected by model B in this new environment. Furthermore, the problem
of the lack of data was solved by applying this method. Figure 7 shows an
example of how uncalibrations were detected by model B in this specific
context.

Figure 7: Uncalibration detection in the basement. In this figure, the
detection of an uncalibration taking place in the basement is shown. In this
case, transfer learning has been applied in order to re-train a model with a
small quantity of data coming from the basement. The model was previously
trained with data coming from the second floor. In the upper subplot, the
error corresponding to a linear uncalibration is shown in blue. The red dashed
line stands for the error zero value. In the middle subplot, the value of the
rejection density for error values exceeding the lower bounds of the
confidence interval is shown. In this case, the maximum value for this
rejection is 0.3. In the lower subplot, the value of the rejection density for
error values exceeding the upper bounds of the confidence interval is shown.
As it can be seen, this density reaches the maximum value starting from a
critical point and remains during the whole uncalibration phenomenon.

The obtained results suggest that transfer learning is a valid technique for
the adaptation of the solution to some other locations.

##  5 Discussion and Conclusions

Regarding the performance of the solution, it can be firstly concluded that
uncalibrations can indeed be detected even for values smaller than the defined
tolerance ranges for most cases in temperature and humidity sensors.

Secondly, it can be concluded that the application of transfer learning
yielded results that were coherent with the expected ones. Besides, these
results empowered the flexibility and scalability of the solution, enabling
the use of the model in a wide variety of contexts such as sensor addition,
integration within new environments and partially solving problems associated
to the lack of data.

Furthermore, transfer learning retrieved good results for very low amounts of
data. The minimum amount of data required is about 10000 samples which is the
data coming from one week of observations at a sampling rate of one sample per
minute. Indeed, the fact that only the last weight matrix was needed to re-
train means that the information associated to all these modifications were
learned on a very abstract level in the NN [21]. This fact suggest, on the one
hand, that offsets have a low impact on the joint dynamics of sensors. Thus,
the joint behavior of the whole set of sensors is ruled by much deeper
relationships than the addition of specific values. On the other hand, the
results obtained for pressure sensors, despite they are not under a common
constant condition, may be explained by this fact. The sudden changes in
pressure values could be then understood as an offset, what implies that the
model has a deep knowledge of the sensor behavior, although these offsets
reduce the accuracy on the detection of uncalibrations.

As a remarkable aspect of this research, it is important to note that the
solution has been trained and tested in a real context, with data coming
directly from production. Indeed, the solution presented in this work is
currently deployed in Azure [22], where the solution performs the
uncalibration detection directly from the information of the sensors from a
digital twin of the building. This digital twin contains the current status of
the physical sensors.

To the best of our knowledge, this is the first time that potential
uncalibrations of a set of sensors are online detected whenever the set point
is unknown. Furthermore, the proposed architecture can be extended by means of
transfer learning to a wide range of different fields.

Regarding the future work, different targets are considered: (1) extending the
results obtained after the application of transfer learning. Thus, we are
interested in training the solution with a generic set of sensors of one type
(such as temperature sensors) and then, testing the performance of the
solution on sensor systems of a different kind, for instance, engine sensor
systems, smart-city sensors systems and power plants sensor systems, among
others; (2) obtaining information about the learning of the solution through
the application of Explainable Artificial Intelligence (XAI). This could
provide useful information about the potential scalability and flexibility of
the solution and the different models to which it could be applied.

## Acknowledgement

This work was partially supported by the Altran Innovation Center for Advance
Manufacturing, the EU H2020 project Virtual IoT Maintenance System (VIMS Grant
ID: 878757), the Spanish grant (with support from the European Regional
Development Fund) MIND-ROB (PID2019-105556GB-C33) and by the EU H2020 project
CHIST-ERA SMALL (PCI2019-111841-2).

## References

  * Suk Kyu Lee and Kim [2017] M. B. Suk Kyu Lee, H. Kim,  Future of iot networks: A survey,  Applied Sciences 1072 (2017) 7(10).
  * Lammel [2015] G. Lammel,  The future of mems sensors in our connected world,  in: 2015 28th IEEE International Conference on Micro Electro Mechanical Systems (MEMS), 2015, pp. 61–64. doi:10.1109/MEMSYS.2015.7050886.
  * Alharbi and Soh [2019] N. Alharbi, B. Soh,  Roles and challenges of network sensors in smart cities,  IOP Conference Series: Earth and Environmental Science 322 (2019) 012002\. doi:10.1088/1755-1315/322/1/012002.
  * Singh et al. [2020] M. Singh, S. Sachan, A. Singh, K. K. Singh, Internet of Things in pharma industry: possibilities and challenges, 2020, pp. 195–216. doi:10.1016/B978-0-12-819593-2.00007-8.
  * Sharma et al. [2010] A. Sharma, L. Golubchik, R. Govindan,  Sensor faults detection methods and prevalence in real-world datasets,  TOSN 6 (2010).
  * Tobar et al. [2011] F. A. Tobar, L. Yacher, R. Paredes, M. E. Orchard,  Anomaly detection in power generation plants using similarity-based modeling and multivariate analysis,  in: Proceedings of the 2011 American Control Conference, IEEE, 2011, pp. 1940–1945.
  * Lamrini et al. [2018] B. Lamrini, A. Gjini, S. Daudin, P. Pratmarty, F. Armando, L. Travé-Massuyès,  Anomaly detection using similarity-based one-class svm for network traffic characterization.,  in: DX@ Safeprocess, 2018.
  * Feremans et al. [2019] L. Feremans, V. Vercruyssen, B. Cule, W. Meert, B. Goethals,  Pattern-based anomaly detection in mixed-type time series,  in: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, Springer, 2019, pp. 240–256.
  * Muruti et al. [2018] G. Muruti, F. A. Rahim, Z. bin Ibrahim,  A survey on anomalies detection techniques and measurement methods,  in: 2018 IEEE Conference on Application, Information and Network Security (AINS), 2018, pp. 81–86. doi:10.1109/AINS.2018.8631436.
  * Chalapathy and Chawla [2019] R. Chalapathy, S. Chawla,  Deep learning for anomaly detection: A survey,  CoRR abs/1901.03407 (2019). URL: http://arxiv.org/abs/1901.03407. arXiv:1901.03407.
  * Ahmad et al. [2017] S. Ahmad, A. Lavin, S. Purdy, Z. Agha,  Unsupervised real-time anomaly detection for streaming data,  Neurocomputing 262 (2017) 134–147.
  * Liu et al. [2019] J. Liu, M. Zhang, H. Wang, W. Zhao, Y. Liu,  Sensor fault detection and diagnosis method for ahu using 1-d cnn and clustering analysis,  Computational intelligence and neuroscience 2019 (2019).
  * Banjanovic-Mehmedovic et al. [2017] L. Banjanovic-Mehmedovic, A. Hajdarevic, M. Kantardzic, F. Mehmedovic, I. Dzananovic,  Neural network-based data-driven modelling of anomaly detection in thermal power plant,  Automatika: časopis za automatiku, mjerenje, elektroniku, računarstvo i komunikacije 58 (2017) 69–79.
  * Amiruddin et al. [2020] A. A. A. M. Amiruddin, H. Zabiri, S. A. A. Taqvi, L. D. Tufa,  Neural network applications in fault diagnosis and detection: an overview of implementations in engineering-related systems,  Neural Computing and Applications (2020) 1–26.
  * Mahalanobis [1936] P. C. Mahalanobis,  On the generalized distance in statistics,  National Institute of Science of India, 1936\.
  * Savitzky and Golay [1964] A. Savitzky, M. J. Golay,  Smoothing and differentiation of data by simplified least squares procedures.,  Analytical chemistry 36 (1964) 1627–1639.
  * Mahalanobis [1936] P. C. Mahalanobis,  On tests and measures of group divergence i. theoretical formulae.,  Proc. Nat. Inst. Sci. India (Calcutta) (1936) 49–55.
  * of Mathematics [2001] E. of Mathematics, Mahalanobis distance, 2001\. URL: https://encyclopediaofmath.org/index.php?title=Mahalanobis_distance.
  * Hawkins and Blakeslee [2004] J. Hawkins, S. Blakeslee, On intelligence, Macmillan, 2004\.
  * Bozinovski [2020] S. Bozinovski,  Reminder of the first paper on transfer learning in neural networks, 1976,  Informatica 44 (2020).
  * Géron [2019] A. Géron, Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build intelligent systems, O’Reilly Media, 2019\.
  * Chappell et al. [2010] D. Chappell, et al.,  Introducing the windows azure platform,  David Chappell & Associates White Paper (2010).

  *[ML]: Machine Learning
  *[AI]: Artificial Intelligence
  *[CNN]: Convolutional Neural Networks
  *[ANN]: Artificial Neural Network
  *[NN]: Neural Network
  *[MLP]: Multi Layer Perceptron
  *[RNN]: Recursive Neural Networks
  *[LSTM]: Long-Short Term Memory
  *[XAI]: Explainable Artificial Intelligence


$((q_{95\%}-q_{50\%})-(q_{50\%}-q_{5\%}))/(q_{95\%}-q_{5\%})$, where
$q_{5\%}$, $q_{50\%}$ and $q_{95\%}$ represent the $5\textsuperscript{th}$,
$50\textsuperscript{th}$, and $95\textsuperscript{th}$ percentiles of the
predictive densities, respectively. The green solid line represents the
computed skewness measure, while the black solid line marks the zero line.