# Towards Cross-Disaster Building Damage Assessment with Graph Convolutional
Networks ††thanks: Citation: Authors. Title. Pages…. DOI:000000/11111.

Ali Ismail
Department of Electrical and Computer Engineering
American University of Beirut
Beirut, Lebanon
<EMAIL_ADDRESS>
Mariette Awad
Department of Electrical and Computer Engineering
American University of Beirut
Beirut, Lebanon
<EMAIL_ADDRESS>
0000-0001-8636-1373 0000-0002-4815-6894

###### Abstract

In the aftermath of disasters, building damage maps are obtained using change
detection to plan rescue operations. Current convolutional neural network
approaches do not consider the similarities between neighboring buildings for
predicting the damage. We present a novel graph-based building damage
detection solution to capture these relationships. Our proposed model
architecture learns from both local and neighborhood features to predict
building damage. Specifically, we adopt the sample and aggregate graph
convolution strategy to learn aggregation functions that generalize to unseen
graphs which is essential for alleviating the time needed to obtain
predictions for new disasters. Our experiments on the xBD dataset and
comparisons with a classical convolutional neural network reveal that while
our approach is handicapped by class imbalance, it presents a promising and
distinct advantage when it comes to cross-disaster generalization.

_K_ eywords Change detection $\cdot$ Building damage $\cdot$ Urban data
science $\cdot$ Disaster informatics $\cdot$ Graph convolutional networks
$\cdot$ Siamese networks $\cdot$ Cross-domain generalization

##  1 Introduction

Change detection (CD) is a remote sensing technique used to detect changes
between images taken at different times. It is widely used in many
applications such as environmental [1] and urban monitoring [2]. Notably, CD
has been instrumental in detecting building damage in disasters’ aftermath [3,
4] Recent CD solutions have heavily relied on convolutional neural networks
(CNN) [5] including those aimed at building damage detection [6]. Disaster-
induced damage is a spatially spreading process therefore buildings in
proximity sustain damage similarly. However, this aspect is ignored by CNNs
which only consider the window defined by the convolution kernel. Graph
convolutional networks (GCN) are an emerging architecture that generalizes
stationary convolutions by taking into account non-euclidean relationships
defined by non-local graph connections [7]. This has lead to an improvement in
many computer vision applications [8, 9] including urban change detection
[10].

Additionally, effective planning of rescue operations requires damage mappings
to be obtained in a short time to minimize lives lost. However, data labeling
and machine learning training can take weeks. This has been mainly addressed
by pretrained generalizable models to infer on new disasters [11, 12, 13].
While GCNs initially were unable to be used on unseen contexts[7], the SAGE
(sample and aggregate) framework has been proposed to learn inductively on
mini-batches of graphs and can be used to predict on unseen graphs [14].

To address the gaps highlighted in the field, we approach building damage CD
as a node classification problem. Our graph construction uses building image
representations as node features and a geometric triangulation to connect each
building to its neighbors as opposed to the state-of-the-art which uses
summary superpixel features for nodes and connects only adjacent superpixels
[10]. To leverage local and neighborhood building features, we present a novel
hybrid architecture composed of a CNN backbone and a SAGE network. We test the
cross-domain generalization of this model using an experimental design similar
to [11]. However, our experiments consider different types of disasters across
multiple regions as opposed to earthquakes occurring in different locations.
We use the xBD dataset for training and validation [4].

The contributions of this paper are: 1) A novel graph node classification
approach to disaster building damage CD, 2) A novel model architecture
combining a CNN and a SAGE GCN to capture both local image features as well as
neighborhood building similarities, 3) A cross-domain generalization of the
proposed approach.

To the best of our knowledge, there is no prior work that studies the
transferability of Graph SAGE models for building damage CD.

The rest of this document is organized as follows. In Section 2, we explain
our methodology. Section 3 details the experiments performed, and result
interpretation while Section 4 concludes the work with future research
directions.

##  2 Methodology

###  2.1 Graph Data Formulation

We create an undirected acyclic graph where each building is connected to its
neighbors to exploit similarities between buildings in proximity. Each node
represents one building defined by a rectangular envelope enclosing the
building polygon. The node features are the concatenation of the pre and post
image crops within the envelope. To unify crop size and reduce the memory
footprint, we resize them to a width and height of $128\times 128$. The edges
are built using a Delaunay triangulation [15] based on the building envelope
centroid pixel coordinates. Edge weights are calculated according the
similarity measure used by [10] to build their adjacency matrix.

In the SAGE paradigm, subgraphs are sampled from the large graph to create
mini-batches. Since our graph is based on geographical proximity, every region
can be represented by a single graph. However, storing this large graph to
sample subgraphs proved to be too memory consuming due to the large number of
buildings and the high dimensionality of the node features. The xBD dataset is
divided into image chips each covering a certain geographical area and
containing a non uniform number of buildings. Therefore, we use this
predefined organization to build a subgraph out of each chip. Additionally,
image chips with only one building and image chips with only one annotated
building or with no labeled buildings were discarded. Figure 1 shows a a view
of an implemented subgraph.

Figure 1: Subgraph implementation.

###  2.2 Model Architecture

We present a novel architecture based on a GCN with a CNN backbone which is
necessary to reliably extract local image features while the downstream GCN is
used to aggregate them with CNN features of connected nodes. The backbone is a
Siamese ResNet34 network with the classification layers removed [16]
initialized with Imagenet weights. For the graph convolutions, we use the SAGE
operator [14]. The feature maps produced by the two ResNet streams are
subtracted and inputted to the SAGE network whose output is the damage
classification. All layers use the ReLu nonlinearity except for the final
output layer which uses a Softmax activation. Figure 2 shows the layout of the
proposed architecture.

Figure 2: Architecture of the proposed model.

##  3 Experiments and Results

###  3.1 Experimental Setup

The experiments were implemented using Python 3.8, Pytorch 1.7.1 and Pytorch
Geometric 1.7.0 [17]. The code111https://gitlab.com/awadailab/sage-project was
executed on a virtual machine running on 8 cores of an AMD EPYC 7551 32-Core
Processor with an Nvidia V100 32GB GPU. The Adam optimizer [18] was used with
the categorical cross entropy loss which included class weights to improve
performance on minority classes [6]. Additionally, we merged the "major
damage" and "destroyed" classes into a single class because of their ambiguity
which would have degraded model performance as done in [11]. The number of
neurons, number of layers, dropout rate, learning rate and batch size were
tuned empirically based on the performance on the test set. The resulting
hyperparameters are two layers with 32 neurons and a dropout rate of 0.5. The
batch size was set to 256 and the learning rate to 0.0003. For evaluation, we
compute the accuracy, macro-averaged F1, weighted F1 and AUC scores.

###  3.2 Experimental Procedures and Results

At first, we train and test on the same disaster type (fire) and we use two
fire disasters to make the training set sufficiently large. The following
experiments test transfer from one disaster (flooding) and from multiple
disasters (flooding, tornado and volcano) to a different disaster (fire).
Lastly, we explore leaking a small number of test samples. As a benchmark, we
adopt a traditional Siamese CNN identical to the Graph SAGE architecture
except that the SAGE layers are replaced with fully connected layers. This
comparison, shown in Table 1, allows us to isolate the contribution of the
graph convolutions.

Table 1: Comparison between the Graph SAGE and the Siamese CNN.

Train | Test/Hold | Siamese CNN | Graph SAGE
---|---|---|---
Socal + Portugal fire | Socal |  | Acc | Macro F1 | Weighted F1 | AUC |  | Acc | Macro F1 | Weighted F1 | AUC
Train | 0.8746 | 0.5994 | 0.9036 | 0.9648 | Train | 0.9110 | 0.5642 | 0.9205 | 0.9347
Test | 0.8439 | 0.5267 | 0.8696 | 0.8606 | Test | 0.8672 | 0.5042 | 0.8768 | 0.8567
Hold | 0.8581 | 0.5449 | 0.8787 | 0.8609 | Hold | 0.8870 | 0.5315 | 0.8902 | 0.8374
Nepal flooding | Socal |  | Acc | Macro F1 | Weighted F1 | AUC |  | Acc | Macro F1 | Weighted F1 | AUC
Train | 0.6921 | 0.6342 | 0.7275 | 0.8764 | Train | 0.5092 | 0.3421 | 0.5487 | 0.5314
Test | 0.4609 | 0.3223 | 0.5911 | 0.6229 | Test | 0.5697 | 0.3245 | 0.6525 | 0.6121
Hold | 0.4356 | 0.315 | 0.5708 | 0.6601 | Hold | 0.6075 | 0.3275 | 0.6727 | 0.4999
Nepal flooding + Joplin tornado \+ Puna volcano | Socal |  | Acc | Macro F1 | Weighted F1 | AUC |  | Acc | Macro F1 | Weighted F1 | AUC
Train | 0.9283 | 0.9008 | 0.9306 | 0.9877 | Train | 0.7042 | 0.2824 | 0.5870 | 0.5199
Test | 0.599 | 0.3832 | 0.6928 | 0.7265 | Test | 0.8744 | 0.3110 | 0.8157 | 0.6319
Hold | 0.6082 | 0.3913 | 0.6966 | 0.6597 | Hold | 0.8542 | 0.3072 | 0.7873 | 0.6445
Nepal flooding + Joplin tornado \+ Puna volcano + 10% Socal Fire | Socal |  | Acc | Macro F1 | Weighted F1 | AUC |  | Acc | Macro F1 | Weighted F1 | AUC
Train | 0.9666 | 0.9522 | 0.9671 | 0.9963 | Train | 0.7094 | 0.2873 | 0.5951 | 0.5394
Test | 0.6659 | 0.4205 | 0.7418 | 0.7173 | Test | 0.8591 | 0.3172 | 0.7959 | 0.6765
Hold | 0.6668 | 0.4308 | 0.7405 | 0.6853 | Hold | 0.8552 | 0.3131 | 0.7899 | 0.6665

The Graph SAGE outperforms the Siamese CNN only on $20$ out of $48$ total
occasions in terms of raw performance. The likely reason is that the Graph
SAGE is being more severely degraded by the class imbalance. The Graph SAGE
seems to outperform the Siamese CNN in terms of accuracy and weighted F1 score
with an average increase of $14.74\%$ and $6.34\%$ on the hold set
respectively. These two metrics ignore the impact of class imbalance. On the
other hand, it falls behind most of the time in terms of macro-averaged F1 and
underperforms systematically in terms of AUC with an average difference of
$5.07\%$ and $5.44\%$ on the hold set respectively. Another possible reason is
that given the geographical distribution of the samples, a large majority of
subgraphs had no damaged buildings.

In terms of cross-domain generalization, the Graph SAGE offers a more
satisfactory behavior. We demonstrate this by computing the difference between
the performance on the training set and the hold set (Figure 3). For the first
experiment, the Siamese CNN presents a lesser drop in performance. However,
for subsequent experiments which are cross-domain, the difference for Graph
SAGE is negative in most cases which means that performance on the hold set
was higher than performance on the training set. Furthermore, the negative
difference increases in magnitude as we move from left to right. Meaning that
generalization improved when augmenting the training data with more diverse
disasters. This implies that larger training data improves generalization to
new domains. When the training data is supplemented with a small number of
test samples, we achieve the best performance. On the other hand, the Siamese
CNN shows a consistent positive generalization gap which is likely due to
overfitting the training data.

Figure 3: Difference between training and hold scores for the Siamese CNN and
Graph SAGE. The x-axis indexes indicate the different experiments with
different train - test configurations: 1) Fire - Fire, 2) Flooding - Fire, 3)
Flooding + Tornado + Volcano - Fire and 4) Flooding + Tornado + Volcano + 10%
Fire - Fire.

##  4 Conclusion

In the wake of disasters, building damage mappings provide essential
intelligence for planning effective rescue and relief efforts. And the sooner
this information is acquired, the more likely it is to rescue injured and
trapped victims.

Like most geographical phenomena, disaster damage propagates such that
buildings in proximity sustain damage in similar patterns. In this work, we
formulated a graph-based building damage CD solution that captures both local
patterns and building neighborhood relationships which are not captured by
current approaches. Additionally, we focus on the cross-disaster
generalization of our model which would enable our solution to be used to
infer on new disasters without training thus accelerating the process.
Although our results are not conclusive in terms of raw performance when
compared to a classical CNN model, we note an advantage in generalizing to
unseen and different disasters.

Based on the insights gained throughout this work, our future research
includes devising more optimal subgraph sampling strategies and using graph-
specific data balancing mechanisms to address the shortcomings of our
approach.

## References

  * [1] Ahram Song, Yeji Kim, and Yongil Kim.  Change Detection of Surface Water in Remote Sensing Images Based on Fully Convolutional Network.  Journal of Coastal Research, 91(sp1):426, August 2019.
  * [2] Chi Zhang, Shiqing Wei, Shunping Ji, and Meng Lu.  Detecting Large-Scale Urban Land Cover Changes from Very High Resolution Remote Sensing Images Using CNN-Based Classification.  ISPRS International Journal of Geo-Information, 8(4):189, April 2019\.
  * [3] Stefan A. Robila.  Use of Remote Sensing Applications and its Implications to the Society.  In 2006 IEEE International Symposium on Technology and Society, pages 1–6, Queens, NY, USA, June 2006. IEEE.
  * [4] Ritwik Gupta, Richard Hosfelt, Sandra Sajeev, Nirav Patel, Bryce Goodman, Jigar Doshi, Eric Heim, Howie Choset, and Matthew Gaston.  xBD: A Dataset for Assessing Building Damage from Satellite Imagery.  arXiv:1911.09296 [cs], November 2019.
  * [5] Lazhar Khelifi and Max Mignotte.  Deep Learning for Change Detection in Remote Sensing Images: Comprehensive Review and Meta-Analysis.  arXiv:2006.05612 [cs], June 2020.
  * [6] Bradley J. Wheeler and Hassan A. Karimi.  Deep Learning-Enabled Semantic Inference of Individual Building Damage Magnitude from Satellite Images.  Algorithms, 13(8):195, August 2020.
  * [7] Thomas N. Kipf and Max Welling.  Semi-Supervised Classification with Graph Convolutional Networks.  arXiv:1609.02907 [cs, stat], February 2017.
  * [8] Lichao Mou, Xiaoqiang Lu, Xuelong Li, and Xiao Xiang Zhu.  Nonlocal Graph Convolutional Networks for Hyperspectral Image Classification.  IEEE Transactions on Geoscience and Remote Sensing, pages 1–12, 2020.
  * [9] Danfeng Hong, Lianru Gao, Jing Yao, Bing Zhang, Antonio Plaza, and Jocelyn Chanussot.  Graph Convolutional Networks for Hyperspectral Image Classification.  IEEE Transactions on Geoscience and Remote Sensing, pages 1–13, 2020.
  * [10] Sudipan Saha, Francesca Bovolo, and Lorenzo Bruzzone.  Semisupervised Change Detection Using Graph Convolutional Network.  IEEE Geoscience and Remote Sensing Letters, 18(4):607 – 611, 2020\.
  * [11] Joseph Z. Xu, Wenhan Lu, Zebo Li, Pranav Khaitan, and Valeriya Zaytseva.  Building Damage Detection in Satellite Imagery Using Convolutional Neural Networks.  arXiv:1910.06444 [cs, eess, stat], October 2019.
  * [12] Vitus Benson and Alexander Ecker.  Assessing out-of-domain generalization for robust building damage detection.  arXiv:2011.10328 [cs], November 2020.
  * [13] Wanting Yang, Xianfeng Zhang, and Peng Luo.  Transferability of Convolutional Neural Network Models for Identifying Damaged Buildings Due to Earthquake.  Remote Sensing, 13(3):504, January 2021.
  * [14] William L. Hamilton, Rex Ying, and Jure Leskovec.  Inductive Representation Learning on Large Graphs.  arXiv:1706.02216 [cs, stat], September 2018.
  * [15] D. T. Lee and B. J. Schachter.  Two algorithms for constructing a Delaunay triangulation.  International Journal of Computer & Information Sciences, 9(3):219–242, June 1980.
  * [16] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition, 2015.
  * [17] Matthias Fey and Jan Eric Lenssen.  Fast Graph Representation Learning with PyTorch Geometric.  arXiv:1903.02428 [cs, stat], April 2019.
  * [18] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.


RNN

    _Recurrent Neural Networks_
LSTM

    _Long-Short Term Memory Networks_
NODE

    _Neural Ordinary Differential Equations_
ODE

    _Ordinary Differential Equations_
RKNN

    _Runge-Kutta Neural Networks_
LTC

    _Liquid Time Constant Neural Networks_
NN

    _Neural Network_
ReLU

    _Rectified Linear Function_
DAE

    _Differential Algebraic Equations_
MSE

    _Mean Squared Error_
TPE

    _Tree-Structured Parzen Estimator_
TSEM

    _Truncated Simulation Error Minimization_
BLA

    _Best Linear Approximation_
RMSE

    _Root Mean Squared Error_
CTS

    _Cascaded Tank System_
BLA

    _Best Linear Approximation_
NLSS

    _Nonlinear State Space Model_
CT

    _Continuous-Time_
DT

    _Discrete-Time_
DFT

    _Discrete Fourier Transform_
SDN

    _State Derivative Normalization_
PEM

    _Prediction Error Minimization_
††thanks: The project RACKET supported this research under grant 01IW20009 by
the German Federal Ministry of Education and Research.

This work has been submitted to IFAC for possible publication.

# State Derivative Normalization for Continuous-Time Deep Neural Networks

Jonas Weigand  Gerben I. Beintema  Jonas Ulmen  Daniel Görges  Roland Tóth
Maarten Schoukens  Martin Ruskowski  Chair of Machine Tools and Control
Systems, RPTU Kaiserslautern, and the German Research Center for Artificial
Intelligence, Kaiserslautern, Germany (e-mail<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Control Systems (CS) Group at the Department of
Electrical Engineering, Eindhoven University of Technology, Netherlands. R.
Tóth is also affiliated to the Systems and Control Laboratory at the Institute
for Computer Science and Control, Budapest, Hungary (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>Institute for
Electromobility, RPTU Kaiserslautern, Germany (e-mail<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

The importance of proper data normalization for deep neural networks is well
known. However, in continuous-time state-space model estimation, it has been
observed that improper normalization of either the hidden state or hidden
state derivative of the model estimate, or even of the time interval can lead
to numerical and optimization challenges with deep learning based methods.
This results in a reduced model quality. In this contribution, we show that
these three normalization tasks are inherently coupled. Due to the existence
of this coupling, we propose a solution to all three normalization challenges
by introducing a normalization constant at the state derivative level. We show
that the appropriate choice of the normalization constant is related to the
dynamics of the to-be-identified system and we derive multiple methods of
obtaining an effective normalization constant. We compare and discuss all the
normalization strategies on a benchmark problem based on experimental data
from a cascaded tanks system and compare our results with other methods of the
identification literature.

###### keywords:

System Identification, Continuous-Time, Neural Ordinary Differential
Equations, Nonlinear State-Space

##  1 Introduction

Machine learning is a powerful tool for time series modeling and system
identification (Schoukens and Ljung, 2019). Models can be categorized into
_Discrete-Time_ (DT) models and _Continuous-Time_ (CT) models such as _Neural
Ordinary Differential Equations_ (NODE) (Chen et al., 2018), _Runge-Kutta
Neural Networks_ (RKNN) (Wang and Lin, 1998), Deep Encoder Networks (Beintema
et al., 2023) or _Liquid Time Constant Neural Networks_ (LTC) (Hasani, 2020).
In control engineering, CT models are often preferred, as they are closely
related to physical models, they are often found more attractive for designing
controllers since shaping performance with such models is more intuitive, and
they can be used under irregular sampling.

Data normalization is a key data preprocessing step in modern machine learning
approaches. The core idea of normalization is to map the input and output data
to a numerically favorable range (e.g. zero-mean and standard deviation of 1)
(Ba et al., 2016). This is important to (_i_) comply with the underlying
assumption in many parameter initialization strategies for _Neural Networks_
(e.g., Xavier initialization (Glorot and Bengio, 2010)) and to (_ii_)
numerically improve the effective numerical range of the gradients for the
nonlinear activation functions (Ioffe and Szegedy, 2015).

Prior research demonstrates the significant enhancement in model performance
through _State Derivative Normalization_ (SDN) as evidenced by (Weigand et
al., 2021; Beintema et al., 2023). Despite these advancements, a comprehensive
exploration and analysis of SDN’s mechanisms remain absent from published
literature. In addition, practical methodologies to estimate the normalization
factor are missing in the literature.

  * •

Consequently, this paper analyses the state, the state derivative, and the
time normalization as well as their coupling. Additional graphical
interpretation and the connection to _Ordinary Differential Equations_ (ODE)
solvers are presented.

  * •

Furthermore, three approaches are proposed to estimate this normalization
factor and these approaches are compared on a well-studied benchmark example.

The remainder of the paper is structured as follows. Section 2 introduces the
considered identification problem. We present multiple interpretations of the
use of a normalization constant in Section 3. Section 4 defines the criteria
for an effective normalization and provides three methods for practical
implementation. Section 5 applies the method to a real-world benchmark
problem, and in Section 6 concluding remarks on the established results are
made.

##  2 The Identification Problem

###  2.1 Data-generating System

Consider a nonlinear system that is represented by a nonlinear state-space
representation

$\displaystyle\dot{x}(t)$ $\displaystyle=f\left(x(t),u(t)\right)$ (1a)
$\displaystyle y_{k}$ $\displaystyle=h(x_{k},u_{k})+w_{k},$ (1b)
$\displaystyle x(0)$ $\displaystyle=x_{0},$ (1c)

where $t\in\mathbb{R}_{\geq 0}$ is the continuous time,
$x(t)\in\mathbb{R}^{n_{\mathrm{x}}}$ is the state associated with (1) with
$n_{\mathrm{x}}\in\mathbb{N}$, $x_{0}\in\mathbb{R}^{n_{\mathrm{x}}}$ being the
initial state, $u(t)\in\mathbb{R}^{n_{\mathrm{u}}}$ is the exogenous input,
and $f:\mathbb{R}^{n_{\mathrm{x}}\times
n_{\mathrm{u}}}\rightarrow\mathbb{R}^{n_{\mathrm{x}}}$ is locally Lipschitz
continuous ensuring that all solutions of $x$ are forward complete. With
respect to the output equation, ${y}(k)\in\mathbb{R}^{n_{\mathrm{y}}}$ stands
for the sampled system output at time moments $kT_{\mathrm{s}}$ with
$T_{\mathrm{s}}>0$ being the sampling interval and $k\in\mathbb{Z}_{\geq 0}$,
$x_{k}=x(kT_{\mathrm{s}})$ and $u_{k}=u(kT_{\mathrm{s}})$, while
$w_{k}\in\mathbb{R}^{n_{\mathrm{y}}}$ is an i.i.d. zero-mean white noise
process with finite variance $\Sigma_{w}$. Such a hybrid process and noise
model is often used in CT nonlinear system identification (Schoukens and
Ljung, 2019).

###  2.2 Identification

For the considered system class, the CT model identification problem can be
expressed for a given data set of measurements

$D_{N}=\\{(u_{0},y_{0}),(u_{1},y_{1}),...,(u_{N-1},y_{N-1})\\},$ (2)

generated by (1), with unknown $w_{k}$, $x(t)$, $\dot{x}(t)$, and initial
state $x_{0}$, as the following optimization problem (a.k.a. $\ell_{2}$
simulation loss minimization)

$\displaystyle\min_{\theta,\hat{x}_{0}}$
$\displaystyle\frac{1}{N}\sum_{k=0}^{N-1}\|y_{k}-\hat{y}_{k}\|_{2}^{2},$ (3)
s.t. $\displaystyle\hat{y}_{k}=h_{\theta}(\hat{x}(kT_{\mathrm{s}})),$
$\displaystyle\dot{\hat{x}}(t)=f_{\theta}(\hat{x}(t),u(t)),$

where $\hat{x}(t)\in\mathbb{R}^{n_{\mathrm{x}}}$ is the model state,
$h_{\theta}$ and $f_{\theta}$ are the output and state-derivative functions
parameterized by $\theta\in\mathbb{R}^{n_{\theta}}$. These two functions are
represented by multi-layer feedforward NN’s within this paper. However, the
presented results also hold when other function approximators such as
polynomials are considered.

##  3 Interpretation of the Normalization

_State Derivative Normalization_ (SDN) The idea is the introduction of a
normalization factor $\tau$, which scales the state-derivative equation of the
state-space model during training:

$\dot{\hat{x}}(t)=\frac{1}{\tau}f_{{}_{\mathrm{NN}}}\\!\left(\hat{x}(t),u(t)\right)$
(4)

The normalization can be implemented as a scalar normalization
$\tau\in\mathbb{R}_{>0}$, or as a vector
$\tau\in\mathbb{R}^{n_{\mathrm{x}}}_{>0}$ corresponding to elements of the
output of the hidden state network $f_{{}_{\mathrm{NN}}}$. In the latter case,
(4) is defined as
$\dot{\hat{x}}(t)=\mathrm{diag}(\frac{1}{\tau_{1}},\ldots,\frac{1}{\tau_{n_{\mathrm{x}}}})f_{\theta}(\hat{x}(t),u(t))$.
This linear scaling can also be interpreted as changing the weight
initialization of the output layer of $f_{{}_{\mathrm{NN}}}(\cdot,\cdot)$.

It is clear that $\tau$ in (4) scales the hidden state derivative. Hence, it
can be used to normalize the state derivatives. Next, we show that the
derivative normalization is inherently coupled to the scaling of the state $x$
and the scaling of the time $t$.

State Normalization. It is also possible to rewrite (4) as:

$\displaystyle{\frac{d(\tau\hat{x}(t))}{dt}}$
$\displaystyle={f_{{}_{\mathrm{NN}}}}(\hat{x}(t),u(t))$ (5a)
$\displaystyle\tilde{x}(t)$ $\displaystyle\triangleq\tau\hat{x}(t)$ (5b)
$\displaystyle{\frac{d\tilde{x}(t)}{dt}}$
$\displaystyle={f_{{}_{\mathrm{NN}}}}(\tilde{x}(t)/\tau,u(t))$ (5c)

In this way, normalization can be interpreted as the normalization of the
magnitude of the state.

Time-Based Normalization. We can rewrite (4) as

$\displaystyle{\frac{d\hat{x}(t)}{d(t/\tau)}}$
$\displaystyle={f_{{}_{\mathrm{NN}}}}(\hat{x}(t),u(t)),$ (6a)
$\displaystyle\tilde{t}$ $\displaystyle\triangleq t/\tau,$ (6b)
$\displaystyle{\frac{d\hat{x}(\tilde{t}\tau)}{d\tilde{t}}}$
$\displaystyle={f_{{}_{\mathrm{NN}}}}(\hat{x}(\tilde{t}\tau),u(\tilde{t}\tau)),$
(6c)

which suggests that the normalization can be viewed as a scaling of time by a
factor $\tau$.

This time rescaling can also be viewed as changing the effective integration
length in ODE solvers. For instance,

$\displaystyle\hat{x}(t)$
$\displaystyle=\hat{x}(0)+\int_{0}^{t}\frac{1}{\tau}f_{{}_{\mathrm{NN}}}(\hat{x}(t^{\prime}),u(t^{\prime}))\
dt^{\prime}$ (7a)
$\displaystyle=\hat{x}(0)+\int_{0}^{t/\tau}\\!\\!\\!f_{{}_{\mathrm{NN}}}(\hat{x}(\tau\tilde{t}^{\prime}),u(\tau\tilde{t}^{\prime}))\
d\tilde{t}^{\prime}$ (7b)

using the substitution of $\tilde{t}^{\prime}=t^{\prime}/\tau$. This
integration by substitution shows that the normalization can be viewed as
rescaling the effective integration time from $t$ to $t/\tau$. This can be
translated directly to a wide range of numerical integration schemes by
rescaling the integration length and step size simultaneously.

Graphical intuition about the time rescaling can be gained from output data of
the _Cascaded Tank System_ (CTS) benchmark in Fig. 1 (Schoukens et al., 2016).
The measured data (black) of the water level is given at a sample rate of
$T_{\mathrm{s}}=4.0\,\mathrm{s}$, the scaled data (red) is transferred to a
sample time of $T_{m}=T_{\mathrm{s}}/\tau=1.0$, while the number of time steps
remains unchanged. The original data is transferred to the scaled time domain,
where the model is trained and evaluated, and the results are transferred back
to the original time grid.

$0$$200$$400$$600$$800$$1{,}000$$1{,}200$$1{,}400$$4$$6$$8$$10$Time
(s)/UnitlessWaterlevel (V)
---
Figure 1: Time scaling of the output measurement of the _Cascaded Tank System_
(CTS) benchmark (Schoukens et al., 2016). Black: Original measurement (unit
seconds). Red: Time-domain scaled data (unitless).

It is known that the stiffness of an ODE system is subject to temporal
variation. However, employing a constant (linear) normalization across all
stiffness regions is essential from a practical model inference point of view.
Moreover, averaging across multiple regions increases the robustness of the
estimation by providing a broader range of suitable normalization factors. Our
subsequent experiments demonstrate that there exists a broad valley of good
normalization factors (Fig. 3). Remarkably, all normalization factors within
this range lead to high-performance models.

The key takeaway of representations (4), (6) and (5) is that SDN collectively
influences the hidden state, the hidden state derivative, and the scaling of
time $t$.

##  4 Estimation of the Normalization Factor

Since the system dynamics are unknown a priori, selecting a normalization
factor $\tau$ that results in a favorable numerical performance of the
training algorithm is challenging a priori, it has to be estimated during
training, from the available data, or by heuristics. We propose three methods:

  * •

make the normalization factor a trainable parameter,

  * •

estimate $\tau$ using cross-validation,

  * •

use a heuristic-based approach starting from a linear approximation of the
system.

Trainable Parameter. The normalization factor $\tau$ can be added to the set
of parameters present in the optimization problem (3). This results in a
simultaneous optimization of $\tau$, $\theta$, and $\hat{x}_{0}$ under the
$\ell_{2}$ simulation loss.

In terms of implementation, $\tau>0$ should hold to fulfil its role as a
magnitude normalization. As most NN libraries do not provide constrained
optimization one can implement

$\hat{\tau}=\mathrm{max}(\epsilon,\tau)$ (8)

with a small positive number $\epsilon\in\mathbb{R}_{>0}$. Empirical evidence
demonstrates that the trainable normalization quickly reaches a stationary
value.

From a pure mathematical perspective, making $\tau$ a trainable parameter may
seem to have little impact because the expressive power of the NN could easily
compensate for it. Practically, the optimizer would need to adapt a large
quantity of parameters instead of a single normalization factor. Additionally,
traditional NN initialization methods assume inputs and outputs are within a
specific range. We consolidated the normalization dependency into a single
trainable parameter $\tau$, making the training process simpler and more
robust, as demonstrated in our experiments.

Overall, this method offers high adaptability across different data sets and
NN configurations, an easy implementation, and demands only little additional
computational load.

Cross-validation. The normalization factor $\tau$ can be considered as a
hyperparameter and estimated using cross-validation methods. For example, (3)
is optimized using a training data set for a chosen set of values of $\tau$, a
validation data set is used subsequently to determine the best performing
$\tau$ value. This method is universally applicable to nonlinear systems, is
easy to implement if a hyperparameter tuning is already available, and, given
sufficient function evaluations, ensures global optimality.

Heuristic approach based on an approximate model. This heuristic is derived
based on a guarantee of the existence of a normalized model formalized by
Beintema et al. (2023) which states;

###### Theorem 1

Given a input trajectory $u(t)$, a non-constant state trajectory $x(t)$ which
satisfies $\dot{x}(t)=f(x(t),u(t))$ as in Eq. 1 for all $t\in[0,L]$ than there
exists a $\tau\in\mathbb{R}^{+}$ and a scalar state transformation
$\gamma\hat{x}(t)=x(t)$ such that both the model state trajectory $\hat{x}(t)$
and model derivative $f_{\textrm{NN}}$ given by
$\dot{\hat{x}}(t)=\tfrac{1}{\tau}f_{\textrm{NN}}(\hat{x}(t),u(t))$ as in Eq. 4
are normalized as

$\displaystyle\mathrm{var}(\hat{x})$
$\displaystyle\triangleq\frac{1}{L}\int_{0}^{L}\frac{1}{n_{\mathrm{x}}}\|\hat{x}(t)\|_{2}^{2}dt=1,$
(9a) $\displaystyle\mathrm{var}(f_{\textrm{NN}}(\hat{x},u))$
$\displaystyle=1.$ (9b)

With

$\displaystyle\gamma$ $\displaystyle=\sqrt{\mathrm{var}(x)}$ (10)
$\displaystyle\tau$
$\displaystyle=\sqrt{\mathrm{var}(x)/\mathrm{var}(\dot{x})}$ (11)

the normalization conditions are satisfied, as shown below

$\displaystyle\mathrm{var}(\hat{x})$
$\displaystyle=\mathrm{var}(x/\gamma)=\mathrm{var}(x)/\gamma^{2}=1,$
$\displaystyle\mathrm{var}(f_{\textrm{NN}}(\hat{x},u))$
$\displaystyle=\mathrm{var}(\tau\dot{\hat{x}})=\mathrm{var}(\tau\dot{x}/\gamma)=\tau^{2}/\gamma^{2}\mathrm{var}(\dot{x})=1.$

We base our heuristic on the relationship of Eq. 11. We make two alterations
(i) the integral over time can be approximated by a summation over the time
samples to make this tractable and (ii) in Eq. 11 we use an approximate model
instead of the system equations.

We have observed that using a linear approximate model based on the _Best
Linear Approximation_ (BLA) of a nonlinear system can get a sufficiently
accurate estimate of the optimal $\tau$. The BLA offers a linear approximation
of a nonlinear system, best in the mean-squared-error sense, based on measured
input-output data (Pintelon and Schoukens, 2012). Estimating a BLA of a
potentially nonlinear system is commonly done using _Prediction Error
Minimization_ (PEM), minimizing a least squares prediction error. Using an
estimated BLA, an estimate of the state and state-derivatives,
$x_{\mathrm{BLA}}$ and $\dot{x}_{\mathrm{BLA}}$, can be obtained to compute
$\tau$ through (11):

$\tau_{\mathrm{BLA}}=\sqrt{\mathrm{var}(x_{\mathrm{BLA}})/\mathrm{var}(\dot{x}_{\mathrm{BLA}})}.$
(12)

This method is especially valuable for dynamic systems that can be reasonably
well represented by a linear model for the considered range of excitation.
Also note that any linear transformation that can be present in the BLA state
and state derivative estimate does not affect the resulting $\tau$ estimate in
(12). Additionally, besides normalization, the BLA can be used as a good NN
weight initialization candidate (Schoukens and Tóth, 2020).

Frequency-domain analysis of this approach offers additional insight into the
optimal choice of the normalization factor. Consider that the input signal has
a periodicity of $N$ samples, then one can decompose the input into its
_Discrete Fourier Transform_ (DFT) components as:

$\displaystyle u_{k}=\frac{1}{N}\sum_{m=0}^{N-1}U_{m}e^{j\frac{2\pi}{N}mk}$
(13) $\displaystyle U_{m}=\sum_{k=0}^{N-1}u_{k}e^{-j\frac{2\pi}{N}mk}$ (14)

Both the variance of the state and state-derivative can be expressed using
these DFT components as:

$\displaystyle\text{var}(x)$
$\displaystyle=\frac{1}{L}\int_{0}^{L}\|x(t)\|_{2}^{2}dt\sim\sum_{m=-\infty}^{\infty}\|X_{m}\|_{2}^{2}$
(15a) $\displaystyle\text{var}(\dot{x})$
$\displaystyle=\frac{1}{L}\int_{0}^{L}\|\dot{x}(t)\|_{2}^{2}dt\sim\sum_{m=-\infty}^{\infty}\omega_{m}^{2}\|X_{m}\|_{2}^{2}$
(15b)

where $\omega_{m}=2\pi\frac{m}{NT_{s}}$.

Since we are considering the BLA, the Fourier components of the state can be
written in terms of the input-to-state frequency response function
$X_{m}=G(j\omega_{m})U_{m}$, assuming a single input system for simplicity.
Therefore, substituting this and (15) into (12) results in

$\tau_{\mathrm{BLA}}=\sqrt{\frac{\sum_{m=-\infty}^{\infty}U_{m}^{H}G^{H}\\!(j\omega_{m})G(j\omega_{m})U_{m}}{\sum_{m=-\infty}^{\infty}\omega_{m}^{2}U_{m}^{H}G^{H}\\!(j\omega_{m})G(j\omega_{m})U_{m}}},$
(16)

where $\cdot^{H}$ denotes the Hermitian operation. Hence, if the signal $u(t)$
only consists of a single sine wave (i.e. only one nonzero $U_{m}$) then
$\tau_{\mathrm{BLA}}=1/\omega_{m}$. Furthermore, if $u(t)$ or $G(j\omega_{m})$
has a finite bandwidth bounded by $\Omega$ then $\tau_{\mathrm{BLA}}\geq
1/\Omega$.

##  5 Experimental Results

Benchmark. We apply the proposed methods to the CTS benchmark (Schoukens et
al., 2016). The setup is depicted in Fig. 2. It consists of two vertically
mounted tanks, where the upper one has a water inflow using a pump, and the
water flows from the upper tank into the lower one. The task is to estimate a
model that can predict the water level of the lower tank given a pump input
sequence. The experiment incorporates an overflow of the tank, which
introduces a hard saturation function.

Figure 2: Picture of the _Cascaded Tank System_ (CTS) (Schoukens et al.,
2016).

Model Configuration. We chose the fixed-step Runge-Kutta 4 ODE solver. The
network ${f_{{}_{\mathrm{NN}}}}(\cdot,\cdot)$ is chosen to consist of linear
layers with matrices $A\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{x}}}$,
$B\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{u}}}$ and two residual hidden
layers with a Leaky ReLU activation function $\sigma(\cdot)$ and $64$ hidden
units for each layer. Weight matrices are denoted as $W_{(\cdot)}$ and bias
terms as $b_{(\cdot)}$ with appropriate dimensions. Bias terms for the linear
layer are disabled. An output network ${g_{{}_{\mathrm{NN}}}}(\cdot,\cdot)$ is
also used with the same structure as the state network
${f_{{}_{\mathrm{NN}}}}(\cdot,\cdot)$. The overall NN model can be written as:

$\displaystyle{\frac{dx(t)}{dt}}$
$\displaystyle={f_{{}_{\mathrm{NN}}}}(x(t),u(t))=Ax(t)+Bu(t)$ (17a)
$\displaystyle+W_{F1}\sigma\left(W_{F2}\sigma\left(W_{F3}\begin{bmatrix}x(t)\\\
u(t)\end{bmatrix}+b_{F3}\right)+b_{F2}\right)$ $\displaystyle y(t)$
$\displaystyle={g_{{}_{\mathrm{NN}}}}(x(t),u(t))=Cx(t)+Du(t)$ (17b)
$\displaystyle+W_{G1}\sigma\left(W_{G2}\sigma\left(W_{G3}\begin{bmatrix}x(t)\\\
u(t)\end{bmatrix}+b_{G3}\right)+b_{G2}\right)$

White box modeling of the CTS would lead to $2$ states. Nevertheless, we set
the number of states $n_{\mathrm{x}}=4$. It is observed that this results in
better learning behavior. This is motivated by the effect of state
augmentation (Dupont et al., 2019), corresponding to the fact that increasing
the state dimension allows to describe the system behavior with less complex
nonlinearities involved (also expressed by the immersion concept in nonlinear
system theory). The initial hidden states are obtained using a _Deep Encoder
Network_ (Beintema et al., 2023) ${e_{{}_{\mathrm{NN}}}}(\cdot,\cdot)$ with
$n_{\mathrm{a}}=n_{\mathrm{b}}=5$. Weight matrices are initialized with a
small random number chosen from a uniform distribution
$\mathcal{U}[-0.01,0.01]$, while bias terms are initialized to zero. A barrier
function $L_{\mathrm{N}}$ is applied to the linear layer $A$ of the state
network to ensure negative definiteness (Weigand et al., 2021). This barrier
minimizes drifts over long simulation horizons and is estimated using the
differentiable Sylvester Criterion. Furthermore, we apply a _Differential
Algebraic Equations_ (DAE) network

$\displaystyle d_{{}_{\mathrm{NN}}}(x(t),u(t))=$ (18) $\displaystyle
W_{D1}\sigma\left(W_{D2}\sigma\left(W_{D3}\begin{bmatrix}x(t)\\\
u(t)\end{bmatrix}+b_{D3}\right)+b_{D2}\right)+b_{D1}$

to account for the hard state saturation, which is trained using an additional
penalty function $L_{\mathrm{D}}$. It does not affect the forward model
evaluation. The optimization problem is given by

$\displaystyle\min_{\theta}$ $\displaystyle L_{\mathrm{N}}+L_{\mathrm{W}}+$
$\displaystyle\sum_{k=\max(n_{\mathrm{a}},n_{\mathrm{b}})}^{K-J}\left(\sum_{j=0}^{J-1}\|y_{k+j}-\hat{y}_{k+j|k}\|_{2}^{2}+L_{\mathrm{D},j}\right),$
s.t. $\displaystyle\hat{y}_{k+j|k}=g_{{}_{\mathrm{NN}}}(x_{k+j|k}),$
$\displaystyle x_{k+j+1|k}$
$\displaystyle=\text{ODE\\_Solve}\left(x_{k+j|k},u_{k+j},\frac{1}{\tau}f_{{}_{\mathrm{NN}}}(\cdot,\cdot),T_{\mathrm{s}}\right),$
$\displaystyle
x_{k|k}=e_{{}_{\mathrm{NN}}}(u_{k-1},...,u_{k-n_{b}},y_{k-1},...,y_{k-n_{a}}),$
$\displaystyle
L_{D,j}=\lambda_{\mathrm{D}}\,\left(d_{{}_{\mathrm{NN}}}(x_{j|j},u_{j|j})\right)^{2},$
$\displaystyle L_{\mathrm{N}}=\begin{cases}0&\textrm{if}\ A\prec 0\\\
\lambda_{\mathrm{N}}&\textrm{otherwise}.\end{cases},$ $\displaystyle
L_{\mathrm{W}}=\lambda_{\mathrm{W}}\|\theta\|_{2}^{2}$

with $\lambda_{\mathrm{N}}=10^{12}$ and $\lambda_{\mathrm{D}}=10^{3}$. The bar
notation $x_{k+j|k}$ reads as ”The simulated state at $x_{k+j}$ starting at
$k$ with initial state $x_{k|k}$”. The notation $\|\cdot\|_{2}^{2}$ stands for
the squared Euclidean norm. The sampling time is $T_{\mathrm{s}}=4$ s.

The implementation is written in Python, using PyTorch. Furthermore, the ADAM
optimizer is applied with unmodified configuration except for the learning
rate, which is set to $0.003$ for the first $1{,}000$ steps, and to $0.0009$
until step $3{,}000$, and to $0.00027$ for all subsequent iterations. The
maximum number of optimization steps is set to $20,000$. Regularization is
obtained using weight decay of $\lambda_{\mathrm{W}}=10^{-8}$ and early
stopping when the best validation error does not improve for $2{,}000$
iterations. We do not apply _Batch-Norm_ or _Dropout_. Input and output data
are z-score normalized. As no explicit validation data set is given for CTS,
we apply the first $512$ time steps of the test data for early stopping. Test
data is not accessed for any other reason. Training is performed in $64$ mini-
batches with a sequence length of $J=128$ steps using _Truncated Simulation
Error Minimization_ (TSEM) (Forgione and Piga, 2021).

Effect of Normalization. We analyze the effect of the normalization factor
$\tau$ given the CTS benchmark, a fixed NN configuration in simulation mode
(free-run simulation/simulation error), and a fixed training pipeline
throughout all experiments. Performance is measured in terms of the _Root Mean
Squared Error_ (RMSE) of the simulation error w.r.t. the test data:

$e_{\mathrm{RMSE}}=\sqrt{\frac{1}{K-\max(n_{\mathrm{a}},n_{\mathrm{b}})}\sum_{k=\max(n_{\mathrm{a}},n_{\mathrm{b}})}^{K-1}\|y_{k}-\hat{y}_{k}\|_{2}^{2}}.$

Fig. 3 displays a grid search for different scalar values of
$T_{\mathrm{s}}/\tau$ fixed prior to the model training on a log scale
($T_{\mathrm{s}}=4\,\mathrm{s}$ for CTS). Each experiment is repeated $20$
times to account for and observe the effect of randomness in the initial
weights. We observe in Fig. 3 that both the performance and variance of the
results get worse for large and small values of $\tau$. Furthermore, Fig. 3
indicates that there exists a desirable optimum. The optimum is different from
the unscaled function estimation with $\tau=1$, which corresponds to
$T_{\mathrm{s}}/\tau=4.0$ for CTS in Fig. 3. The median RMSE for
$T_{\mathrm{s}}/\tau=4.0$ is $0.75(V)$ and the mean RMSE is $1.98(V)$,
considerably larger than the normalized result. The model outputs obtained
with the proposed estimators for $\tau$ are displayed in Figure 4.

0.00010.0020.030.10.52.34.09.540$0$$0.2$$0.4$$0.6$$0.8$$1$Normalized time
$T_{s}/\tau$ (unitless)Error $e_{RMSE}$ (V) Figure 3: $200$ independent
experiments, each with the same configuration except for the random weight
initialization and a fixed scalar normalization factor. We tested $10$
different normalization factors repeated $20$ times each. The box plot
displays the median, lower quartile, upper quartile, minimum and maximum
values. Note that experiments with a normalization $T_{\mathrm{s}}/\tau=40$
sometimes lead to unstable results, with RMSE $>10^{9}$.

Using the cross-validation method to estimate the normalization, we observe in
Fig. 3 good results between $T_{\mathrm{s}}/\tau=0.031$ and
$T_{\mathrm{s}}/\tau=2.276$, with a best RMSE at $T_{\mathrm{s}}/\tau=0.543$.
For the trainable normalization method, we define a normalization vector
$\tau$ as a trainable parameter and implement (8). It is not optimized with
hyperparameter tuning but jointly trained with all NN weights, using the same
learning rate. The normalization is initialized to $T_{\mathrm{s}}/\tau=0.1$.
After training $20$ models, we obtain an average normalized time constant
$T_{\mathrm{s}}/\tau=0.072\pm 0.048\,(\textrm{mean}\pm\textrm{std})$. Using
the BLA and (12), we estimate a normalization factor
$T_{\mathrm{s}}/\tau=0.054$, which matches the results in Fig. 3.

$0$$500$$1{,}000$$1{,}500$$2{,}000$$2{,}500$$3{,}000$$3{,}500$$4{,}000$$4$$6$$8$$10$Time
(s)Waterlevel (V)TargetTrainedCross-ValBLA
---
Figure 4: Simulation results on the CTS benchmark (Schoukens et al., 2016).
Results of the best models obtained with the trained normalization factor, the
cross-validation method, and the _Best Linear Approximation_ are displayed.

Comparison to Literature. Since several black-box identification methods have
been applied to the CTS in the literature, we can compare our method to the
black-box approaches with the best performance (in terms of RMSE on the test
data). This comparison is provided in Table 1. The following approaches are
compared: (Relan et al., 2017) estimates a BLA and develops an unstructured
_Nonlinear State Space Model_ (NLSS) with different initialization schemes. A
nonparametric Volterra series model is estimated in (Birpoutsoukis et al.,
2018). A nonlinear state-space model based on Gaussian processes is applied in
(Svensson and Schön, 2017). As a direct comparison, continuous-time NN on this
modeling problem has been applied in (Beintema et al., 2023; Forgione and
Piga, 2021; Mavkov et al., 2020; Weigand et al., 2021). These works emphasize
initial state estimation, fitting criteria, and model stability differently.
It can be observed that the proposed approach outperforms all other black-box
identification approaches.

Table 1: Results for the Cascaded Tank Benchmark. | test data
---|---
|  $e_{RMSE}$ (V)
Best Linear Approximation | 0.75
Truncated Volterra Model | 0.54
State-space with GP-inspired Prior | 0.45
Integrated Neural Networks | 0.41
Soft-constrained Integration Method | 0.40
Stable Runge-Kutta Neural Network | 0.39
Nonlinear State Space Model | 0.34
Truncated Simulation Error Minimization | 0.33
Deep Subspace Encoder | 0.22
ours (SDN, Trained Parameter, best model) | 0.2151
ours (SDN, Trained Parameter, mean $\pm$ std) | 0.2977 $\pm$ 0.1259
ours (SDN, Cross-Validation, best model) | 0.2054
ours (SDN, Cross-Validation, mean $\pm$ std) | 0.2777 $\pm$ 0.052
ours (SDN, BLA, best model) | 0.2253
ours (SDN, BLA, mean $\pm$ std) | 0.2633 $\pm$ 0.0284

##  6 Conclusion

We have shown the importance of proper state normalization when considering
continuous-time modeling with state-space NN’s. This is handled by introducing
a normalization constant in front of the state derivative network. We have
provided a state domain, a state derivative domain, and a time domain
interpretation of this concept, and showed the beneficial effects of such a
normalization. To estimate the appropriate normalization constant, three
approaches based on a trainable parameter, cross-validation, and BLA have been
proposed to ensure the practical applicability of the normalization. Based on
simulation studies, we have shown that the proposed methodologies enable to
improve model estimation with NN-based state-space modeling methods.

## References

  * Ba et al. (2016) Ba, J.L., Kiros, J.R., and Hinton, G.E. (2016).  Layer normalization.  URL http://arxiv.org/pdf/1607.06450v1.
  * Beintema et al. (2023) Beintema, G.I., Schoukens, M., and Tóth, R. (2023).  Continuous-time identification of dynamic state-space models by deep subspace encoding.  In _The Eleventh International Conference on Learning Representations_.
  * Birpoutsoukis et al. (2018) Birpoutsoukis, G., Csurcsia, P.Z., and Schoukens, J. (2018).  Efficient multidimensional regularization for volterra series estimation.  _Mechanical Systems and Signal Processing_ , 104, 896–914.  10.1016/j.ymssp.2017.10.007.
  * Chen et al. (2018) Chen, R.T.Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. (2018).  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31.  URL http://arxiv.org/pdf/1806.07366v5.
  * Dupont et al. (2019) Dupont, E., Doucet, A., and Teh, Y.W. (2019).  Augmented neural ODEs.  _Advances in neural information processing systems_ , 32.  URL http://arxiv.org/pdf/1904.01681v3.
  * Forgione and Piga (2021) Forgione, M. and Piga, D. (2021).  Continuous-time system identification with neural networks: Model structures and fitting criteria.  _European Journal of Control_ , 59, 69–81.  10.1016/j.ejcon.2021.01.008.
  * Glorot and Bengio (2010) Glorot, X. and Bengio, Y. (2010).  Understanding the difficulty of training deep feedforward neural networks.  _International Conference on Artificial Intelligence and Statistics (AISTATS)_.
  * Hasani (2020) Hasani, R. (2020).  _Interpretable Recurrent Neural Networks in Continuous-time Control Environments_.  Ph.d. thesis, TU Wien.
  * Ioffe and Szegedy (2015) Ioffe, S. and Szegedy, C. (2015).  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , 448–456. PMLR.
  * Mavkov et al. (2020) Mavkov, B., Forgione, M., and Piga, D. (2020).  Integrated neural networks for nonlinear continuous-time system identification.  _IEEE Control Systems Letters_ , 1.  10.1109/LCSYS.2020.2994806.
  * Pintelon and Schoukens (2012) Pintelon, R. and Schoukens, J. (2012).  _System identification: a frequency domain approach_.  John Wiley & Sons.
  * Relan et al. (2017) Relan, R., Tiels, K., Marconato, A., and Schoukens, J. (2017).  An unstructured flexible nonlinear model for the cascaded water-tanks benchmark.  _IFAC-PapersOnLine_ , 50(1), 452–457.  10.1016/j.ifacol.2017.08.074.
  * Schoukens and Ljung (2019) Schoukens, J. and Ljung, L. (2019).  Nonlinear system identification: A user-oriented road map.  _IEEE Control Systems Magazine_ , 39(6), 28–99.
  * Schoukens et al. (2016) Schoukens, M., Mattsson, P., Wigren, T., and Noel, J.P. (2016).  Cascaded tanks benchmark combining soft and hard nonlinearities.  _Workshop on Nonlinear System Identification Benchmarks, Brussels, Belgium_.
  * Schoukens and Tóth (2020) Schoukens, M. and Tóth, R. (2020).  On the initialization of nonlinear LFR model identification with the best linear approximation.  _IFAC-PapersOnLine_ , 53(2), 310–315.  10.1016/j.ifacol.2020.12.142.
  * Svensson and Schön (2017) Svensson, A. and Schön, T.B. (2017).  A flexible state–space model for learning nonlinear dynamical systems.  _Automatica_ , 80, 189–199.  10.1016/j.automatica.2017.02.030.
  * Wang and Lin (1998) Wang, Y.J. and Lin, C.T. (1998).  Runge-kutta neural network for identification of dynamical systems in high accuracy.  _IEEE Transactions on Neural Networks_ , 9(2), 294–307.
  * Weigand et al. (2021) Weigand, J., Deflorian, M., and Ruskowski, M. (2021).  Input-to-state stability for system identification with continuous-time runge–kutta neural networks.  _International Journal of Control_ , 1–17.  10.1080/00207179.2021.1978555.

  *[DT]: Discrete-Time
  *[CT]: Continuous-Time
  *[NODE]: Neural Ordinary Differential Equations
  *[RKNN]: Runge-Kutta Neural Networks
  *[LTC]: Liquid Time Constant Neural Networks
  *[SDN]: State Derivative Normalization
  *[SDN’s]: State Derivative Normalization
  *[ODE]: Ordinary Differential Equations
  *[NN’s]: Neural Network
  *[CTS]: Cascaded Tank System
  *[NN]: Neural Network
  *[BLA]: Best Linear Approximation
  *[PEM]: Prediction Error Minimization
  *[DFT]: Discrete Fourier Transform
  *[DAE]: Differential Algebraic Equations
  *[TSEM]: Truncated Simulation Error Minimization
  *[RMSE]: Root Mean Squared Error
  *[NLSS]: Nonlinear State Space Model


# Determining Neutrino mass hierarchy in an extended Left-Right Model

Prativa Pritimita∗ Indian Institute of Technology Bombay, Mumbai, India
Nitali Dash  School of Physics, Regional Institute of Education (NCERT),
Bhubaneswar, India  Sudhanwa Patra  Indian Institute of Technology Bhilai,
Chhatisgarh, India  Urjit A. Yajnik  Indian Institute of Technology Bombay,
Mumbai, India

###### Abstract

We derive the lower bound on absolute scale of lightest neutrino mass for
normal hierarchy and inverted hierarchy pattern of light neutrinos by studying
the new physics contributions to charged lepton flavour violating decays and
neutrinoless double beta decay in the framework of a TeV scale left-right
symmetric model. Neutrino mass is generated in the model via type-II seesaw
dominance with the help of a heavy sterile neutrino. This scenario allows
large light-heavy neutrino mixing and the mixing facilitates new channels for
neutrinoless double beta decay and lepton flavour violating decays. We express
all the model parameters in terms of oscillation parameters in order to
constrain light neutrino mass scale and extract information on neutrino mass
hierarchy.

##  I Introduction

Left-Right Symmetric Model (LRSM) Mohapatra:1974gc ; Pati:1974yy ;
Senjanovic:1975rk ; Senjanovic:1978ev  was first proposed to explain parity
restoration at high scale. The theory can also explain the origin of neutrino
mass through seesaw mechanism since it has a right-handed neutrino. Many forms
of LRSM are possible. In generic LRSM, the symmetry breakings happen with the
help of scalar triplets and bidoublet which leads to the generation of
neutrino mass through type-I plus type-II seesaw mechanism.

$M_{\nu}\ =\ -M_{D}M_{R}^{-1}M_{D}^{T}+M_{L}\ \equiv\ M_{\nu}^{\rm
I}+M_{\nu}^{\rm II}$

Here $M_{D}$ is the Dirac neutrino mass induced by the bidoublet vacuum
expectation value (VEV), while $M_{R}$ and $M_{L}$ are the Majorana masses of
the right and left-handed neutrinos respectively induced by the triplet VEVs.
For phenomenological purposes, usually one of the contributions is assumed to
be dominant with observable ramifications for different experiments. For
type-I seesaw dominance, it is assumed $M_{L}\to 0$ and the light neutrino
mass depends on the Dirac mass matrix $M_{D}$. However, in this case all the
light-heavy neutrino mixing effects are suppressed for TeV-scale parity
restoration. Type-II seesaw dominance can be realized either by assuming
$M_{D}\to 0$ or with very high scale of parity restoration. Here the light and
heavy neutrino mass matrices are directly proportional to each other and lead
to a more predictive scenario. But in this case also, the light-heavy neutrino
mixing effects on lepton number or lepton flavor violating observables are
lost.

We have extended the generic LRSM only in the fermion sector by adding a
sterile neutrino $S_{L}$ per generation which helps in exactly cancelling the
type-I seesaw term and achieving natural type-II seesaw dominance
Pritimita:2016fgr ; Dash:2021pbx . This scenario also allows large light-heavy
neutrino mixing which has important and non-trivial phenomenological
consequences. This mixing in the neutrino sector facilitates new channels for
lepton number and lepton flavour violating decays. We have expressed these new
physics contributions to cLFV and $0\nu\beta\beta$ decay in terms of the
observed light neutrino oscillation parameters and lightest neutrino mass. As
a result of this, we can derive constraints on the lightest neutrino mass from
the non-observation of these rare processes. Even though many new channels are
possible for these decays in the model, in order to highlight the
contributions of right-handed heavy neutrino and sterile neutrino to LFV
decays and $0\nu\beta\beta$ decay, we have focussed only on the diagrams
mediated by them and ignored other possible channels.

##  II The Model Framework

The left-right symmetric theory is based on the gauge group
$\mathcal{G}_{LR}\equiv SU(3)_{c}\times SU(2)_{L}\times SU(2)_{R}\times
U(1)_{B-L}$. The particle content of our model can be found in Table.1. The
fermion sector contains all the usual quarks and leptons plus one extra
sterile neutrino $S_{L}$. The scalar sector contains the doublets $H_{L}$,
$H_{R}$, triplets $\Delta_{L},\Delta_{R}$ and bidoublet $\Phi$. In the model,
both doublets and triplets are used for the left-right symmetry breaking.

Table 1: Particle content of the model | Fields | $SU(3)_{c}$ | $SU(2)_{L}$ | $SU(2)_{R}$ | $B-L$
---|---|---|---|---|---
| $q_{L}$ | 3 | 2 | 1 | 1/3
| $q_{R}$ | 3 | 1 | 2 | 1/3
| $\ell_{L}$ | 1 | 2 | 1 | $-1$
| $\ell_{R}$ | 1 | 1 | 2 | $-1$
| $S$ | 1 | 1 | 1 | 0
| $\Phi$ | 1 | 2 | 2 | 0
| $H_{L}$ | 1 | 2 | 1 | $-1$
| $H_{R}$ | 1 | 1 | 2 | $-1$
| $\Delta_{L}$ | 1 | 3 | 1 | 2
| $\Delta_{R}$ | 1 | 1 | 3 | 2

With these fermions and scalars, the interaction lagrangian for leptons can be
written as,

$\displaystyle-\mathcal{L}_{Yuk}=\overline{\ell_{L}}\left[Y_{3}\Phi+Y_{4}\widetilde{\Phi}\right]\ell_{R}+f\left[\overline{(\ell_{L})^{c}}\ell_{L}\Delta_{L}+\overline{(\ell_{R})^{c}}\ell_{R}\Delta_{R}\right]$
$\displaystyle+F\overline{(\ell_{R})}H_{R}S^{c}_{L}+F^{\prime}\overline{(\ell_{L})}H_{L}S_{L}+\mu_{S}\overline{S^{c}_{L}}S_{L}+\mbox{h.c.}$
(1)

After spontaneous symmetry breaking, the neutral lepton mass matrix in the
basis $\left(\nu_{L},N^{c}_{R},S_{L}\right)$ can be written as,

$\mathcal{M}_{\nu}=\left(\begin{array}[]{ccc}M_{L}&M_{D}&0\\\
M^{T}_{D}&M_{R}&M^{T}\\\ 0&M&0\end{array}\right)\,,$ (2)

where $M_{D}$ is the Dirac neutrino mass matrix connecting $\nu_{L}$ and
$N_{R}$, $M$ is the mixing matrix in the $N_{R}-S_{L}$ sector, $M_{L}$
$(M_{R})$ is the Majorana mass matrix for left-handed (right-handed) active
neutrinos generated dynamically by non-zero VEV of scalar triplet $\Delta_{L}$
$(\Delta_{R})$. Assuming the mass hierarchy, $M_{R}\gg M>M_{D}\gg M_{L}$, the
light and heavy neutrino mass formulae can be written as,

$\displaystyle m_{\nu}=M_{L}$ (3) $\displaystyle m_{N}=M_{R}$ (4)
$\displaystyle m_{S}\simeq MM^{-1}_{R}M^{T}$ (5)

For detail derivation of neutrino masses and mixing in the model, please refer
Dash:2021pbx .

##  III LFV in LRSM

In our model, LFV decays can be mediated by heavy right-handed neutrinos
$N_{R}$, sterile neutrinos $S_{L}$, charged scalar triplets
$\Delta^{\pm\pm}_{L,R},\Delta^{\pm}_{L,R}$ due to purely left-handed currents
involving $W_{L}$ and due to purely right-handed currents involving $W_{R}$.
However we focus only on those contributions which involve large active-
sterile neutrino mixing, i.e. due to the neutrinos $N_{R}$ and $S_{L}$ in
order to constrain light neutrino masses from LFV decays. We ignore other
contributions by imposing the limiting conditions on masses of $W_{R}$ boson
and scalar trilpets $\Delta_{L,R}$ as, $M_{W_{R}}\gg M_{W_{L}}$ and
$M_{\Delta_{L,R}},M_{H_{1}}\gg M_{N,S}$.

Figure 1: Feynman diagram for lepton flavor violating process $\mu\to e\gamma$
due to the exchange of mass eigenstates of heavy neutrinos $N_{i}$ and
$S_{i}$.

The effective Lagrangian relevant for the lepton flavor violating process
$\mu\to e\gamma$ in our present work can be expressed as,

$\displaystyle\mathcal{L}_{\mu\rightarrow
e}=-\frac{eg^{2}}{4(4\pi)^{2}m_{W_{L}}^{2}}m_{\mu}\overline{e}\sigma_{\mu\nu}(\mathcal{G}^{\gamma}_{L}P_{L}+\mathcal{G}^{\gamma}_{R}P_{R})\mu
F^{\mu\nu}$
$\displaystyle-\frac{\alpha_{W}^{2}}{2m_{W_{L}}^{2}}\sum_{q}\left\\{\overline{e}\gamma_{\mu}\left[\mathcal{W}^{q}_{L}P_{L}+\mathcal{W}^{q}_{R}P_{R}\right]\mu\;\overline{q}\gamma^{\mu}q\right\\}+{\rm
h.c.}$ (6)

where $\sigma_{\mu\nu}\equiv\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$ and the
form factors are $\mathcal{G}^{\gamma}_{L,R}$ and $\mathcal{W}^{u,d}_{L,R}$.
The relevant contributions to $\mu\to e\gamma$ is given by

$\displaystyle i\mathcal{M}(\mu\rightarrow e\gamma)=\frac{e\alpha_{W}}{8\pi
m_{W_{L}}^{2}}\epsilon_{\gamma}^{\mu}\overline{e}\mu\left[\left(q^{2}\gamma_{\mu}-q_{\mu}\not{q}\right)\right.$
$\displaystyle\left.\left(\mathcal{F}^{\gamma}_{L}P_{L}+\mathcal{F}^{\gamma}_{R}P_{R}\right)-im_{\mu}\sigma_{\mu\nu}q^{\nu}\left(\mathcal{G}^{\gamma}_{L}P_{L}+\mathcal{G}^{\gamma}_{R}P_{R}\right)\right]$
(7)

where $\mathcal{F}^{\gamma}_{L,R}$ and $\mathcal{G}^{\gamma}_{L,R}$ are the
anapole and dipole form factors.

Figure 2: Branching ratio of the lepton flavor violating process, $\mu\to
e\gamma$ as a function of lightest neutrino mass $m_{1}$ for NH and $m_{3}$
for IH. The blue (NH) and red (IH) coloured regions display the model
predictions on $\mu\rightarrow e\gamma$ due to the exchange of heavy right
handed neutrino ($N_{R}$) (upper plot) and heavy sterile neutrino ($S_{L}$)
(lower plot).

The analytic expression for the branching ratio for the lepton flavor
violating process $\mu\to e\gamma$ shown in Feynman diagram in Fig.1 due to
mediation of heavy neutrinos ($N$ and $S$) is given by

$\text{Br}_{\mu\to
e\gamma}=\frac{\alpha_{W}^{3}s_{W}^{2}}{256\pi^{2}}\frac{m_{\mu}^{4}}{M_{W_{L}}^{4}}\frac{m_{\mu}}{\Gamma_{\mu}}|\mathcal{G}_{\gamma}^{\mu
e}|^{2}\,,$ (8)

where $\alpha_{W}=1/29.0$ is weak fine struture constant, $m_{\mu}=105$ MeV is
the muon mass, $M_{W_{L}}$ is the SM $W$-boson mass,
$s_{W}\equiv\sin\theta_{W}$ and $\Gamma_{\mu}=2.996\times 10^{-19}$ GeV is the
total decay width of the muon. The important parameter,$G_{\gamma}^{\mu e}$,
for deriving constraints on light neutrino masses, is of the following form,

$\displaystyle G_{\gamma}^{\mu
e}=\bigg{|}\sum_{i=1}^{3}\bigg{\\{}{\mbox{V}^{\nu N}_{\mu
i}}^{*}{\mbox{V}^{\nu
N}_{ei}}\mathcal{G}_{\gamma}\left(x_{N_{i}}\right)+{\mbox{V}^{\nu S}_{\mu
i}}^{*}{\mbox{V}^{\nu
S}_{ei}}\mathcal{G}_{\gamma}\left(x_{S_{i}}\right)\bigg{\\}}\bigg{|}^{2}$ (9)

where $x_{N_{i}}=m_{N_{i}}^{2}/M_{W_{L}}^{2}$,
$x_{S_{i}}=m_{S_{i}}^{2}/M_{W_{L}}^{2}$ and the form of loop function is given
by

$\displaystyle G_{\gamma}(x)$ $\displaystyle=$
$\displaystyle-\frac{x(2x^{2}+5x-1)}{4(1-x)^{3}}-\frac{3x^{3}}{2(1-x)^{4}}\ln
x\,$ (10)

The form of loop function $\mathcal{G}_{\gamma}(x)$ and its dependance with
change of lightest neutrino mass $m_{1}$ (for NH) and $m_{3}$ (for IH) can be
found in the appendix of the main manuscript. The other parameters
$\mbox{V}^{\nu N}$ and $\mbox{V}^{\nu S}$ are mixing matrices representing
mixing of light neutrinos with $N_{R}$ and $S_{L}$, respectively.

Fig.2 shows the variation of $\text{Br}_{\mu\to e\gamma}$ as a function of
lightest neutrino mass with contributions coming from purely $N_{R}$
(presented in the upper plot) and $S_{L}$ (presented in the lower plot). The
blue color (NH) and red color (IH) regions are model prediction on $\mu\to
e\gamma$ within $3-\sigma$ allowed range of neutrino mixing angles and mass
squared differences. In x-axis, $m_{1}$ for NH ($m_{3}$ for IH) represents
absolute value of light neutrino mass. The current experimental limit and
future sensitivity by MEG ($\text{Br}_{\mu\to e\gamma}\leq 4.2\times
10^{-13}$) MEG:2013oxv , MEG upgrade ($\text{Br}_{\mu\to e\gamma}\leq
5.0\times 10^{-14}$) Baldini:2013ke  and PRISM/PRIME ($\text{Br}_{\mu\to
e\gamma}\leq 1.0\times 10^{-16}$) Kuno:2005mm  on branching ratio for $\mu\to
e\gamma$ are presented in dashed horizontal lines. The vertical shaded regions
are PLANCK bound on lightest neutrino mass with 95$\%$ C.L. Tian:2020tur . the
experimental bounds are satisfied by the model predicted branching ratio and
from these plots one can derive bound on lightest neutrino mass as follows.

  * •

Saturating the MEG bound ($\text{Br}_{\mu\to e\gamma}\leq 4.2\times
10^{-13}$), one can derive the limit on lightest neutrino mass in the range of
meV scale in IH case while most of the parameters are ruled out for NH case.

  * •

Saturating the MEG Upgrade bound ($<5.0\times 10^{-14}\,$) one can find
lightest neutrino mass less than few meV for NH case and $4$ meV for IH case.
However, most of the paremeter spaces satisfying PRISM/PRIME ($<1.0\times
10^{-16}\,$) bound are lying in the quasi-degenerate (QD) pattern of light
neutrino masses which already ruled out by PLANCK data.

##  IV LNV in LRSM

We have also discussed how the light-heavy neutrino mixing leads to sizable
new contributions to neutrinoless double beta decay in the model. We emphasize
on left-handed current effects due to the exchange of heavy neutrinos $N_{R}$
and $S_{L}$.

Figure 3: Feynman diagrams for $0\,\nu\,\beta\beta$ transition due to purely
left-handed currents with the exchange of light neutrino $\nu_{i}$ and heavy
neutrinos $N_{j}$ and $S_{k}$ in flavour basis.

The inverse half-life of $0\nu\beta\beta$ transition for a given isotope in
terms of Effective Majorana mass can be written as,

$\displaystyle\left[T_{1/2}^{0\nu}\right]^{-1}$ $\displaystyle=$
$\displaystyle G^{0\nu}_{01}\bigg{|}\frac{{\cal
M}^{0\nu}_{\nu}}{m_{e}}\bigg{|}^{2}\bigg{[}\big{|}m^{\nu}_{ee}\big{|}^{2}+\big{|}m^{N}_{ee}\big{|}^{2}+\big{|}m^{S}_{ee}\big{|}^{2}\bigg{]}$
(11) $\displaystyle=$ $\displaystyle G^{0\nu}_{01}\bigg{(}\frac{{\cal
M}^{0\nu}_{\nu}}{m_{e}}\bigg{)}^{2}\cdot|m^{\rm eff}_{\beta\beta}|^{2}\,.$

Here $G^{0\nu}$ is the phase-space factor, $\mathcal{M}^{0\nu}_{\nu}$ is the
Nuclear Matrix Element and $|m^{\rm eff}_{\beta\beta}|^{2}$ is the sum of
contributions from light active neutrinos $\nu_{L}$, heavy right-handed
neutrinos $N_{R}$ and sterile neutrinos $S_{L}$.

Table 2: $G^{0\nu}_{01}$ and NMEs Meroni:2012qf  Isotope |  $G^{0\nu}_{01}$ $[{\rm yrs}^{-1}]$ | ${\cal M}^{0\nu}_{\nu}$ | ${\cal M}^{0\nu}_{N}$
---|---|---|---
76Ge | $5.77\times 10^{-15}$ | 2.58–6.64 | 233–412
136Xe | $3.56\times 10^{-14}$ | 1.57–3.85 | 164–172

Figure 4: Upper plot: New physics contribution to effective Majorana mass as a function of lightest neutrino mass, $m_{1}$ ($m_{3}$) for NH (IH) with the exchange of heavy neutrino $N_{R}$ and sterile neutrino $S_{L}$. Lower plot: Contributions of $N_{R}$ and $S_{L}$ to half-life vs lightest neutrino mass. Table 3: Limits on the half-life of $0\nu\beta\beta$. Experiment | Limit
---|---
HM | $1.9\times 10^{25}\ {\rm yrs}$
GERDA | $2.1\times 10^{25}\ {\rm yrs}$
Combined 76Ge | $3.0\times 10^{25}\ {\rm yrs}$
GERDA Phase II | $5.2\times 10^{25}\ {\rm yrs}$
EXO | $1.6\times 10^{25}\ {\rm yrs}$
KamLAND-Zen | $1.9\times 10^{25}\ {\rm yrs}$
Combined 136Xe | $3.4\times 10^{25}\ {\rm yrs}$
nEXO | $6.6\times 10^{27}\ {\rm yrs}$

Fig.4 shows the effective Majorana mass parameter in the upper plot and half-
life in the lower plot due to the contributions of $N_{R}$ and $S_{L}$ as a
function of lightest neutrino mass $m_{1}$(NH) and $m_{3}$(IH). The horizontal
blue bands stand for the improved limits on $0\nu\beta\beta$ decay on
effective Majorana mass and half-life with the combined results of GERDA
GERDA:2013vls  and KamLAND-Zen KamLAND-Zen:2012mmx  experiments. The vertial
magenta region is disfavoured by Planck-2018 Tian:2020tur  and KATRIN data
KATRIN:2019yun  on sum of light neutrino masses. The contributions of $N_{R}$
and $S_{L}$ to effective Majorana mass parameter and half-life predictions
makes both NH(blue band) and IH(orange band) patterns of light neutrino masses
sensitive to the current experimental bounds. In these plots the largest value
of right-handed neutrino mass is fixed at $M_{N_{R}}=1$ TeV, the Majorana
phases and Dirac CP-phase are varied between $0\to\pi$, and other neutrino
oscillation parameters are taken in their allowed $3\sigma$ range. From these
plots it is observed that, when the contributions of $N_{R}$ and $S_{L}$ are
considered, the allowed values of lightest neutrino mass are found to be in
the range of $10-25$ meV for $m_{1}$(NH) and $25-40$ meV for $m_{3}$(IH). This
is done by saturating the effective Majorana mass with the current GERDA and
KamLAND-Zen bounds. It remains the same for half-life.

##  V Correlation between $m_{ee}$ and sum of light neutrino masses $\sum
m_{i}$

We have also discussed the model predictions on effective Majorana mass by
changing the sum of light neutrino masses $\Sigma m_{i}$. We use the following
limits on sum of light neutrino masses SDSS:2004kqt ; Costanzi:2014tna ;
Palanque-Delabrouille:2014jca  for displaying the allowed region of $m_{ee}$
in Fig.5.

$\displaystyle m_{\Sigma}<\mbox{84\,meV}\quad\quad~{}(1\sigma~{}\mbox{C.L.})$
$\displaystyle m_{\Sigma}<\mbox{146\,meV}\quad\quad(2\sigma~{}\mbox{C.L.})$
$\displaystyle m_{\Sigma}<\mbox{208\,meV}\quad\quad(3\sigma~{}\mbox{C.L.})$
(12)

In both the plots of Fig.5 the horizontal band represents the experimental
bounds on effective Majorana mass $m_{ee}$ by EXO+KamLAND-Zen while the
vertical dashed line shows the bound on sum of light neutrino masses from
cosmology. The plot in the upper-panel shows the variation of effective
Majorana mass parameter due to standard mechanism with sum of light neutrino
masses. It shows that both NH (green band) and IH (red band) patterns of light
neutrino masses are not sensitive to the current experimetal bounds on
$m_{ee}$ and also disfavoured by cosmology. The plot in the lower panel shows
the variation of effective Majorana mass parameter due to new physics
contributions (contributions of $N_{R}$ and $S_{L}$) with sum of light
neutrino masses. It shows that NH (purple band) pattern of light neutrino
masses saturates the experimental bound on $m_{ee}$ as well as the cosmology
bound on sum of light neutrino masses whereas IH (blue band) pattern is
disfavoured by cosmology. Thus the model with type-II seesaw dominance gives
an important result on the hierarchy of light neutrino masses when the
contributions of heavy and sterile neutrinos are considered.

Figure 5: Allowed region of effective Majorana mass parameter ($|m_{ee}|$) as
a function of sum of light neutrino masses ($\Sigma m_{i}$) for standard
mechanism (upper plot) and new contributions from $N_{R}$ and $S_{L}$ mediated
diagrams (lower plot).

##  VI Conclusion

We considered the mechanism of natural type-II seesaw dominance that allows
large light-heavy neutrino mixing in a TeV scale left-right symmetric model.
This becomes possible with the addition of a sterile neutrino to the particle
content. The model also generates large light-heavy neutrino mixing that gives
new contributions to LFV process of $\mu\to e\gamma$ and neutrinoless double
beta decay. As a result of these new contributions, the new branching ratios
can be accessible to present as well as planned experiments. All the physical
masses and mixing of neutral leptons are expressed in terms of light neutrino
masses and PMNS mixing matrix. Thus the new physics contributions to LFV and
LNV processes arising from heavy and sterile neutrinos depend upon light
neutrino mass. As a result, the bound on absolute scale of light neutrino mass
can be derived by saturating the experimental limits on these decay processes,
which we find to be in the meV range. We also find that the new contributions
to $0\nu\beta\beta$ decay can saturate the experimental bound for
$m_{1}>0.001$eV for NH pattern of light neutrino masses.

## References

  * (1) R. N. Mohapatra and J. C. Pati, Phys. Rev. D 11 (1975), 2558
  * (2) J. C. Pati and A. Salam, Phys. Rev. D 10 (1974), 275-289
  * (3) G. Senjanovic and R. N. Mohapatra, Phys. Rev. D 12 (1975), 1502
  * (4) G. Senjanovic, Nucl. Phys. B 153 (1979), 334-364
  * (5) P. Pritimita, N. Dash and S. Patra, JHEP 10 (2016), 147
  * (6) N. Dash, S. Patra, P. Pritimita and U. A. Yajnik,
  * (7) J. Adam et al. [MEG], Phys. Rev. Lett. 110 (2013), 201801
  * (8) A. M. Baldini, F. Cei, C. Cerri, S. Dussoni, L. Galli, M. Grassi, D. Nicolo, F. Raffaelli, F. Sergiampietri and G. Signorelli, et al.
  * (9) Y. Kuno, Nucl. Phys. B Proc. Suppl. 149(2005), 376-378
  * (10) S. X. Tian, Phys. Rev. D 102 (2020) no.6, 063509
  * (11) A. Meroni, S. T. Petcov and F. Simkovic, JHEP 02 (2013), 025
  * (12) M. Agostini et al. [GERDA], Phys. Rev. Lett. 111 (2013) no.12, 122503
  * (13) A. Gando et al. [KamLAND-Zen], Phys. Rev. Lett. 110 (2013) no.6, 062502
  * (14) M. Aker et al. [KATRIN], Phys. Rev. Lett. 123 (2019) no.22, 221802
  * (15) U. Seljak et al. [SDSS], Phys. Rev. D 71 (2005), 103515
  * (16) M. Costanzi, B. Sartoris, M. Viel and S. Borgani, JCAP 10 (2014), 081
  * (17) N. Palanque-Delabrouille, C. Yèche, J. Lesgourgues, G. Rossi, A. Borde, M. Viel, E. Aubourg, D. Kirkby, J. M. LeGoff and J. Rich, et al. JCAP 02 (2015), 045