# Development Of Automated Cardiac Arrhythmia Detection Methods Using Single
Channel ECG Signal

Arpita Paul
Department of Electronics and Telecommunication
Engineering, IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Avik Kumar Das
Department of Electronics and Telecommunication
Engineering, IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Manas Rakshit
Department of Electronics and Communication Engineering
NIT, Agartala, Tripura, India
<EMAIL_ADDRESS>
Ankita Ray Chowdhury
Department of Electronics and Telecommunication
Engineering, IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Susmita Saha
Department of Electronics and Telecommunication
Engineering, IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Hrishin Roy
Department of Electronics and Telecommunication
Engineering, IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Sajal Sarkar
Department of Electronics Science
Seth Anandram Jaipuria College, Kolkata, India
<EMAIL_ADDRESS>
Dongiri Prasanth
Department of Mining Engineering
IIEST, Shibpur, West Bengal, India
<EMAIL_ADDRESS>
Eravelli Saicharan
Department of Mining Engineering
IIEST, Shibpur , West Bengal, India
<EMAIL_ADDRESS>

###### Abstract

Arrhythmia, an abnormal cardiac rhythm, is one of the most common types of
cardiac disease. Automatic detection and classification of arrhythmia can be
significant in reducing deaths due to cardiac diseases. This work proposes a
multi-class arrhythmia detection algorithm using single channel
electrocardiogram (ECG) signal. In this work, heart rate variability (HRV)
along with morphological features and wavelet coefficient features are
utilized for detection of 9 classes of arrhythmia. Statistical, entropy and
energy-based features are extracted and applied to machine learning based
random forest classifiers. Data used in both works is taken from 4 broad
databases (CPSC and CPSC extra, PTB-XL, G12EC and Chapman-Shaoxing and Ningbo
Database) made available by Physionet. With HRV and time domain morphological
features, an average accuracy of 85.11%, sensitivity of 85.11%, precision of
85.07% and F1 score of 85.00% is obtained whereas with HRV and wavelet
coefficient features, the performance obtained is 90.91% accuracy, 90.91%
sensitivity, 90.96% precision and 90.87% F1 score. The detailed analysis of
simulation results affirms that the presented scheme effectively detects broad
categories of arrhythmia from single-channel ECG records. In the last part of
the work, the proposed classification schemes are implemented on hardware
using Raspberry Pi for real time ECG signal classification.

_Keywords_ Arrhythmia detection $\cdot$ ECG $\cdot$ Heart rate variability
$\cdot$ Physionet challenge 2021 database $\cdot$ Random forest $\cdot$ SMOTE
$\cdot$ Stationary wavelet transform $\cdot$ Raspberry Pi

##  1 Introduction

The group of heart disorders, commonly called cardiovascular diseases (CVDs),
have become the cause of an increasing number of premature deaths worldwide
Tsao et al. (2022). Since 1990, prevalent cases of CVDs all over the world
have doubled and the number of deaths due to the same has increased from 12.1
million in 1990 to 18.6 million in 2019 Roth et al. (2020). According to the
world health organization (WHO), 32 % of all global deaths in 2019 were from
CVDs b2 (2021). Over three-quarters of these deaths took place in poor
countries with low doctor-to-patient ratios and inadequate medical
infrastructure. Efficient and automated computer-aided diagnosis of CVDs can
significantly reduce the burden on already strained healthcare systems,
leading to timely detection and treatment of patients that results in reducing
the number of deaths due to chronic heart diseases Faust and Ng (2016).

ECG is the most extensively used non-invasive method for the clinical
detection of CVDs Hong et al. (2022). ECG signal is difficult to analyze
visually for long-term application due to its non-stationary nature. That may
lead to missed or late diagnosis of life-threatening heart ailments. Automated
signal processing with a machine learning algorithm can be developed to
process ECG data in real-time for accurate and prompt detection of cardiac
diseases with less human effort and error Bertsimas et al. (2021), Ran et al.
(2022).

Arrhythmia is a heart condition characterized by an irregular heart rate where
the heart beats either too slow or too fast. It occurs due to improper
electrical impulses that coordinate the heartbeats. Rhythms are classified
into different categories as per their origin like sinus rhythm, atrial
rhythm, atrioventricular (AV) node rhythm, ventricle rhythms etc. Each
category of rhythms is further sub-divided into several classes as per
characteristics of rhythms Walraven (2010).

Methods for diagnosing arrhythmia using classification approaches based on a
variety of ECG parameters have been proposed in a number of promising research
García et al. (2016); Asgari et al. (2015a); Rahul et al. (2021); Mjahad et
al. (2017); Acharya et al. (2016); Malik et al. (2021); Elhaj et al. (2016).
Most of the works have classified individual beats into arrhythmia classes.
Some works have focused on single arrhythmia type detection García et al.
(2016),Asgari et al. (2015a) while other works have classified more than one
class of arrhythmia Rahul et al. (2021); Mjahad et al. (2017); Acharya et al.
(2016); Malik et al. (2021); Elhaj et al. (2016). In García et al. (2016),
relative wavelet energy from wavelet decomposition of T-Q segments in the ECG
cycle is used for the detection of atrial fibrillation (AF). In Asgari et al.
(2015a), an AF detection method is proposed using features from Stationary
Wavelet Transform (SWT) decomposition coefficients and support vector machine
as the classifier. The technique described in Rahul et al. (2021), offers
training several classifiers to distinguish between premature atrial
contraction (PAC), premature ventricular contraction (PVC), and normal beats
based on the R-R interval and other statistical features. In Mjahad et al.
(2017), distinguishing between ventricular fibrillation and ventricular
tachycardia is done using time-frequency representations of ECG signal. The
method in Acharya et al. (2016) recognizes and categorizes four groups of ECG
beats using a set of nonlinear features, including Shannon entropy, fuzzy
entropy, approximate entropy, permutation entropy, sample entropy, etc. When
analyzing patient-specific based arrhythmia classification, the method
described in Malik et al. (2021) employs self-organized operational neural
networks to detect and categorize five groups of arrhythmia beats. In Elhaj et
al. (2016), five categories of heartbeat classification have been done using
wavelet transform coefficients and independent component analysis (ICA).

Most of the previous works in the literature have used either only one or very
limited databases and the performance of such methods are tested only in a
small and homogeneous population. Moreover, the classification ability of the
existing approaches is limited to only very few rhythm types. Broad-range
classification of multiple classes has not been attempted much in the
available literature. Majority of the works in arrhythmia classification have
focused on classifying single heartbeat classification only, not arrhythmic
episode detection. In doing so, a significant amount of inter-beat information
is lost which restricts the classification to only a few classes of
arrhythmia. Considering the above-mentioned aspect, the primary aim of this
work is to propose a multi-class rhythm classification scheme to detect
arrhythmia from ECG records in widely collected different unbalanced
databases. The contributions of the paper are as follows:

  * •

Development of robust multi-class ECG arrhythmia classification schemes which
can effectively detect broad categories of arrhythmia like normal sinus rhythm
(NSR), sinus arrhythmia (SA), sinus bradycardia (SB), sinus tachycardia
(STACH), atrial fibrillation (AF), atrial flutter (AFL), premature atrial
contraction (PAC), 1st Degree AV block (1AVB) and premature ventricular
contraction (PVC).

  * •

Utilization of broadly acquired ECG records combined from four standard
databases of Physionet challenge 2021 which makes the model more robust to
different methods of data acquisition and population demography.

  * •

Extraction of robust heart rate variability and time domain features based on
scientific knowledge of characteristics of ECG signal variations caused by
each of the arrhythmia classes. Also, deriving the feature set from
established facts used for actual pathological diagnosis by physicians which
provides the classifier model interpretability.

  * •

Extraction of effective features using SWT-based sub-band signal
decomposition. Deployment of heart rate variability features along with
wavelet coefficient features for efficient classification of multi-class
cardiac arrhythmia.

  * •

Hardware implementation of the ECG classification algorithms using Raspberry
Pi for real time ECG signal analysis.

The rest of the paper is as follows: Section 2 describes the detailed
information on ECG records and databases which are utilized in this work.
Sections 3, 4, 5 and 6 describes the proposed multi-class arrhythmia
classification methodology followed by hardware implementation of the proposed
methodology in Section 7. Section 8 discusses all the results obtained from
the methodologies presented in the preceding sections. Finally, the conclusion
is summarised in Section 9.

##  2 Test Database

The present work uses four open source databases, such as; CPSC database and
CPSC-Extra database Liu et al. (2018), PTB-XL database Wagner et al. (2020),
the Georgia 12-lead ECG challenge (G12EC) database and Chapman-Shaoxing and
Ningbo database Zheng et al. (2020a), Zheng et al. (2020b), made available by
Physionet for the computing in cardiology challenge 2021 Reyna et al. (2021).

Each database contains a variety of 12-channel ECG records having different
arrhythmia conditions. In this work, the lead II ECG signal is utilized for
simulation purposes. Each of the signals is sampled at a rate of 500 Hz. The
ECG records with PAC and supraventricular premature beats (SVPB) are grouped
into a single class of PAC. Similarly, records of PVC and ventricular
premature beats (VPB) are merged as the same class PVC. A view of each
arrhythmia category of ECG signals is presented in Figure 1.

Figure 1: ECG records having different arrhythmia conditions.

The databases are used to accumulate single channel ECG records which form the
test dataset. Due to the degradation of the performance of the classifier, the
test dataset can be highly imbalanced with 33 records of PVC and 14,993
entries of NSR. Hence, the synthetic minority oversampling technique (SMOTE)
and random under-sampling are used to balance the dataset. Initially, SMOTE is
applied to increase the data of minority classes Chawla et al. (2002),
Sivapalan et al. (2022). Further, random under-sampling is used to reduce the
data in the majority class by randomly eliminating the records to obtain a
balancing in data Batista et al. (2004). A total of 31059 ECG records (3451
from each class) are considered in this work. The detailed information on each
class of ECG records is described in Table 1.

Table 1: Number of ECG records of each arrhythmia class Arrhythmia class | Original ECG records | After SMOTE | Final ECG records
---|---|---|---
NSR | 14993 | 14993 | 3451
SA | 1361 | 3451 | 3451
SB | 9353 | 9353 | 3451
STACH | 3451 | 3451 | 3451
AF | 1500 | 3451 | 3451
AFL | 1526 | 3451 | 3451
PAC | 538 | 3451 | 3451
1AVB | 754 | 3451 | 3451
PVC | 33 | 3451 | 3451
NSR - Normal sinus rhythm, SA - Sinus Aarrhythmia, SB - Sinus bradycardia,
STACH - Sinus tachycardia, AF - Atrial fibrillation, AFL - Atrial flutter, PAC
- Premature atrial contraction, 1AVB - 1st Degree AV Block, PVC - Premature
ventricular contraction

##  3 Proposed ECG Arrhythmia Classification Methodology

The proposed arrhythmia detection and classification methodology is
implemented in 3 parts - In the first part, heart rate variability and time
domain morphological features from ECG signals are extracted to train a
classifier. To improve the results obtained from this classification approach,
in the second part, heart rate variability features are paired with wavelet
coefficient features obtained from stationary wavelet decomposition of the ECG
signal. This gives significantly improved classification performance for each
of the 9 classes. In the last part, both the classification methods are
implemented in a hardware embedded device to make it viable for real time ECG
signal analysis and arrhythmia detection,

Flowchart of the proposed arrhythmia classification work is presented in
Figure 2. The detailed description of each step is discussed in the following
subsections.

Figure 2: Flowchart of the proposed arrhythmia detection scheme.

###  3.1 ECG Signal Pre-processing

The ECG records are often corrupted with baseline wander, muscle artifacts
etc. These noises mask the clinical components in the ECG signal which results
in poor classification performance Satija et al. (2018), Chatterjee et al.
(2020). Most of the clinical information of ECG signals is in the frequency
range of 1-150 Hz Walraven (2010). For the first part of the work, the ECG
signals are filtered through a band-pass filter with a cut-off frequency of
1-150 Hz, which will remove low out-band noise (baseline wander) as well as
high frequency noise (muscle artifacts). In the second part, for the wavelet
coefficient features, the coefficient sets with frequency range of baseline
wander and muscle artifacts are discarded and not used for feature extraction.
Powerline interference is removed in both parts of the work using a 50 Hz
notch filter. For time domain features, it is crucial to properly delineate
the local components (P wave, QRS complex and T wave) in an ECG signal. In
this work, ECGDeli, an open-source ECG delineation toolbox is used for
accurate delineation of the local components in ECG signals Pilia et al.
(2021). The onset, offset and peak of P waves, QRS complex and T waves are
reliably detected by this toolbox as shown in Fig 3

Figure 3: Local signal component delineation in ECG using ECGDeli toolbox.

##  4 Heart rate variability and time domain feature extraction

As per the description of arrhythmia in Walraven (2010), it can be concluded
that heart rate (HR), regularity of HR, P wave morphology, PR interval (PRI)
and QRS complex morphology are the important factors to be considered while
looking for the signs of arrhythmia in ECG records. Distinct variations are
caused by different types of arrhythmia in HR, P wave morphology, PRI timing
and QRS complex morphology. These characteristic variations are analyzed by
physicians to detect arrhythmia from ECG recordings and accurate assessment of
its type. Hence, in the first part of this work, keeping in consideration the
characteristic variations attributed to different types of arrhythmia,
detected fiducial points are used to extract the time domain and morphological
features which can be grouped under four categories - HR variability features,
P wave features, PRI features and QRS complex features. The features are a
combination of non-linear, higher-order statistical, entropy and energy-based
features so as to best capture the morphological and timing variations in the
ECG for different arrhythmia conditions. All categories of features are listed
in Table 2. The details regarding each category of features are described in
the following subsections.

Table 2: Category wise feature details.

Feature Category | Features
---|---
HRV Feature |  | Maximum HR ($HR_{max}$), Minimum HR ($HR_{min}$), Mean HR ($HR_{mean}$), Std. of HR ($HR_{std}$),
---
Max deviation of HR ($HR_{MaxDev}$), Mean of absolute difference of HR
($HR_{MAD}$),
Kurtosis of HR ($HR_{kurt}$), Skewness of HR ($HR_{skew}$), Approximate
entropy of HR ($HR_{ApEn}$),
Shannon entropy of HR ($HR_{ShEn}$), Permutation entropy of HR ($HR_{PeEn}$)
P wave Feature |  | P wave peak ($P_{peak}$), P wave width ($P_{width}$), Max Deviation of P wave ($P_{MaxDev}$),
---
P wave energy ($P_{energy}$), Correlation of P waves ($P_{Corr}$), Spectral
Entropy of P wave ($P_{SpEn}$),
Kurtosis of P wave ($P_{kurt})$, Skewness of P wave ($P_{skew}$), Atrial HR
($P_{atrialHR}$)
PRI Feature |  | Mean PR Interval ($PRI_{mean}$), Std. of PR Interval ($PRI_{std}$), Maximum PR Interval ($PRI_{max}$),
---
Minimum PR Interval ($PRI_{min}$), Max Deviation of PR Interval
($PRI_{MaxDev}$)
QRS Complex Feature |  | QRS Width ($QRS_{width}$), Correlation of QRS Complex ($QRS_{Corr}$), QRS complex energy ($QRS_{energy}$)
---
Spectral Entropy of QRS Complex ($QRS_{SpEn}$), Sample Entropy of QRS
($QRS_{SaEn}$),
Kurtosis of QRS complex ($QRS_{kurt})$, Skewness of QRS complex ($QRS_{skew}$)

###  4.1 Heart Rate Variability (HRV) Features

The analysis of HR is an imperative step in detecting arrhythmia from ECG
signals. Heart rate variation may contain signs of cardiac disease that are
already present or warnings of impending cardiac diseases. HR regularity and
variability analysis is a crucial step for determining heart rate rhythm as
well as arrhythmia Walraven (2010). In this work, a total of 11 features are
extracted from heart rate for analyzing the HR regularity and variability. HR
is obtained from the interval between two consecutive R peaks. R peaks in this
work are readily obtained from ECG deli. If $RR_{i}$ is the $i^{th}$ RR
interval, then for sampling frequency 500 Hz, the HR is calculated as:

$HR_{i}=\frac{60*500}{RR_{i}}\>bpm$ (1)

Mean, standard deviation, maximum and minimum HR are some of the linear
features calculated. Maximum deviation in HR is a feature calculated as the
difference in maximum and minimum HR. From the difference between consecutive
HR, the mean absolute difference of HR is also calculated.

$HR_{MeanAbsDiff}=\frac{1}{n}\sum_{i=1}^{n}|HR_{i+1}-HR_{i}|$ (2)

Higher order statistics feature such as kurtosis and skewness of HRV signal
are also considered and calculated as follows:

$HR_{Kurt}=\frac{\frac{1}{n}\sum_{i=1}^{n}(HR_{i}-HR_{mean})^{4}}{[\frac{1}{n}\sum_{i=1}^{n}(HR_{i}-HR_{mean})^{2}]^{2}}$
(3)

$HR_{Skew}=\frac{\frac{1}{n}\sum_{i=1}^{n}(HR_{i}-HR_{mean})^{3}}{[\frac{1}{n}\sum_{i=1}^{n}(HR_{i}-HR_{mean})^{2}]^{\frac{3}{2}}}$
(4)

Non-linear entropy based features such as approximate entropy ($HR_{ApEn}$),
permutation entropy ($HR_{PeEn}$) and Shannon entropy ($HR_{ShEn}$) are
calculated from HRV signal. These features are determined as follows:

For a time series of length $L$, window size of $m$ and threshold of $r$

$C_{i}^{m}(r)=\frac{n_{im}(r)}{L+m-1}$ (5)
$\emptyset^{m}(r)=\frac{\sum_{i=1}^{L-m+1}ln[C_{i}^{m}(r)]}{L-m+1}$ (6)

where, $C_{i}^{m}(r)$ is probability of similarity at threshold $r$,
$n_{im}(r)$ is no. of segments similar to $i^{th}$ segment with threshold $r$
and $\emptyset^{m}(r)$ is segment value.
Approximate entropy is determined as

$HR_{ApEn}=\emptyset^{m}(r)-\emptyset^{m+1}(r)$ (7)

Permutation entropy ($PeEn$) is a measure of complexity taking into account
temporal order of the successive points in a time series Riedl et al. (2013).
In $PeEn$ calculation total time series is presented into a group of pattern.
If $p_{k}$ is probability of occurrence of the $k_{th}$ pattern and $K$ is the
length of each pattern then $PeEn$ can be expressed as:

$HR_{PeEn}=-\sum_{k=1}^{K!}p_{k}log_{2}(p_{k})$ (8)
$HR_{ShEn}=-\sum_{i=1}^{n}HR_{norm}log_{2}(HR_{norm})$ (9)

where

$HR_{norm}=\frac{HR}{max(HR)}$ (10)

###  4.2 P wave features

P wave in ECG signal is the representation of the atrial activity of the heart
and signifies depolarisation of atria. Many of the arrhythmia types (AF, AFL,
PAC) affect the morphology of P waves in different ways Walraven (2010).
Hence, the morphological information from P wave features can be used
effectively to distinguish these arrhythmia classes. Amplitude, width, energy
and maximum deviation of P waves are extracted as features which capture
essential information like the presence, shape and morphology of P waves. The
correlation coefficient of the P wave analyzes the similarity of P waves in
consecutive beats. This feature can effectively detect the changes in P wave
morphology from beat to beat within an ECG segment in case of an ectopic beat.
Other non-linear, statistical features like spectral entropy, skewness, and
kurtosis are also calculated for P wave signal components which capture the
difference of P wave information for different arrhythmia conditions. Atrial
HR is also calculated from the P-P interval which is the time between two
consecutive P peaks. As multiple P waves exist throughout the entire ECG
segment hence mean and standard deviation are calculated for each of the
above-mentioned features. A total group of 18 features are extracted from the
P wave signal component.

###  4.3 PR interval features

PRI is the time duration between the onset of a P wave to the onset of the QRS
complex. It signifies the start of atrial depolarisation to the start of
ventricular depolarization. It is an important marker for the atrial activity
and conduction of electrical impulses through the AV node. The PRI value may
be prolonged for the rhythm generated other than sinus or heart block
condition. Hence, information on PRI can be used to distinguish between
arrhythmia classes. A total 5 PRI features are extracted such as mean PRI
($PRI_{mean}$), standard deviation of PRI ($PRI_{std}$), maximum PRI
($PRI_{max}$), minimum PRI ($PRI_{min}$) and maximum deviation of PRI
($PRI_{MaxDev}$).

###  4.4 QRS features

QRS complex represents the ventricular activity of the heart. The ventricular
rhythm arrhythmia like PVC affects the morphology of the QRS complex as
compared to normal sinus rhythm. Hence, effective features can be extracted
from QRS morphology for better classification. Different morphology, entropy
and higher order statics features such as QRS width, correlation coefficient
of consecutive QRS, spectral entropy, sample entropy, kurtosis and skewness
are extracted from QRS complexes. Similar to P waves multiple QRS complexes
exist throughout the entire ECG signal hence mean and standard deviation are
calculated for each of the above-mentioned features. A total group of 14
features are extracted from QRS complexes.

##  5 Wavelet coefficient feature extraction

To improve the methodology proposed in the first part, in the second part, an
arrhythmia detection scheme is devised where in addition to heart rate
features, wavelet domain features are obtained from coefficients of wavelet
decomposition of the ECG signals using SWT.

###  5.1 SWT based signal decomposition

Wavelet transform is a signal processing tool that allows the decomposition of
the signal into different time and frequency scales where each scale allows
analysis of various signal properties and characteristics. This tool for
analyzing non-stationary signals is useful and simple for identifying subtle
variations in the signal morphology over the scales of interest Asgari et al.
(2015b). A series of high and low-pass filters are used to analyze high and
low-frequency components of the signal. At each level, convolution of the
input signal with high pass filters gives the detail coefficients $D_{n}$ and
convolution with low pass filters gives approximation coefficients $A_{n}$.

SWT is a time-invariant discrete wavelet transform method. At each
decomposition level, SWT coefficients have the same number of samples as the
original signal, thus preserving the temporal information of the signal and
overcoming the problem of repeatability and robustness which exists with
discrete wavelet transformAsgari et al. (2015b). SWT of a signal $x[n]$ gives
the coefficients $c_{i,j}$ obtained using equation Merah et al. (2015) -

$c_{i,j}=\sum_{n\in z}x[n]\psi^{*}_{i,j}[n]$ (11)

where $\psi_{i,j}$ is a translated and scaled version of the mother wavelet
$\psi_{0,0}$

$\psi_{i,j}[k]=2^{-(i/2)}\psi_{0,0}(2^{-i}(k-j))$ (12)

To implement $L$-level SWT on a signal, the length of the signal should be a
multiple of $2^{L}$. Signals with lesser samples are zero-padded to make the
signal length a multiple of $2^{L}$. At each successive decomposition level,
the impulse response of the high and low pass filters are upsampled by a
factor of 2 giving a coefficient series with the same temporal resolution as
the original signal.

In this work, the ECG records are sampled at 500 Hz with a maximum frequency
component of 250 Hz. Considering the frequency ranges of the QRS complex, P
wave and T wave, a 7-level SWT decomposition is chosen to be applied to the
ECG records using Symlet-7 wavelet. Symlet-7 wavelet is chosen as the mother
wavelet because of its close similarity to ECG signal morphology and is
extensively used in different ECG signal processing-based works Ansari et al.
(2017). The frequency range of each decomposition level is shown in Figure 4.
Detail coefficients D3, D4, D5 and D6 correspond to the frequency range of the
QRS complex while D6 and D7 correspond to the P and T waves of the ECG signal.
Thus, D3, D4, D5, D6 and D7 are considered for wavelet coefficient feature
calculation in the subsequent steps of the algorithm.

Figure 4: Frequency range information of SWT decomposition sub-band levels.

Figure 5: Original ECG with decomposed signals (a) Original ECG ; (b)-(h)
Detail coefficients D1 - D7 ; (i) Approximate coefficient A7.

###  5.2 HRV features

In addition to the 11 HR features from the Subsection 4.1, 5 more features
were added in this scheme to increase the arrhythmia classification accuracy
as HR features contribute the most in differentiating multiple classes of
arrhythmia. These new features are - standard deviation of absolute
difference, coefficient of variance (CoV) of HR, Higuchi fractal dimension,
Hjorth mobility and Hjorth complexity. CoV of HR is the ratio of the standard
deviation of HR to the mean of HR. Higuchi fractal dimension is a
computational method to determine changes in a signal from the measure of its
complexity Chinara et al. (2021). For a finite set of time series observations
taken at regular intervals $X(1),X(2),X(3)....X(N)$, a new time sub series can
be constructed by taking

$X_{k}^{m}:X(m),X(m+k),X(m+2k),...,X(m+[\frac{N-m}{k}].k)$ (13)

where $m=1,2,3,..,k$ is the start of every sub series and $k$ is the interval.
This gives a total of $k$ new sub series. The length of each of the sub series
is given by

$L_{m}(k)=\frac{{\lfloor{\sum_{i=1}^{[\frac{N-m}{k}]}x_{m+id}-x_{m+(i-1)d}}}\rfloor\frac{N-1}{\frac{N-m}{k}.k}}{k}$
(14)

On plotting the average value of the length over $k$ sub series against $k$ in
a double logarithmic scale, the data falls on a straight line with slope-D.
The value $D$ is Higuchi fractal dimension for the given time series Higuchi
(1988).

Hjorth mobility ($HR_{HM}$) is a Hjorth descriptor that describes the mean
frequency of a signal. Hjorth complexity ($HR_{HC}$) denotes an estimation of
the signal bandwidth from the ratio of the peak value to the harmonic content
of the signal Hjorth (1970). These two features are calculated as:

$HR_{HM}=\frac{\sigma^{\prime}_{x}}{\sigma_{x}}$ (15)

$HR_{HC}=\frac{\frac{\sigma^{\prime\prime}_{x}}{\sigma^{\prime}_{x}}}{\frac{\sigma^{\prime}_{x}}{\sigma_{x}}}$
(16)

where, $\sigma_{x}$ is the variance of the time series, $\sigma^{\prime}_{x}$
is the 1st derivative of the variance and $\sigma^{\prime\prime}_{x}$ is the
$2^{nd}$ derivative of the variance.

###  5.3 Wavelet coefficient features

Detail coefficients D3, D4, D5, D6 and D7 correspond to the frequency range of
clinically significant information in ECG signal. 10 features are extracted
from each of these four sets of coefficients to give a total of 50 features.
The feature set includes statistical, entropy and energy based features.
Statistical features such as mean, standard deviation, skewness and kurtosis
are obtained from each set of detail coefficients. Entropy features of
approximate entropy, Shannon entropy, permutation entropy and log energy
entropy ($LEEn$), when applied to wavelet coefficients, capture the
complexity, regularity and uncertainty in the wavelet decomposed subbands of
the ECG signal. Two energy based features namely: relative wavelet energy
($RWE$) and mean wavelet energy ($MWE$) are also considered. The calculation
of $LEEn$, $RWE$ and $MWE$ are as follows Kumar et al. (2018):

$LE_{En}=\sum_{i=1}^{N}log(x_{i}^{2})$ (17)

where, $x_{j}$ is the $i^{th}$ sample and $N$ is the length of the sub-band
signal.

$RWE=\frac{\sum_{i=1}^{N}C_{j}(i)^{2}}{\sum_{j=1}^{L}\sum_{i=1}^{N}C_{j}(i)^{2}}$
(18)

$MWE=\frac{\sum_{i=1}^{N}C_{j}(i)^{2}}{N}$ (19)

where, $N$ is the total number of coefficients in $j^{th}$ level and $L$ is
the total number of decomposition levels.

Table 3: Set of wavelet coefficient features Wavelet features
---
Mean ($CD_{mean}$)
Standard deviation ($CD_{std}$)
Skewness ($CD_{skew}$)
Kurtosis ($CD_{kurt}$)
Approximate entropy ($CD_{ApEn}$)
Shannon entropy ($CD_{ShEn}$)
Permutation entropy ($CD_{PeEn}$)
Log energy entropy ($CD_{LE_{En}}$)
Relative wavelet energy ($CD_{MWE}$)
Mean wavelet energy ($CD_{RWE}$)

##  6 Machine Learning based Arrhythmia Classifier

In this section, a machine learning based classifier is utilized for the
automated detection of arrhythmia. The supervised machine learning based
approaches require feature matrix along with labels. In part one of the work,
a total of 48 features (11 HRV, 18 P wave, 5 PRI and 14 QRS features) made up
the feature set and in next part, there were 66 total features (16 HRV and 50
wavelet features) in the feature set. Feature sets obtained in the previous
sections are supplied to the classifier. There are several machine learning-
based techniques and each of them is extensively used in many literature.
Random forest (RF) is one of the popular ensemble classifiers. It is
extensively employed in many classification problem including medical signal
and image processing Kung et al. (2020), Panayides et al. (2020). RF is an
ensemble classier which is made from bootstrap aggregation of multiple
decision trees. Each decision tree independently generates an output as per
the input feature matrix supplied to it. Finally, the net resultant output can
be found by applying a voting strategy on the outputs from multiple trees. The
aggregation of voting makes RF more effective classifier and less susceptible
to outliers and noises Shaikhina et al. (2019). A comparative analysis on
different machine learning classifiers is studied in the following subsection
and it can be checked that the RF performs better compared to others. Hence,
RF is used here as the arrhythmia classifier. Hyper-parameter tuning is done
separately for both features sets to obtain a set of optimal parameters for
the RF classifier giving the best performance.

##  7 Hardware implementation of the proposed classification schemes

To validate the proposed arrhythmia detection and classification schemes, this
section presents the hardware implementation using a Raspberry Pi 4 model B. A
programmable system on a chip like Raspberry Pi can handle the complexity of
computation while keeping the power consumption low. Raspberry Pi is a single
board computer working on a Linux based operating system and takes real time
input data which can then be used for a multitude of applications. Apart from
this low cost Pi, a computer for feature extraction from dataset and a display
monitor to view the output from the Raspberry Pi is used as hardware
components. The ECG dataset consisting of 31,059 10-sec ECG records is given
as input to the feature extraction algorithm. Feature extraction in both the
proposed schemes is done using Matlab 2022b. The obtained feature set is given
as input to the Raspberry Pi where the classifier is deployed. The classifier
model is written in Python programming language using scikit-learn library.
The classifier is trained on the feature set and the trained model then can be
used to detect and classify arrhythmic rhythms in real time ECG signals.

Figure 6: Hardware implementation with Raspberry Pi 4B

Figure 7: Feature extraction is fed to Raspberry Pi 4B (left) for feature
extraction and a display monitor (right) to show output obtained from
classifier loaded on Raspberry Pi 4B

##  8 Results and Discussion

Performance of both the proposed cardiac arrhythmia detection scheme is
presented in this section. Training and testing of the classifier is done
through 10-fold cross-validation. In 10-fold cross-validation, the entire
dataset is divided into 10 subsets. At each time, one subset is used for
testing purposes and the remaining nine subsets are used for training the
algorithm. This process is repeated for all 10 sub-sets and finally, the
average values of the performance metrics are reported. The performance of the
proposed arrhythmia detection method is evaluated using standard metrics:
accuracy ($Acc$), sensitivity ($Se$), precision ($+P$) and $F1$ score. These
are calculated from parameters: true positives ($TP$), true negative ($TN$),
false positive ($FP$) and false negative ($FN$).

$Accuracy(Acc)=\frac{TP+TN}{TP+TN+FP+FN}$ (20)

$Sensitivity(Se)=\frac{TP}{TP+FN}$ (21)

$Precision(+P)=\frac{TP}{TP+FP}$ (22)

$F1Score=\frac{2*TP}{2*TP+FP+FN}$ (23)

The detailed 10-fold cross-validation performance result of classification
with HRV, P wave, PRI and QRS complex morphological features are presented in
Table 4 . Average performance obtained is $Acc$ of 85.11%, $Se$ of 85.11%, $P$
of 85.07% and $F1$ score of 85%.

Table 4: 10- fold Cross-validation performance results with heart rate and time domain features Fold | Acc(%) | Se(%) | +P(%) | F1(%)
---|---|---|---|---
1 | 85.64 | 85.64 | 85.65 | 85.57
2 | 85.12 | 85.12 | 85.02 | 85.00
3 | 84.86 | 84.87 | 84.81 | 84.77
4 | 84.83 | 84.84 | 84.79 | 84.68
5 | 85.57 | 85.57 | 85.55 | 85.46
6 | 85.47 | 85.47 | 85.31 | 85.26
7 | 85.44 | 85.44 | 85.46 | 85.34
8 | 85.02 | 85.02 | 84.98 | 84.93
9 | 84.93 | 84.92 | 84.95 | 84.86
10 | 84.21 | 84.21 | 84.15 | 84.09
Average | 85.11 | 85.11 | 85.07 | 85.00

Class-wise classification performance of the described method is presented in
Table 5. Maximum classification $F1$ score of 96.36% is obtained for STACH.
Lowest $F1$ score of 75.31% and 71.60% is obtained for the classes of AF and
AFL.

Table 5: Class wise performance of arrhythmia classification with heart rate and time domain features Class | Se(%) | +P(%) | F1(%)
---|---|---|---
NSR | 79.98 | 92.00 | 85.57
SA | 86.58 | 86.56 | 86.57
SB | 97.19 | 89.04 | 92.93
STACH | 97.74 | 95.01 | 96.36
AF | 75.02 | 75.59 | 75.31
AFL | 70.65 | 72.58 | 71.60
PAC | 79.83 | 79.76 | 79.80
1AVB | 83.14 | 83.74 | 83.44
PVC | 96.35 | 91.70 | 93.97

To have a more in depth analysis of the feature set with HRV features and
three categories of time domain and morphological features, performance of
different subsets of the feature set is compared in Table 6. It can be
observed that the HRV feature has the superior effectiveness of detecting
arrhythmia with a maximum classification $Acc$ of 82.81%. It solidifies the
importance of HR features for arrhythmia classification. PRI features alone
are the least accurate in multi-class arrhythmia classification as the PRI
segment mostly signifies the conduction time of impulse through the AV node.
Similarly, P wave features and QRS complex features also provide less
classification accuracy compared to the HRV feature set when used
individually. Adding P-wave, PRI and QRS complex features along with HRV
features makes the classification of nine classes more robust and accurate. It
is also established that each of the four category features adds more
robustness to the algorithm and makes it efficient in multi-class arrhythmia
classification. This gives important insights into the clinical relevance of
different ECG signal components.

Table 6: Performance comparison with different subsets of HRV and time domain features. Feature Category | Acc(%) | Se(%) | P(%) | F1(%)
---|---|---|---|---
HRV | 82.81 | 82.81 | 82.64 |  82.61
P Wave | 73.10 | 73.10 | 72.87 | 72.69
PRI | 43.65 | 43.65 | 42.83 | 43.06
QRS Complex | 59.19 | 59.19 | 59.58 | 58.81
HRV + P Wave | 80.48 | 80.48 | 80.18 | 80.15
HRV + P Wave + PRI | 82.46 | 82.46 | 82.43 | 82.33
HRV + P Wave + PRI + QRS | 85.11 | 85.11 | 85.07 | 85.00

Next, results of classification with HRV and wavelet features are presented in
Table 7. 10-fold cross-validation with the proposed scheme gives an average
$Acc$ of 90.91% , $Se$ of 90.91%, $+P$ of 90.96% and $F1$ score of 90.87%.
Both the methods show similar results for each validation subset which
suggests the robust classification performance of the proposed schemes.

Table 7: Detailed cross-validation results with heart rate and wavelet features Fold | Acc(%) | Se(%) | +P(%) | F1(%)
---|---|---|---|---
1 | 90.98 | 90.98 | 91.09 | 90.98
2 | 91.21 | 91.21 | 91.25 | 91.18
3 | 90.18 | 90.18 | 90.20 | 90.13
4 | 91.07 | 91.07 | 91.05 | 91.02
5 | 90.62 | 90.62 | 90.66 | 90.59
6 | 91.36 | 91.37 | 91.34 | 91.32
7 | 91.69 | 91.69 | 91.80 | 91.66
8 | 90.62 | 90.62 | 90.68 | 90.59
9 | 90.59 | 90.59 | 90.67 | 90.53
10 | 90.78 | 90.78 | 90.80 | 90.73
Average | 90.91 | 90.91 | 90.96 | 90.87

The performance of the proposed scheme for individual arrhythmia class is
described in Table 8. The method shows a maximum classification $F1$ score of
98.00% for PVC. The $F1$ score of AFL arrhythmia class is marginally small at
80.81% compared to other classes. The ECG morphology characteristics of AF and
AFL arrhythmia classes are quite similar, hence for these two classes, both
the methods give larger false detection compared to other classes. This
results in a degradation of performance for AFL arrhythmia classification.

Table 8: Class-wise performance result with heart rate and wavelet features Class | Se(%) | +P(%) | F1(%)
---|---|---|---
NSR | 87.54 | 92.08 | 89.75
SA | 91.97 | 88.09 | 89.99
SB | 97.74 | 93.07 | 95.35
STACH | 98.12 | 95.49 | 96.78
AF | 79.27 | 87.11 | 83.01
AFL | 82.61 | 79.08 | 80.81
PAC | 88.03 | 90.23 | 89.12
1AVB | 85.71 | 88.80 | 87.23
PVC | 99.57 | 96.49 | 98.00

Proper selection of the mother wavelet in SWT-based decomposition is a crucial
task. Generally, a wavelet is chosen whose shape is nearly similar to the
morphology of the ECG cycle. Different mother wavelets such as Daubechies,
Coiflet, Symlet, Bi-orthogonal etc are extensively used in several ECG signal
processing-related works. In this work, a comparative performance study has
been carried out using different mother wavelets. As presented in Table 9, the
Symlet wavelet of order 7 (Symlet 7) shows a better performance of accuracy of
90.91%.

Table 9: Performance comparison of the proposed scheme for different mother wavelets Wavelet | Acc(%) | Se(%) | +P(%) | F1(%)
---|---|---|---|---
Haar | 89.97 | 89.97 | 89.98 | 89.90
Daubechies 6 | 90.41 | 90.41 | 90.40 | 90.35
Coiflet 2 | 90.86 | 90.86 | 90.89 | 90.82
Biorthogonal 4.4 | 90.90 | 90.90 | 90.93 | 90.86
Symlet 7 | 90.91 | 90.91 | 90.96 | 90.87

The arrhythmia classification performance of HRV features along with different
combinations of wavelet coefficient features is presented in Table 10. It can
be observed that the HRV features alone have an effective classification
accuracy of 82.81%. Wavelet features of D3, D4, D5, D6 and D7 coefficient set
give an accuracy of 83.04% without including HRV features. Performance
gradually increases on combining the HRV features and wavelet coefficient
feature sets. This justifies the effectiveness of HRV features along with
wavelet coefficient features for the detection of multi-class cardiac
arrhythmia.

Table 10: Performance comparison of proposed scheme for different feature combinations Feature Category | Acc(%) | Se(%) | +P(%) | F1(%)
---|---|---|---|---
HRV | 82.81 | 82.81 | 82.64 | 82.61
D3 + D4 + D5 + D6 + D7 | 83.04 | 83.03 | 82.88 | 82.73
HRV + D3 | 89.45 | 89.45 | 89.43 | 89.40
HRV + D3 + D4 | 89.57 | 89.57 | 89.54 | 89.52
HRV + D3 + D4 + D5 | 89.62 | 89.62 | 89.63 | 89.58
HRV + D3 + D4 + D5 + D6 | 90.45 | 90.45 | 90.47 | 90.41
HRV+ D3 + D4 + D5 + D6 + D7 | 90.91 | 90.91 | 90.96 | 90.87

The machine learning-based classifier for arrhythmia detection is an
imperative component in this work. Output of the classification is highly
dependent on the proper selection of a machine learning classifier. In this
section, a comparative analysis is studied on the performance of both schemes
for different machine learning classifiers. The techniques are tested on three
other extensively used classifiers namely: K- nearest neighbours (KNN),
decision tree (DT), support vector machine (SVM) and RF. As presented in Table
11, maximum accuracy is obtained with RF classifier. Aggregation of the voting
concept in RF improves the classification performance compared to other
classifiers.

Table 11: Performance comparison with different machine learning classifiers Classifier | HRV and time domain features | HRV and wavelet features
---|---|---
Acc(%) | Se(%) | +P(%) | F1(%) | Acc(%) | Se(%) | +P(%) | F1(%)
KNN | 72.97 | 72.97 | 73.01 | 72.77 | 79.29 | 79.28 | 79.16 | 79.08
Decision Tree | 74.22 | 74.22 | 74.29 | 74.23 | 81.76 | 81.76 | 81.80 | 81.75
SVM | 81.30 | 81.30 | 81.21 | 81.13 | 85.76 | 85.76 | 85.81 | 85.64
Random Forest | 85.11 | 85.11 | 85.07 | 85.00 | 90.91 | 90.91 | 90.96 | 90.87

##  9 Conclusion

In this work, multi-class cardiac arrhythmia detection schemes are proposed.
In the first part, HRV features are incorporated together with time domain
statistical, entropy and higher order statistical features of P wave, PR
interval and QRS complex for the effective classification of cardiac
arrhythmia using single-channel ECG records. A set of total 48 features are
applied to a machine learning-based random forest classifier. The detailed
10-fold cross-validation results show that the proposed multi-class cardiac
arrhythmia detection algorithm effectively classifies nine rhythms with an
average Acc of 85.11%, Se of 85.11%, P of 85.07% and F1 score of 85.00%. In
the second part, wavelet coefficient features are used along with HRV features
to successfully classify different arrhythmia types. Initially, stationary
wavelet transform is applied to decompose the ECG signal into different sub-
band levels. Considering the frequency localization property of wavelet
transform, matching with the frequency of ECG local components, the detail
coefficients D3, D4, D5, D6, and D7 are further processed for feature
extraction. A set of 66 features based on timing information, entropy, higher-
order statistics, and energy are extracted and applied to RF classifier.
Detailed cross-validation results show that the proposed multi-class cardiac
arrhythmia detection algorithm can effectively classify nine rhythm type with
average Acc of 90.91%, Se of 90.91%, +P of 90.96% and F1 score of 90.87%. For
both parts of the work, ECG records of broadly distributed four standard
databases of the Physionet Challenge 2021 are combined to prepare a test
database having nine classes of arrhythmia. Lastly, both the classification
schemes are implemented on Raspberry Pi. It’s low power consumption, light
weight and compact design makes it suitable for application in real time
monitoring and processing of ECG signals. A close observation of the
simulation results affirm that the proposed schemes can effectively be
utilized in an advanced automated cardiac disease monitoring system.

##  10 Acknowledgment

The authors would like to thank Dr. Ankita Pramanik, Assistant Professor,
Electronics and Telecommunication Engineering Department, IIEST, Shibpur, for
her guidance and support.

## References

  * Tsao et al. [2022] Connie W Tsao, Aaron W Aday, Zaid I Almarzooq, Alvaro Alonso, Andrea Z Beaton, Marcio S Bittencourt, Amelia K Boehme, Alfred E Buxton, April P Carson, Yvonne Commodore-Mensah, et al.  Heart disease and stroke statistics—2022 update: a report from the american heart association.  _Circulation_ , 145(8):e153–e639, 2022.
  * Roth et al. [2020] Gregory A Roth, George A Mensah, Catherine O Johnson, Giovanni Addolorato, Enrico Ammirati, Larry M Baddour, Noël C Barengo, Andrea Z Beaton, Emelia J Benjamin, Catherine P Benziger, et al.  Global burden of cardiovascular diseases and risk factors, 1990–2019: update from the gbd 2019 study.  _Journal of the American College of Cardiology_ , 76(25):2982–3021, 2020.
  * b2 [2021] Jun 2021.  URL https://www.who.int/news-room/fact-sheets/detail/cardiovascular-diseases-(cvds).
  * Faust and Ng [2016] Oliver Faust and Eddie YK Ng.  Computer aided diagnosis for cardiovascular diseases based on ecg signals: A survey.  _Journal of Mechanics in Medicine and Biology_ , 16(01):1640001, 2016.
  * Hong et al. [2022] Shenda Hong, Wenrui Zhang, Chenxi Sun, Yuxi Zhou, and Hongyan Li.  Practical lessons on 12-lead ecg classification: meta-analysis of methods from physionet/computing in cardiology challenge 2020.  _Frontiers in Physiology_ , 12:2505, 2022.
  * Bertsimas et al. [2021] Dimitris Bertsimas, Luca Mingardi, and Bartolomeo Stellato.  Machine learning for real-time heart disease prediction.  _IEEE Journal of Biomedical and Health Informatics_ , 25(9):3627–3637, 2021.
  * Ran et al. [2022] Shaolin Ran, Xiaoyun Yang, Ming Liu, Yong Zhang, Cheng Cheng, Hongling Zhu, and Ye Yuan.  Homecare-oriented ecg diagnosis with large-scale deep neural network for continuous monitoring on embedded devices.  _IEEE Transactions on Instrumentation and Measurement_ , 71:1–13, 2022.
  * Walraven [2010] Gail Walraven.  _Basic arrhythmias_.  Pearson Education, Upper Saddle River, NJ, 7th ed edition, 2010.  ISBN 9780135002384.
  * García et al. [2016] Manuel García, Juan Ródenas, Raúl Alcaraz, and José J Rieta.  Application of the relative wavelet energy to heart rate independent detection of atrial fibrillation.  _computer methods and programs in biomedicine_ , 131:157–168, 2016.
  * Asgari et al. [2015a] Shadnaz Asgari, Alireza Mehrnia, and Maryam Moussavi.  Automatic detection of atrial fibrillation using stationary wavelet transform and support vector machine.  _Computers in biology and medicine_ , 60:132–142, 2015a.
  * Rahul et al. [2021] Jagdeep Rahul, Marpe Sora, Lakhan Dev Sharma, and Vijay Kumar Bohat.  An improved cardiac arrhythmia classification using an rr interval-based approach.  _Biocybernetics and Biomedical Engineering_ , 41(2):656–666, 2021.
  * Mjahad et al. [2017] Azeddine Mjahad, A Rosado-Muñoz, Manuel Bataller-Mompeán, JV Francés-Víllora, and JF Guerrero-Martínez.  Ventricular fibrillation and tachycardia detection from surface ecg using time-frequency representation images as input dataset for machine learning.  _Computer methods and programs in biomedicine_ , 141:119–127, 2017.
  * Acharya et al. [2016] U Rajendra Acharya, Hamido Fujita, Muhammad Adam, Oh Shu Lih, Tan Jen Hong, Vidya K Sudarshan, and Joel EW Koh.  Automated characterization of arrhythmias using nonlinear features from tachycardia ecg beats.  In _2016 IEEE international conference on systems, man, and cybernetics (SMC)_ , pages 000533–000538. IEEE, 2016.
  * Malik et al. [2021] Junaid Malik, Ozer Can Devecioglu, Serkan Kiranyaz, Turker Ince, and Moncef Gabbouj.  Real-time patient-specific ecg classification by 1d self-operational neural networks.  _IEEE Transactions on Biomedical Engineering_ , 69(5):1788–1801, 2021.
  * Elhaj et al. [2016] Fatin A Elhaj, Naomie Salim, Arief R Harris, Tan Tian Swee, and Taqwa Ahmed.  Arrhythmia recognition and classification using combined linear and nonlinear features of ecg signals.  _Computer methods and programs in biomedicine_ , 127:52–63, 2016.
  * Liu et al. [2018] Feifei Liu, Chengyu Liu, Lina Zhao, Xiangyu Zhang, Xiaoling Wu, Xiaoyan Xu, Yulin Liu, Caiyun Ma, Shoushui Wei, Zhiqiang He, et al.  An open access database for evaluating the algorithms of electrocardiogram rhythm and morphology abnormality detection.  _Journal of Medical Imaging and Health Informatics_ , 8(7):1368–1373, 2018.
  * Wagner et al. [2020] Patrick Wagner, Nils Strodthoff, Ralf-Dieter Bousseljot, Dieter Kreiseler, Fatima I Lunze, Wojciech Samek, and Tobias Schaeffter.  Ptb-xl, a large publicly available electrocardiography dataset.  _Scientific data_ , 7(1):154, 2020.
  * Zheng et al. [2020a] Jianwei Zheng, Jianming Zhang, Sidy Danioko, Hai Yao, Hangyuan Guo, and Cyril Rakovski.  A 12-lead electrocardiogram database for arrhythmia research covering more than 10,000 patients.  _Scientific data_ , 7(1):48, 2020a.
  * Zheng et al. [2020b] Jianwei Zheng, Huimin Chu, Daniele Struppa, Jianming Zhang, Sir Magdi Yacoub, Hesham El-Askary, Anthony Chang, Louis Ehwerhemuepha, Islam Abudayyeh, Alexander Barrett, et al.  Optimal multi-stage arrhythmia classification approach.  _Scientific reports_ , 10(1):2898, 2020b.
  * Reyna et al. [2021] Matthew A Reyna, Nadi Sadr, Erick A Perez Alday, Annie Gu, Amit J Shah, Chad Robichaux, Ali Bahrami Rad, Andoni Elola, Salman Seyedi, Sardar Ansari, et al.  Will two do? varying dimensions in electrocardiography: the physionet/computing in cardiology challenge 2021.  In _2021 Computing in Cardiology (CinC)_ , volume 48, pages 1–4. IEEE, 2021.
  * Chawla et al. [2002] Nitesh V Chawla, Kevin W Bowyer, Lawrence O Hall, and W Philip Kegelmeyer.  Smote: synthetic minority over-sampling technique.  _Journal of artificial intelligence research_ , 16:321–357, 2002.
  * Sivapalan et al. [2022] Gawsalyan Sivapalan, Koushik Kumar Nundy, Soumyabrata Dev, Barry Cardiff, and Deepu John.  Annet: A lightweight neural network for ecg anomaly detection in iot edge sensors.  _IEEE Transactions on Biomedical Circuits and Systems_ , 16(1):24–35, 2022.
  * Batista et al. [2004] Gustavo EAPA Batista, Ronaldo C Prati, and Maria Carolina Monard.  A study of the behavior of several methods for balancing machine learning training data.  _ACM SIGKDD explorations newsletter_ , 6(1):20–29, 2004.
  * Satija et al. [2018] Udit Satija, Barathram Ramkumar, and M Sabarimalai Manikandan.  A review of signal processing techniques for electrocardiogram signal quality assessment.  _IEEE reviews in biomedical engineering_ , 11:36–52, 2018\.
  * Chatterjee et al. [2020] Shubhojeet Chatterjee, Rini Smita Thakur, Ram Narayan Yadav, Lalita Gupta, and Deepak Kumar Raghuvanshi.  Review of noise removal techniques in ecg signals.  _IET Signal Processing_ , 14(9):569–590, 2020\.
  * Pilia et al. [2021] Nicolas Pilia, Claudia Nagel, Gustavo Lenis, Silvia Becker, Olaf Dössel, and Axel Loewe.  Ecgdeli-an open source ecg delineation toolbox for matlab.  _SoftwareX_ , 13:100639, 2021.
  * Riedl et al. [2013] Maik Riedl, Andreas Müller, and Niels Wessel.  Practical considerations of permutation entropy: A tutorial review.  _The European Physical Journal Special Topics_ , 222(2):249–262, 2013.
  * Asgari et al. [2015b] Shadnaz Asgari, Alireza Mehrnia, and Maryam Moussavi.  Automatic detection of atrial fibrillation using stationary wavelet transform and support vector machine.  _Computers in biology and medicine_ , 60:132–142, 2015b.
  * Merah et al. [2015] Mostefa Merah, TA Abdelmalik, and BH Larbi.  R-peaks detection based on stationary wavelet transform.  _Computer methods and programs in biomedicine_ , 121(3):149–160, 2015.
  * Ansari et al. [2017] Sardar Ansari, Negar Farzaneh, Marlena Duda, Kelsey Horan, Hedvig B Andersson, Zachary D Goldberger, Brahmajee K Nallamothu, and Kayvan Najarian.  A review of automated methods for detection of myocardial ischemia and infarction using electrocardiogram and electronic health records.  _IEEE reviews in biomedical engineering_ , 10:264–298, 2017\.
  * Chinara et al. [2021] Suchismitha Chinara et al.  Automatic classification methods for detecting drowsiness using wavelet packet transform extracted time-domain features from single-channel eeg signal.  _Journal of neuroscience methods_ , 347:108927, 2021.
  * Higuchi [1988] Tomoyuki Higuchi.  Approach to an irregular time series on the basis of the fractal theory.  _Physica D: Nonlinear Phenomena_ , 31(2):277–283, 1988.
  * Hjorth [1970] Bo Hjorth.  Eeg analysis based on time domain properties.  _Electroencephalography and clinical neurophysiology_ , 29(3):306–310, 1970.
  * Kumar et al. [2018] Mohit Kumar, Ram Bilas Pachori, and U Rajendra Acharya.  Automated diagnosis of atrial fibrillation ecg signals using entropy features extracted from flexible analytic wavelet transform.  _Biocybernetics and Biomedical Engineering_ , 38(3):564–573, 2018.
  * Kung et al. [2020] Bo-Han Kung, Po-Yuan Hu, Chiu-Chang Huang, Cheng-Che Lee, Chia-Yu Yao, and Chieh-Hsiung Kuan.  An efficient ecg classification system using resource-saving architecture and random forest.  _IEEE Journal of Biomedical and Health Informatics_ , 25(6):1904–1914, 2020.
  * Panayides et al. [2020] Andreas S Panayides, Amir Amini, Nenad D Filipovic, Ashish Sharma, Sotirios A Tsaftaris, Alistair Young, David Foran, Nhan Do, Spyretta Golemati, Tahsin Kurc, et al.  Ai in medical imaging informatics: current challenges and future directions.  _IEEE journal of biomedical and health informatics_ , 24(7):1837–1857, 2020.
  * Shaikhina et al. [2019] Torgyn Shaikhina, Dave Lowe, Sunil Daga, David Briggs, Robert Higgins, and Natasha Khovanova.  Decision tree and random forest models for outcome prediction in antibody incompatible kidney transplantation.  _Biomedical Signal Processing and Control_ , 52:456–462, 2019.


# Periods of Hilbert Modular forms, Kronecker series and Cohomology

YoungJu Choie  YoungJu Choie Department of MathematicsPohang University of
Science and Technology (POSTECH) Pohang, Republic of Korea<EMAIL_ADDRESS>
###### Abstract.

Generalizing a result of [24, 8] about elliptic modular forms, we give a
closed formula for the sum of all Hilbert Hecke eigenforms over a totally real
number field with strict class number $1$, multiplied by their period
polynomials, as a single product of the Kronecker series.

###### Key words and phrases:

Hilbert modular form, parabolic cohomology, period polynomial, Kronecker
series

###### 2000 Mathematics Subject Classification:

11F41, 11F50, 11F60, 11F67

This work was partially supported by NRF 2018R1A4A1023590 and NRF
2017R1A2B2001807

##  1\.  Introduction

Based on Bol’s result [3] Eichler initiated a theory of the periods of
integrals so that an automorphic form of the first or second kind leads to a
cohomology class in the mapping of a Fuchsian group into a polynomial module
and the (converse) correspondence of each such cohomology class leads to an
automorphic form in one complex variable [9]. Shimura extended this theory by
showing that the structure of an abelian variety in certain cases can be also
given to the periods of such integrals and showed critical values of the
$L$-functions attached to elliptic modular forms can be computed explicitly
using the cohomology group [20]. This method was developed by Manin [16] who
proved an algebraic theorem for the periods of elliptic cusp forms for the
full modular group and studied $p$-adic properties of the algebraic factors in
$L$-functions. Kohnen-Zagier [14] further extended this theory to elliptic
modular forms including Eisenstein series [14] and studied forms whose period
polynomials have arithmetically interesting rational structure relating to
Bernoulli numbers, binary quadratic forms, zeta-functions of real quadratic
fields, modular forms of half-integral weight and Hilbert modular forms.
Hence, period polynomials, which allow us to compute the critical values of
$L$-function of modular forms at once, give a rich source of relations between
modular forms and arithmetic.

The period polynomial of an elliptic cusp form $f(\tau)=\sum_{\ell\geq
1}a_{f}(\ell)q^{\ell}\,\,(\tau\in\mathbb{H}=$ upper half plane, $q=e^{2\pi
i\tau})$ of weight $k$ on $SL_{2}(\mathbb{Z})$ is the polynomial of degree
$k-2$ defined by

(1.1) $\displaystyle r_{f}(X)=\int_{0}^{i\infty}f(\tau)(\tau-X)^{k-2}d\tau$

or equivalently by

$r_{f}(X)=-\sum_{n=0}^{k-2}\frac{(k-2)!}{(k-2-n)!}\frac{L(f,n+1)}{(2\pi
i)^{n+1}}X^{k-2-n},$

where $L(f,s)=\sum_{n\geq 1}\frac{a_{f}(n)}{n^{s}}\,(Re(s)\gg 0).$ The maps
$f\rightarrow r_{f}^{ev}$ and $f\rightarrow r_{f}^{od}$ assigning to $f$ the
even and odd parts of $r_{f}$ are both injective with known images from the
Eichler-Shimura-Manin theory.

When $f$ is a Hecke eigenform then one has the two-variable polynomial

$\displaystyle
r_{f}(X,Y):=\frac{r_{f}^{ev}(X)r_{f}^{od}(Y)+r_{f}^{od}(X)r_{f}^{ev}(Y)}{(2i)^{k-3}<f,\,f>}\in\mathbb{Q}_{f}[X,Y]$

where $\mathbb{Q}_{f}$ is the field generated by Fourier coefficients of $f$
over $\mathbb{Q}.$

Zagier [24] found the following attractive formula :

$\displaystyle\frac{(XY-1)(X+Y)}{X^{2}Y^{2}}T^{-2}+\sum_{k\geq
2}\sum_{f\in\mathcal{B}_{k}}r_{f}(X,Y)f(\tau)\frac{T^{k-2}}{(k-2)!}$
$\displaystyle=$ $\displaystyle
F_{\tau}(T,-XYT)F_{\tau}(XT,YT),\,F_{\tau}(u,v)=\frac{\theta^{{}^{\prime}}(0)\theta(u+v)}{\theta(u)\theta(v)}$

where
$\theta(u)=\sum_{n\in\mathbb{Z}}(-1)^{n}q^{\frac{1}{2}(n+\frac{1}{2})^{2}}e^{(n+\frac{1}{2})u}$
is a Jacobi theta function and $\mathcal{B}_{k}$ is a set of all Hecke
eigenforms of weight $k$ on $SL_{2}(\mathbb{Z}).$

The identity by Zagier (1) relates a generating function, which contains all
Hecke eigenforms together with all critical values, to the Jacobi form
$F_{\tau}(u,v).$ Such expansions with respect to the variable $T$ give an
algorithm to compute Hecke eigenforms (see [24] for more details).

It took almost 30 years to see that such an identity (1) is not accidental but
also exists for a general group $\Gamma_{0}(N)$ (see [8]). Now it seems
natural to ask if one can get such a relation for general automorphic forms.
In this paper we attempt to get such a formula, namely, an identity between a
generating function of periods of Hilbert modular forms over totally real
number fields with strict class number one and Jacobi forms (see Theorem 2.2).
The function $F_{\tau}(u,v)$ in (1) was introduced by Kronecker in a more
general form (see page 70 in [23]) and several properties of that have been
explored by Zagier [24]. The essential property of $F_{\tau}(u,v)$ is that it
can be identified as a sum of derivatives of Eisenstein series (called the
”Kuznetsov lifting”) and using this fact we are able to extend Zaiger’
identity to that for a totally real number field. The main result of this
paper shows the first connection between the Kronecker series and the critical
$L$\- values of Hilbert modular forms over the totally real number fields. It
also gives a systematic way to compute Hilbert Hecke eigenforms and the
special values of $L$ -functions by taking the expansions of the Kronecker
series.

This paper is organized as follows: in section $2$ we state (Main) Theorem 2.2
after introducing necessary notations. In section $3$ the analog of
Eisenstein-Kronecker series over a totally real number field and the
rationality of period polynomials of Hilbert modular forms are discussed.
Section $4$ gives detailed proof of Main Theorem. Finally, we give a comment
on a connection between parabolic cohomology and a period theory of Hilbert
modular forms. Also, we discuss a possible application of the Kronecker series
to evaluate the special $L$-values of a general automorphic form as a
conclusion.

Acknowledgement I would like to thank the referees for numerous helpful
comments and suggestions which greatly improved the exposition of this paper.

##  2\. Notations and statement of Main Theorem

###  2.1. Notation

  * •

$\mathbb{F}:$ a totally real number field of degree $t$ with discriminant $D$
and strict class number $1$

  * •

$\mathcal{O}:$ the ring of integers of $\mathbb{F}$ containing a unit of
negative norm

  * •

$\mathcal{O}^{*}:$ the group of units of $\mathcal{O}$

  * •

$\mathcal{O}^{*,+}:$ the group of totally positive units of $\mathcal{O}$

  * •

$U^{+}=\\{\epsilon^{2}\,:\,\epsilon\in\mathcal{O}^{*,+}\\}$

  * •

$\alpha\succ 0:\,$ a totally positive element

  * •

$\alpha_{1}\rightarrow\alpha_{2},\cdots,\alpha_{t}$ for the conjugation

  * •

$\mathcal{N}(\alpha)=\prod_{j=1}^{t}\alpha_{j}$ the norm

  * •

$tr(\alpha)=\sum_{j=1}^{t}\alpha_{j}$ the trace

  * •

$\mathfrak{D}:$ the different of $\mathbb{F}$

  * •

$\zeta_{\mathbb{F}}(s)=\sum_{{c}\subset\mathcal{O}}\frac{1}{\mathcal{N}({c})^{s}},$
where $s\in\mathbb{C}$ and the integral ideal $c$

  * •

$|\mathbf{r}|=\sum_{i=1}^{t}r_{i},\,\mathbf{r}+\mathbf{r}^{\prime}=(r_{1}+r_{1}^{\prime},\cdots,r_{t}+r_{t}^{\prime}),\,$
$\Gamma(\mathbf{r}+\mathbf{1})=\mathbf{r}!=r_{1}!\cdots
r_{t}!,\,\left(\begin{smallmatrix}\mathbf{r}\\\
\mathbf{r}^{\prime}\end{smallmatrix}\right)=\left(\begin{smallmatrix}r_{1}\\\
r_{1}^{\prime}\end{smallmatrix}\right)\left(\begin{smallmatrix}r_{2}\\\
r_{2}^{\prime}\end{smallmatrix}\right)\cdots\left(\begin{smallmatrix}r_{t}\\\
r_{t}^{\prime}\end{smallmatrix}\right)$ for
$\mathbf{r}=(r_{1},r_{2},\cdots,r_{t}),\mathbf{r}^{\prime}=(r_{1}^{\prime},r_{2}^{\prime},\cdots,r_{t}^{\prime})\in\mathbb{Z}_{\geq
0}^{t}$

  * •

$z^{\mathbf{r}}={z_{1}}^{r_{1}}{z_{2}}^{r_{2}}\cdots{z_{t}}^{r_{t}},\,tr(\mathbf{m}z)=\sum_{j=1}^{t}m_{j}z_{j},\,\mathcal{N}(z)=\prod_{j=1}^{t}z_{i}$
for $z=(z_{1},z_{2},\cdots,z_{t})\in\mathbb{C}^{t}$ and
${\mathbf{m}}\in\mathbb{F}$

  * •

$\mathbb{H}^{t}:$ the $t$-copies of complex upper half plane $\mathbb{H}$

  * •

$\tau=(\tau_{1},\cdots,\tau_{t})=x+\sqrt{-1}y\in\mathbb{H}^{t},x=(x_{1},\cdots,x_{t})\in\mathbb{R}^{t},y=(y_{1},\cdots,y_{t})\in(\mathbb{R}^{+})^{t},\,q=\prod_{j=1}^{t}q_{j},\,q_{j}=e^{2\pi
i\tau_{j}},\,1\leq j\leq t.$

  * •

$\sigma=(\sigma_{1},\cdots,\sigma_{t})\in\Gamma=SL_{2}(\mathcal{O})^{t}:$ an
element in the Hilbert modular group

  * •

The action of the group $\Gamma$, which is embedded into
$SL_{2}(\mathbb{R})\times\cdots\times SL_{2}(\mathbb{R}),$ on $\mathbb{H}^{t}$
is given by linear fractional transformations

$\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\tau=\frac{a\tau+b}{c\tau+d}=\bigl{(}\frac{a_{1}\tau_{1}+b_{1}}{c_{1}\tau_{1}+d_{1}},\cdots,\frac{a_{t}\tau_{t}+b_{t}}{c_{t}\tau_{t}+d_{t}}\bigr{)},\tau=(\tau_{1},\cdots,\tau_{t})\in\mathbb{H}^{t}$

  * •

For a holomorphic function $\chi$ on $\mathbb{H}^{t},$

$\chi^{(\ell)}(\tau)=\frac{\partial^{|\ell|}}{\partial{\tau}^{\ell}}\chi(\tau):=\frac{\partial^{|\ell|}}{\partial{\tau_{1}}^{\ell_{1}}\cdots\partial{\tau_{t}}^{\ell_{t}}}\chi(\tau),\,\,\,\forall\ell=(\ell_{1},\cdots,\ell_{t})\in\mathbb{Z}^{t}_{\geq
0}$

  * •

$\mathbb{D}^{\ell}\bigl{(}\chi(\tau)\bigr{)}:=\chi^{(\ell)}(\tau)$

  * •

$S_{\mathbf{k}}\subset M_{\mathbf{k}}:$ the space of Hilbert cusp form
$\subset$ the space of Hilbert modular form on $\Gamma$ with a parallel weight
$\mathbf{k}=(k,\cdots,k),$ even $k\geq 2.$

  * •

$\mathcal{B}^{0}_{k}\subset\mathcal{B}_{k}:$ a basis, consisting of all
normalized Hecke eigenforms, of $S_{\mathbf{k}}\subset M_{\mathbf{k}},$
respectively

  * •

$\mathbb{Q}_{f}:$ the field spanned by Fourier coefficients of $f$ over
$\mathbb{Q}$

  * •

$\bigl{<}f\,,\,g\bigr{>}:=\int_{\Gamma\backslash\mathbb{H}^{t}}f(\tau)\overline{g(\tau)}\frac{dx\,dy}{y^{2}},\,$
the Petersson inner product for $f\in S_{\mathbf{k}},g\in M_{\mathbf{k}}$

  * •

For a function $f$ on $\mathbb{H}^{t}$ and
$\mathbf{\ell}=(\ell_{1},\cdots,\ell_{t})\in\mathbb{Z}^{t},$

$\displaystyle(f|_{\ell}\sigma)(z):=(cz+d)^{-\mathbf{\ell}}f(\frac{az+b}{cz+d}),\sigma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\Gamma$

###  2.2. Statement of Main Theorem

Take a cusp form $f(\tau)=\sum_{\mathfrak{D}^{-1}\ni\nu\succ
0}a_{f}(\nu)e^{2\pi itr(\nu\tau)}$ in $S_{\mathbf{k}}$ and consider the
complete $L$-function of $f:$ for $s\in\mathbb{C},$

$\Lambda(f,s):=\int_{\mathbb{R}_{+}^{t}/U^{+}}f(iy)y^{s-1}\,dy=D^{s}(2\pi)^{-ts}\Gamma(s)^{t}L(f,s),$

where $\,L(f,s)=\sum_{\mathfrak{D}^{-1}/U^{+}\ni\nu\succ
0}\frac{a_{f}(\nu)}{\mathcal{N}(\nu)^{s}}\,\,\,(Re(s)\gg 0).$ It is well-known
that $\Lambda(f,s)$ has an analytic continuation and functional equation [4,
11]

$\displaystyle\Lambda(f,s)=(-1)^{\frac{tk}{2}}\Lambda(f,k-s).$

Consider the following polynomials in $X=(X_{1},\cdots,X_{t})$, called the
even (odd) period polynomial associated to $f:$

$\displaystyle R_{f}^{ev}(X)$ $\displaystyle:=$
$\displaystyle\sum_{\tiny{\begin{array}[]{cc}0\leq n\leq k-2\\\ n\equiv
0\pmod{2}\end{array}}}\frac{\Gamma({k}-{1})^{t}}{\Gamma(n+1)^{t}\Gamma(k-n-1)^{t}}R_{{k-2-n}}(f)\mathcal{N}(X)^{{n}},$
$\displaystyle R_{f}^{od}(X)$ $\displaystyle:=$
$\displaystyle\sum_{\tiny{\begin{array}[]{cc}0<n<k-2\\\ n\equiv
1\pmod{2}\end{array}}}\frac{(-1)^{nt}\Gamma({k}-{1})^{t}}{\Gamma(n+1)^{t}\Gamma(k-n-1)^{t}}R_{{k-2-n}}(f)\mathcal{N}(X)^{{n}},$
$\displaystyle
R_{f}(X):=(-1)^{t}\bigl{(}R_{f}^{ev}(X)+R_{f}^{od}(X)\bigr{)},\,\,$

where

$\displaystyle
R_{{n}}(f):=\int_{{\mathbb{R}_{+}}^{t}/U^{+}}f(\tau)\tau^{n}d\tau=i^{t(n+1)}\Lambda(f,n+1).$

Using the functional equation of $\Lambda(f,s)$ we get

$\displaystyle R_{{k-2-n}}(f)=(-1)^{t(n+1)}R_{{n}}(f)$

and so we get

$\mathcal{N}(X)^{k-2}R_{f}(-\frac{1}{X})=(-1)^{t}R_{f}(X).$

Let $f$ be a primitive (Hilbert) Hecke eigenform and consider the polynomial
of the $2t$-variables in $X=(X_{1},\cdots,X_{t})$ and $Y=(Y_{1},\cdots,Y_{t})$

$\displaystyle
R_{f}(X,Y):=(-1)^{t}\frac{R^{ev}_{f}(X)R^{od}_{f}(Y)+R^{ev}_{f}(Y)R^{od}_{f}(X)}{D^{k-\frac{1}{2}}\,(2i)^{t(k-3)}<f,\,f>}\in\mathbb{C}[X,Y].$

It transforms under $\sigma\in Gal(\mathbb{C}/\mathbb{Q})$ by
$R_{\sigma(f)}=\sigma(R_{f})$ so that $R_{f}(X,Y)$ has coefficients in the
number field $\mathbb{Q}_{f}$ generated by the Fourier coefficients of $f.$
Summing over the basis $\mathcal{B}_{k}^{0},$ consisting of all normalized
Hecke eigenforms of $S_{\mathbf{k}},$ the following function

(2.3) $\displaystyle
C_{k}^{cusp}(X,Y;\tau):=\sum_{f\in\mathcal{B}_{k}^{0}}R_{f}(X,Y)f(\tau)$

is in $\mathbb{Q}[[q]][X,Y]$ for each even integer $k\geq 2.$ Further we
extend the definition of $R_{f}(X,Y)$ (see section 3.2) to non-cusp forms and
include the Eisenstein series in the sum (2.3). Then we define

(2.4) $\displaystyle
C_{k}(X,Y;\tau):=\sum_{f\in\mathcal{B}_{k}}R_{f}(X,Y)f(\tau).$

######  Example 2.1.

Take $\mathbb{F}=\mathbb{Q}(\sqrt{5}).$

  1. (1)

$\displaystyle C_{2}(X,Y;\tau)$
$\displaystyle=\sum_{f\in\mathcal{B}_{k}}R_{f}(X,Y)f(\tau)=2^{4}\cdot 3\cdot
5\frac{\bigl{(}\mathcal{N}(X)+\mathcal{N}(Y)\bigr{)}\bigl{(}\mathcal{N}(XY)+1\bigr{)}}{\mathcal{N}(XY)}G_{\mathbb{F},2}(\tau),$

with the normalized Eisenstein series of weight $(k,k)$ on $\Gamma$ given by

$\displaystyle
G_{\mathbb{F},k}(\tau)=\frac{\zeta_{\mathbb{F}}(1-k)}{2^{2}}+\sum_{\mathcal{D}^{-1}\ni\nu\succ
0}\sigma_{k-1}(\nu\mathcal{D})e^{2\pi
itr(\nu\tau)},\sigma_{r}(\mathfrak{n})=\sum_{\mathfrak{c}|\mathfrak{n}}\mathcal{N}(\mathfrak{c})^{r}.$

  2. (2)

Let $\mathbb{F}=\mathbb{Q}(\sqrt{5})$ and take a unique cusp form $f$ of
weight $8$ on $\Gamma.$ Using the example in [1] we get

$\frac{R_{f}^{ev}(X)R_{f}^{od}(Y)}{5^{\frac{15}{2}}(2i)^{10}<f,f>}=c(1+\frac{361}{2^{2}\cdot
5}X^{2}+\frac{361}{2^{2}\cdot 5}X^{4}+X^{6})(Y+\frac{2}{3}Y^{3}+Y^{5})$

up to rational constant multiple $c.$

Combining all these functions into a single generating function to define

(2.5) $\displaystyle C(X,Y;\tau;T)$
$\displaystyle:=\frac{(\mathcal{N}(X)+\mathcal{N}(Y))(\mathcal{N}(XY)+(-1)^{t})}{\mathcal{N}(XYT)^{2}}+\sum_{k\geq
2}C_{k}(X,Y;\tau)\frac{\mathcal{N}(T)^{{k-2}}}{\Gamma(k-1)^{t}}.$

On the other hand, consider

(2.6) $\displaystyle\,\,F_{\tau}(u,v):=(-2)^{t}\sum_{k\geq
0}\widetilde{{G}_{\mathbb{F},{k}}}(\tau,\frac{uv}{2\pi
i})(\mathcal{N}(u)^{{k-1}}+\mathcal{N}(v)^{{k-1}}),u,v\in\mathbb{C}^{t},$

where

(2.9)
$\displaystyle\widetilde{{G}_{\mathbb{F},{k}}}(\tau,\lambda):=\biggl{\\{}\begin{array}[]{cc}\sum_{\ell=(\ell_{1},\cdots,\ell_{t})\in\mathbb{Z}^{t}_{\geq
0}}\frac{{\lambda}^{\ell}}{{\ell}!({\ell}+\mathbf{{k}-{1}})!}\mathbb{D}^{\ell}\bigl{(}{G}_{\mathbb{F},{k}}(\tau)\bigr{)}&\mbox{
if $k\geq 2$}\\\ \frac{1}{2^{t}}&\mbox{if
$k=0$}\end{array}\biggr{\\}},\lambda\in\mathbb{C}^{t}$

and a normalized Hilbert Eisenstein series $G_{\mathbb{F},{k}}(\tau)$ (p 20 in
[11]) defined by

$\displaystyle E_{\mathbb{F},{k}}(\tau)$ $\displaystyle:=$
$\displaystyle\frac{D^{\frac{1}{2}-k}(2\pi
i)^{tk}}{\Gamma(k)^{t}}\bigl{(}\frac{1}{2^{t}}\zeta_{\mathbb{F}}(1-k)+\sum_{\tiny{\begin{array}[]{cc}\nu\in\mathfrak{D}^{-1}\\\
\nu\succ 0\end{array}}}\sigma_{k-1}(\nu\mathfrak{D})e^{2\pi
itr(\nu\tau)}\bigr{)}$ $\displaystyle:=$
$\displaystyle\frac{D^{\frac{1}{2}-k}(2\pi
i)^{tk}}{\Gamma(k)^{t}}G_{\mathbb{F},{k}}(\tau),\,\sigma_{r}(\mathfrak{n})=\sum_{\mathfrak{c}|\mathfrak{n}}\mathcal{N}(\mathfrak{c})^{r}.$

Now we state the main result of this paper :

######  Theorem 2.2.

(Main Theorem) Let $C(X,Y;\tau;T)$ be the generating function of the periods
of Hilbert modular forms given in (2.5). Then we have

  1. (1)

$C(X,Y;\tau;T)\in\frac{1}{\mathcal{N}(XYT)^{2}}{\mathbb{Q}}[X,Y][[q,T]].$

  2. (2)

$C(X,Y;\tau;T)=F_{\tau}(T,-XYT)\,F_{\tau}(XT,YT).$

######  Remark 2.1.

$\widetilde{{G}_{\mathbb{F},{k}}}(\tau,\lambda)$ in (2.9) is the ”Kuznetsov
lifting” of the Hilbert Eisenstein series $G_{\mathbb{F},\mathbf{k}}(\tau).$
Its modular transformation property is known (see Theorem 2 in [6]) : for any
$\left(\begin{smallmatrix}a&b\\\ c&d\end{smallmatrix}\right)\in\Gamma,k\geq
2,$

$\mathcal{N}(c\tau+d)^{-{k}}e^{-tr(\frac{c\lambda}{c\tau+d})}\widetilde{{G}_{\mathbb{F},{k}}}(\frac{a\tau+b}{c\tau+d},\frac{\lambda}{c\tau+d})=\widetilde{{G}_{\mathbb{F},{k}}}(\tau,\lambda)$

and its generating function $\mathcal{F}_{\tau}(u,v):=(-2)^{t}\sum_{k\geq
2}\widetilde{{G}_{\mathbb{F},{k}}}(\tau,\frac{uv}{2\pi
i})(\mathcal{N}(u)^{{k-1}}+\mathcal{N}(v)^{{k-1}})$ behaves as a Jacobi-like
form [10, 7] with a modular transformation property

$\mathcal{F}_{\frac{a\tau+b}{c\tau+d}}(\frac{u}{c\tau+d},\frac{v}{c\tau+d})=\mathcal{N}(c\tau+d)\,e^{tr(\frac{cuv}{c\tau+d})}{\mathcal{F}}_{\tau}(u,v),\forall\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\Gamma.$

##  3\.  Algebraicity and Period of Hecke eigen forms

###  3.1. Algebraicity

The study of period relation for automorphic forms was started by Shimura. He
showed the existence of relations up to factors in $\overline{\mathbb{Q}}^{*}$
in many instances and made a general conjecture relating periods of Hilbert
modular varieties and their compact analogs, that is, the quaternionic modular
varieties [17, 18]. There is a weaker conjecture, which gives a relation
between a product of two periods, called the quadratic periods, may be
interpreted, up to algebraic factors, as Petersson inner products. This was
proved by M. Harris [12] under a certain technical condition. More precisely,
for each $m,0\leq m\leq k-2,$ $\Lambda(f,m+1)$ is called the critical values.

######  Theorem 3.1.

(Theorem 4.3 in [19]) Let $f$ be a Hilbert Hecke eigenform of weight
$\mathbf{k}=(k,\cdots,k)$ over a totally real number field $\mathbb{F}$ of
degree $t$ and $\sigma\in Gal(\overline{\mathbb{Q}}/\mathbb{Q}).$

  1. (1)

For each $r\in\mathbb{Z}^{t}/2\mathbb{Z}^{t}$ and for $f^{\sigma},\sigma\in
Gal(\overline{\mathbb{Q}}/\mathbb{Q}),$ there exist nonzero complex numbers
$\omega_{f}^{r}$ such that $(\frac{L(f,m)}{(2\pi
i)^{tm}\omega_{f}^{r}})^{\sigma}=\frac{L(f^{\sigma},m)}{(2\pi
i)^{tm}\omega_{f^{\sigma}}^{r}},$ for any integer $m$ such that $0<m<k.$

  2. (2)

$\frac{L(f,m)}{(2\pi i)^{tm}\omega_{f}^{r}}\in\mathbb{Q}_{f}.$

  3. (3)

If $p=(p_{1},\cdots,p_{t}),r=(r_{1},\cdots,r_{t})$ with $p_{i}+r_{i}\equiv
1\pmod{2},1\leq\forall i\leq t,$ we have $\frac{w_{f}^{p}\cdot
w_{f}^{r}}{<f\,,f>}\in\mathbb{Q}_{f}$ and $\bigl{(}\frac{w_{f}^{p}\cdot
w_{f}^{r}}{<f\,,f>}\bigr{)}^{\sigma}=\frac{w_{f^{\sigma}}^{p}\cdot
w_{f^{\sigma}}^{r}}{<f^{\sigma}\,,f^{\sigma}>}.$

###  3.2.  Period of non-cusp forms

Since $\mathcal{B}_{k}$ in (2.4) contains non-cusp forms one needs to explain
” period function” corresponding to a non-cusp form $f.$ Take a non-cusp form
$f(\tau)=\sum_{0\preceq\nu\in\mathfrak{D}^{-1}}a_{f}(\nu)e^{2\pi
itr(\nu\tau)}$ in $M_{\mathbf{k}}$ and consider

$\Lambda(f,s):=\int_{(\mathbb{R}_{+})^{t}/U^{+}}\bigl{(}f(iy)-a_{f}(0)\bigr{)}y^{{s}-1}dy=D^{s}(2\pi)^{-ts}\Gamma(s)^{t}L(f,s),s\in\mathbb{C}.$

It has a meromorphic continuation to $\mathbb{C}$ and satisfies a functional
equation $\Lambda(f,s)=(-1)^{\frac{tk}{2}}\Lambda(f,k-s),$ but now has simple
poles of residue as $-a_{f}(0)$ and $(-1)^{kt}a_{f}(0)$ , up to a constant
multiple, at $s=0$ and $s=k,$ respectively.

Define

$\displaystyle R_{f}(X):=\frac{(-1)^{t}\sqrt{D}\cdot
a_{f}(0)}{({k-1})^{t}}(\mathcal{N}(X)^{{k-1}}+(-1)^{t}\mathcal{N}(X)^{-{1}})$
$\displaystyle+$
$\displaystyle\sum_{n=0}^{k-2}(-1)^{\frac{t(k+n-1)}{2}}\frac{\Gamma({k-1})^{t}}{\Gamma(n+1)^{t}\Gamma(k-n-1)^{t}}\Lambda(f,k-1-n)\mathcal{N}(X)^{n}.$

The assumption that $\mathbb{F}$ has the strict class number $1$ implies that
the space of Hilbert modular forms is a direct sum (see [4], p 12)

$M_{\mathbf{k}}=S_{\mathbf{k}}\oplus<G_{\mathbb{F},k}>$

and, so (2.4) becomes

$C_{k}(X,Y;\tau)=\sum_{f\in\mathcal{B}_{k}^{0}}R_{f}(X,Y)f(\tau)+R_{G_{\mathbb{F},k}}(X,Y)G_{\mathbb{F},k}(\tau),$

where $R_{G_{\mathbb{F},k}}(X,Y),$ defined by

(3.2) $\displaystyle
R_{G_{\mathbb{F},k}}(X,Y)=(-1)^{t}\frac{R_{G_{\mathbb{F},k}}^{ev}(X)R_{G_{\mathbb{F},k}}^{od}(Y)+R_{G_{\mathbb{F},k}}^{ev}(Y)R_{G_{\mathbb{F},k}}^{od}(X)}{D^{k-\frac{1}{2}}(2i)^{t(k-3)}<G_{\mathbb{F},k},G_{\mathbb{F},k}>},$

is a symmetrized sum of the product of period polynomials of the normalized
Hecke Eisenstein series $G_{\mathbb{F},{k}}(\tau)$ given as followings:

######  Proposition 3.2.

Take

$w_{G_{\mathbb{F},{k}}}^{-}=\frac{\sqrt{D}\Gamma(k-1)^{t}}{2^{t}},\,\,w_{G_{\mathbb{F},{k}}}^{+}=\frac{D^{k-\frac{3}{2}}\zeta_{\mathbb{F}}(k-1)}{(2\pi
i)^{t(k-1)}}w_{G_{\mathbb{F},{k}}}^{-}$

and

$\displaystyle p_{{k}}^{+}(X)=\mathcal{N}(X)^{k-2}+(-1)^{t},\,$ $\displaystyle
p_{{k}}^{-}(X)=\sum_{-1\leq n\leq k-1,n\equiv
1\pmod{2}}\frac{\zeta_{\mathbb{F}}(1-(n+1))\zeta_{\mathbb{F}}(n+2-k)}{\Gamma(n+1)^{t}\Gamma(k-n-1)^{t}}{\mathcal{N}(X)}^{n}.$

For $k\geq 2$ the period function of $G_{\mathbb{F},{k}}(\tau)$ is given by

(3.3) $\displaystyle
R_{G_{\mathbb{F},{k}}}(X)=(-1)^{t}\bigl{(}w_{G_{\mathbb{F},{k}}}^{-}\cdot
p_{{k}}^{-}(X)+w_{G_{\mathbb{F},{k}}}^{+}\cdot p_{{k}}^{+}(X)\bigr{)}$

so that

$R_{G_{\mathbb{F},{k}}}^{ev}(X)=w_{G_{\mathbb{F},{k}}}^{+}\cdot
p_{{k}}^{+}(X)\mbox{ and
}R_{G_{\mathbb{F},{k}}}^{od}(X)=w_{G_{\mathbb{F},{k}}}^{-}\cdot
p_{{k}}^{-}(X).$

######  Remark 3.3.

  1. (1)

Like in the case of an elliptic modular form (see [14]), the period function
$R_{G_{\mathbb{F},{k}}}(X)$ is $\frac{1}{\mathcal{N}(X)}$ times a polynomial :

$R_{G_{\mathbb{F},{k}}}(X)\in\frac{1}{\mathcal{N}(X)}\mathbb{C}[X].$

  2. (2)

[11] Note that
$\frac{D^{n-\frac{1}{2}}\zeta_{\mathbb{F}}(n)\Gamma(n)^{t}}{(2\pi
i)^{tn}}=\frac{\zeta_{\mathbb{F}}(1-n)}{2^{t}}$ and $\zeta_{\mathbb{F}}(-n)=0$
for any positive even integer $n.$

  3. (3)

    1. (a)

[26] Let $\mathbb{F}$ be a real quadratic field with discriminant $D.$ It is
known that $\zeta_{\mathbb{F}}(1-n)=B_{n}B_{n,\chi}$ for even positive integer
$n.$ $B_{r}$ and $B_{r,\chi}$ are the $r$th Bernoulli number
$(B_{0}=1,B_{1}=-\frac{1}{2},B_{2}=\frac{1}{6},\cdots)$ and the $r$th twisted
Bernoulli number
$(B_{1,\chi}=\frac{1}{D}\sum_{a=1}^{D}\chi(a)a,B_{2,\chi}=\frac{1}{D}\sum_{a=1}^{D}\chi(a)a^{2}-\sum_{a=1}^{D}\chi(a)a,B_{3,\chi}=\cdots),$
respectively. Here $\chi\pmod{D}$ is a primitive character defined by
$\chi(\cdot)=\bigl{(}\frac{D}{\cdot}\bigr{)}.$

    2. (b)

(open problem)  It is well known that the generating functions of $B_{j}$ and
$B_{j,\chi}$ are

$\sum_{n=0}^{\infty}B_{n}\frac{t^{n}}{n!}=\frac{te^{t}}{e^{t}-1}\mbox{\, and
}\sum_{m\geq
0}B_{m,\chi}\frac{t^{m}}{m!}=\sum_{j=1}^{D}\frac{\chi(j)te^{jt}}{e^{Dt}-1}.$

Similarly, it will be interesting to express the following generating function

$\sum_{m\geq 2}^{\infty}B_{m}B_{m,\chi}\frac{t^{m}}{m!}=\sum_{m\geq
2}^{\infty}\zeta_{\mathbb{F}}(1-m)\frac{t^{m}}{m!},$

as elementary functions.

Proof  of Proposition 3.2 : The period polynomial of $G_{\mathbb{F},k}(\tau)$
can be computed from the definition in (3.2) : using
$\Lambda(G_{\mathbb{F},k},n+1)=\frac{D^{n+1}\Gamma(n+1)^{t}}{(2\pi)^{t(n+1)}}\zeta_{\mathbb{F}}(n+1)\zeta_{\mathbb{F}}(n+2-k),$
we have

$\displaystyle
R_{G_{\mathbb{F},k}}(X)=\frac{(-1)^{t}\sqrt{D}\zeta_{\mathbb{F}}(1-k)}{2^{t}(k-1)^{t}}(\mathcal{N}(X)^{k-1}+(-1)^{t}\mathcal{N}(X)^{-1})$
$\displaystyle+\frac{(-1)^{t}{D}^{k-1}\Gamma(k-1)^{t}\zeta_{\mathbb{F}}(k-1)}{2^{t}(2\pi
i)^{t(k-1)}}\bigl{(}\mathcal{N}(X)^{k-2}+(-1)^{t}\bigr{)}$
$\displaystyle+\sum_{n=1}^{k-3}(-1)^{\frac{t(k+n-1)}{2}}\frac{\Gamma({k-1})^{t}D^{k-n-1}}{(2\pi)^{t(k-n-1)}\Gamma(n+1)^{t}}\zeta_{\mathbb{F}}(1-(n+1))\zeta_{\mathbb{F}}(k-1-n)\mathcal{N}(X)^{{n}}$

(using the functional equation of $\zeta_{\mathbb{F}}(s)$ of Remark 3.3-(2))

$\displaystyle=\frac{(-1)^{t}{D}^{k-1}\Gamma(k-1)^{t}\zeta_{\mathbb{F}}(k-1)}{2^{t}(2\pi
i)^{t(k-1)}}\bigl{(}\mathcal{N}(X)^{k-2}+(-1)^{t}\bigr{)}$
$\displaystyle+\frac{(-1)^{t}\sqrt{D}\Gamma(k-1)^{t}}{2^{t}}\sum_{-1\leq
n\equiv 1\pmod{2}\leq
k-1}\frac{\zeta_{\mathbb{F}}(1-(n+1))\zeta_{\mathbb{F}}(n+2-k)}{\Gamma(n+1)^{t}\Gamma(k-n-1)^{t}}\mathcal{N}(X)^{n}$
$\displaystyle=(-1)^{t}\bigl{(}\omega_{G_{\mathbb{F},k}}^{+}p_{{k}}^{+}(X)+\omega_{G_{\mathbb{F},k}}^{-}p_{{k}}^{-}(X)\bigr{)}.$

This completes a proof. ∎

####  3.2.1. Petersson scalar product of Eisenstein Series

The Petersson scalar product of a $SL_{2}(\mathbb{Z})$-invariant function had
been defined by Zagier [25] using Rankin-Selberg method. Similarly we have

######  Proposition 3.4.

$<G_{\mathbb{F},{k}}(\tau),G_{\mathbb{F},{k}}(\tau)>=\frac{\Gamma(k-1)^{t}\zeta_{\mathbb{F}}(k-1)}{(4\pi)^{t(k-1)}}\frac{\zeta_{\mathbb{F}}(1-k)}{2^{t}}$

Proof  of Proposition 3.4 : Following the method in [25] the Petersson norm of
the Hilbert Eisenstein series $G_{\mathbb{F},{k}}(\tau)$ can be computed as

$\displaystyle<G_{\mathbb{F},{k}}(\tau),G_{\mathbb{F},{k}}(\tau)>=(-1)^{\frac{tk}{2}}(4\pi)^{-tk}\Gamma(k)^{t}\cdot\zeta^{*}_{\mathbb{F}}(k)\zeta^{*}_{\mathbb{F}}(2-k)$

where
$\zeta_{\mathbb{F}}^{*}(s):=D^{\frac{s}{2}}\pi^{-\frac{ts}{2}}\Gamma(\frac{s}{2})^{t}\zeta_{\mathbb{F}}(s)=\zeta_{\mathbb{F}}^{*}(1-s)$
(see p 57 [26]). Using the identities
$\Gamma(\frac{k}{2})\Gamma(\frac{k-1}{2})=\Gamma(k-1)\sqrt{\pi}2^{-(k-2)}$ and
$\zeta_{\mathbb{F}}(k)=D^{-k+\frac{1}{2}}\frac{(2\pi
i)^{tk}}{2^{t}\Gamma(k)^{t}}\zeta_{\mathbb{F}}(1-k)$ we get the Petersson norm
of the Eisenstein series $G_{\mathbb{F},k}.$ ∎

Now write the function $C_{k}(X,Y;\tau)$ in (2.4) as a sum of
$C^{cusp}_{k}(X,Y;\tau)$ in (2.3) and
$C^{Eis}_{k}(X,Y;\tau):=R_{G_{\mathbb{F},k}}(X,Y)G_{\mathbb{F},k}(\tau):$

$C_{k}(X,Y;\tau)=C^{cusp}_{k}(X,Y;\tau)+C^{Eis}_{k}(X,Y;\tau).$

######  Proposition 3.5.

We have

  1. (1)

$C^{Eis}_{k}(X,Y;\tau)=(-1)^{t}\frac{2^{t}\Gamma(k-1)^{t}}{\zeta_{\mathbb{F}}(1-k)}(p^{+}_{k}(X)p_{k}^{-}(Y)+p^{+}_{k}(Y)p_{k}^{-}(X))G_{\mathbb{F},k}(\tau).$

  2. (2)

$C(X,Y;\tau;T)=F_{\tau}(XT,YT)F_{\tau}(T,-XYT)\mbox{ as
$\tau\rightarrow(i\infty,\cdots,i\infty)$}$

Proof  of Proposition 3.5 :

  1. (1)

From (3.2) recall that

$R_{G_{\mathbb{F},k}}(X,Y)=(-1)^{t}\frac{R_{G_{\mathbb{F},k}}^{ev}(X)R_{G_{\mathbb{F},k}}^{od}(Y)+R_{G_{\mathbb{F},k}}^{ev}(Y)R_{G_{\mathbb{F},k}}^{od}(X)}{D^{k-\frac{1}{2}}(2i)^{t(k-3)}<G_{\mathbb{F},k},G_{\mathbb{F},k}>}.$

So, Proposition 3.2 and Proposition 3.4 imply that

$\displaystyle
C^{Eis}_{k}(X,Y;\tau)=R_{G_{\mathbb{F},k}}(X,Y)G_{\mathbb{F},k}(\tau)$
$\displaystyle=\frac{\omega^{+}_{G_{\mathbb{F},k}}\omega^{-}_{G_{\mathbb{F},k}}\bigl{(}p^{+}_{k}(X)p^{-}_{k}(Y)+p^{+}_{k}(Y)p^{-}_{k}(X)\bigr{)}}{D^{k-\frac{1}{2}}(2i)^{t(k-3)}<G_{\mathbb{F},k},\,G_{\mathbb{F},k}>}G_{\mathbb{F},k}(\tau)$
$\displaystyle=(-1)^{t}\frac{2^{t}\Gamma(k-1)^{t}}{\zeta_{\mathbb{F}}(1-k)}(p^{+}_{k}(X)p_{k}^{-}(Y)+p^{+}_{k}(Y)p_{k}^{-}(X))G_{\mathbb{F},k}(\tau).$

  2. (2)

Using Proposition 3.5 part (1) the value of $C(X,Y;\tau;T)$ as
$\tau\rightarrow(i\infty,\cdots,i\infty)$ is

$\displaystyle C(X,Y;(i\infty,\cdots,i\infty);T)$
$\displaystyle=\frac{(\mathcal{N}(X)+\mathcal{N}(Y))(\mathcal{N}(XY)+(-1)^{t})}{\mathcal{N}(XYT)^{2}}+\sum_{k\geq
2}C_{k}^{Eis}(X,Y;(i\infty,\cdots,i\infty))\frac{\mathcal{N}(T)^{k-2}}{\Gamma(k-1)^{t}}$
$\displaystyle=$
$\displaystyle\frac{(\mathcal{N}(X)+\mathcal{N}(Y))(\mathcal{N}(XY)+(-1)^{t})}{\mathcal{N}(XYT)^{2}}+(-1)^{t}\sum_{k\geq
2}\bigl{(}p_{{k}}^{+}(X)p_{{k}}^{-}(Y)+p_{{k}}^{+}(Y)p_{{k}}^{-}(X)\bigr{)}{\mathcal{N}(T)^{k-2}}$

since $G_{\mathbb{F},k}(i\infty)=\frac{\zeta_{\mathbb{F}}(1-k)}{2^{t}}.$ On
the other hand, a direct computation shows that

$\displaystyle
F_{\tau}(T,-XYT)F_{\tau}(XT,YT)|_{\tau\rightarrow(i\infty,\cdots,i\infty)}$
$\displaystyle=\frac{(\mathcal{N}(X)+\mathcal{N}(Y))(\mathcal{N}(XY)+(-1)^{t})}{\mathcal{N}(XYT)^{2}}+(-1)^{t}\sum_{k\geq
2}\bigl{(}p^{+}_{k}(X)p^{-}_{k}(Y)+p^{+}_{k}(Y)p^{-}_{k}(X)\bigr{)}\mathcal{N}(T)^{k-2}$

This completes the proof of Proposition 3.5.

∎

##  4\. Proofs

###  4.1. Proof  of  Theorem 2.2 (Main Theorem)

$(1)$ Using Theorem 3.1 with a proper choice of the Petersson norm $<f,\,f>,$
we see that

$\displaystyle
R_{f}(X,Y):=(-1)^{t}\frac{R^{ev}_{f}(X)R^{od}_{f}(Y)+R^{ev}_{f}(Y)R^{od}_{f}(X)}{D^{k-\frac{1}{2}}\,(2i)^{t(k-3)}<f,\,f>}\in\mathbb{Q}_{f}[X,Y],f\in
S_{\mathbf{k}}.$

With an action of $\sigma\in Gal(\mathbb{C}/\mathbb{Q}_{f})$ by
$R_{\sigma(f)}=\sigma(R_{f})$ we see that

$\displaystyle
C(X,Y;\tau;T)=\frac{(\mathcal{N}(X)+\mathcal{N}(Y))(\mathcal{N}(XY)+(-1)^{t})}{\mathcal{N}(XYT)^{2}}$
$\displaystyle+$ $\displaystyle\sum_{k\geq
2}\sum_{f\in\mathcal{B}_{k}}R_{f}(X,Y)f(\tau)\frac{\mathcal{N}(T)^{{k-2}}}{\Gamma(k-1)^{2}}\in\frac{1}{\mathcal{N}(XYT)^{2}}\mathbb{Q}[X,Y][[q,T]].$

This proves rationality of $C(X,Y;\tau;T).$

$(2)$ To prove Theorem 2.2 part (2) write the Taylor expansion (2.6)

$\displaystyle F_{\tau}(u,v)$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}(u)}+\frac{1}{\mathcal{N}(v)}$
$\displaystyle+$ $\displaystyle(-2)^{t}\sum_{k\geq
2}\sum_{\ell\in\mathbb{Z}_{\geq
0}^{t}}\frac{\mathbb{D}^{\ell}\bigl{(}G_{\mathbb{F},k}(\tau)\bigr{)}}{(2\pi
i)^{\ell}\ell!(\ell+\mathbf{k}-\mathbf{1})!}(u^{\ell}v^{\ell+\mathbf{k}-1}+u^{\ell+\mathbf{k}-1}v^{\ell})$

or write it as

(4.1) $\displaystyle
F_{\tau}(u,v)=\sum_{\mathbf{h}=(h,\cdots,h),\mathbf{\ell}=(\ell_{1},\cdots,\ell_{2})}g_{{h},\mathbf{\ell}}(\tau)(u^{\ell}v^{\ell+\mathbf{h}-\mathbf{1}}+u^{\ell+\mathbf{h}-\mathbf{1}}v^{\ell})$

with

$g_{{h},\mathbf{\ell}}(\tau)=\left\\{\begin{array}[]{ccll}\frac{(-2)^{t}}{\,(2\pi
i)^{\ell}\Gamma(\mathbf{\ell}+\mathbf{1})\Gamma(\mathbf{\ell}+\mathbf{h})}\mathbb{D}^{\ell}(G_{\mathbb{F},{h}}(\tau)),&\mbox{\,
if $h\geq 2,\ell_{i}\geq 0,i=1,2$}\\\ \frac{1}{2^{t}},&\mbox{\, if
${h}=0,\ell=(0,\cdots,0)$}\\\ 0,&\mbox{\, otherwise}\end{array}\right\\}.$

Next let

$F_{\tau}(T,-XYT)F_{\tau}(XT,YT)=\sum_{{k}\geq
0}b_{\mathbf{k}}(X,Y;\tau)\frac{\mathcal{N}(T)^{{k}-{2}}}{\Gamma(k-1)^{2}}.$

Since we have already checked that

$C(X,Y;(i\infty,\cdots,i\infty),T)=F_{\tau}(T,-XYT)F_{\tau}(XT,YT)|_{\tau\rightarrow(i\infty,\cdots,i\infty)}$

in Proposition 3.5, it is enough to confirm that

$\frac{\bigl{<}b_{\mathbf{k}}(X,Y;\cdot),f(\cdot)\bigr{>}}{\bigl{<}f,\,f\bigr{>}}=R_{f}(X,Y)\mbox{\,\,
for each $f\in\mathcal{B}^{0}_{k}$}.$

From the expression (4.1) we see

(4.6) $\displaystyle b_{\mathbf{k}}(X,Y;\tau)$ $\displaystyle=$
$\displaystyle\sum_{\tiny{\begin{array}[]{ccc}\ell,\mathbf{h},\ell^{\prime},\mathbf{h}^{\prime}\\\
\mathbf{h}+\mathbf{h}^{\prime}+{2}(\mathbf{\ell}+\mathbf{\ell}^{\prime})=\mathbf{k}\\\
\mathbf{h}=(h,\cdots,h),\mathbf{h}^{\prime}=(h^{\prime},\cdots,h^{\prime})\end{array}}}g_{{h},\mathbf{\ell}}(\tau)g_{{h}^{\prime},\mathbf{\ell}^{\prime}}(\tau)$
$\displaystyle\times[(-XY)^{\mathbf{\ell}+\mathbf{h}-\mathbf{1}}+(-XY)^{\mathbf{\ell}}][X^{\mathbf{\ell}^{\prime}}Y^{\mathbf{\ell}^{\prime}+\mathbf{h}^{\prime}-\mathbf{1}}+X^{\mathbf{\ell}^{\prime}+\mathbf{h}^{\prime}-\mathbf{1}}Y^{\mathbf{\ell}^{\prime}}]$

The coefficients of $\mathcal{N}(X)^{{p}}\mathcal{N}(Y)^{{q}}$ with $q$ or $p$
equal to $-1$ or to $k-1$ involve only $G_{\mathbb{F},k}(\tau)$ and have
already been treated in Proposition 3.5. Also the coefficients of
$\mathcal{N}(X)^{p}\mathcal{N}(Y)^{q}$ in (4.6) is invariant under
$q\leftrightarrow k-2-p$ and $q\leftrightarrow p$ so that we may assume $0\leq
p<q\leq\frac{k-2}{2}.$ For such $p,q$ the coefficient of
$\mathcal{N}(X)^{{p}}\mathcal{N}(Y)^{{q}}$ in (4.6) equals

$\displaystyle\sum_{\tiny{\begin{array}[]{cc}\ell,\ell^{\prime}\succeq-\mathbf{1}\\\
\ell+\ell^{\prime}=\mathbf{p}=(p,\cdots,p)\end{array}}}(-1)^{|\ell|}g_{{k}-{p}-{q}-{1},\mathbf{\ell}}(\tau)g_{{q}-{p}+{1},\mathbf{\ell}^{\prime}}(\tau)$
$\displaystyle=$
$\displaystyle\sum_{\tiny{\begin{array}[]{cc}\ell,\ell^{\prime}\succeq-\mathbf{1}\\\
\ell+\ell^{\prime}=\mathbf{p}=(p,\cdots,p)\end{array}}}2^{2t}\frac{(-1)^{|\ell|}\mathbb{D}^{\ell}(G_{\mathbb{F},k-p-q-1})\mathbb{D}^{\ell^{\prime}}(G_{\mathbb{F},q-p+1})}{\ell!(\ell+\mathbf{k-p-q}-2)!{\ell^{\prime}}!(\mathbf{q}-\ell^{\prime})!}$
$\displaystyle=$ $\displaystyle\frac{2^{2t}}{\Gamma({q}+\
{1})^{t}\Gamma({k}-{q}-{1})^{t}}[G_{\mathbb{F},{q}-{p}+{1}},G_{\mathbb{F},{k}-{q}-{p}-{1}}]^{Hil}_{\mathbf{p}}.$

Here, $[\cdot\,\,\cdot]^{Hil}_{\mathbf{p}}$ denotes the
$\mathbf{p}=(p,\cdots,p)$th Rankin-Cohen bracket (see Corollary 1 in [6])
defined by

$\displaystyle[G_{\mathbb{F},{q}-{p}+{1}},G_{\mathbb{F},{k}-{q}-{p}-{1}}]^{Hil}_{\mathbf{p}}:=$
$\displaystyle\frac{1}{(2\pi
i)^{tp}}\sum_{\tiny{\begin{array}[]{ccc}0\leq{\ell}_{i}\leq{p}\\\
\ell=(\ell_{1},\cdots,\ell_{t})\\\
\ell+\ell^{\prime}=\mathbf{p}\end{array}}}\frac{(-1)^{|\mathbf{\ell}|}\Gamma({q}+{1})^{t}\Gamma({k}-{q}-{1})^{t}}{{\ell^{\prime}}!(\mathbf{{k}+{\ell^{\prime}}-{q}-{p-2}})!{\ell}!(\mathbf{{q}-{\ell}})!}\mathbb{D}^{\mathbf{\ell}}(G_{\mathbb{F},{q}-{p}+{1}})\mathbb{D}^{{\ell^{\prime}}}(G_{\mathbb{F},{k}-{q}-{p}-{1}}).$

On the other hand we recall the following result (Theorem 3 in [6]) :

######  Theorem 4.1.

[6] Suppose that $f(\tau)=\sum_{\mathcal{D}^{-1}\ni\nu\succ
0}a_{f}(\nu)e^{2\pi itr(\nu\tau)}\in S_{\mathbf{k}}$ and
$g(\tau)=\sum_{\mathcal{D}^{-1}\ni\nu\succeq 0}a_{g}(\nu)e^{2\pi
itr(\nu\tau)}\in M_{\mathbf{k}_{2}}$ with $k=k_{1}+k_{2}+2p>2.$ Then

$\frac{D^{\frac{1}{2}-k_{1}}(2\pi
i)^{tk_{1}}}{\Gamma(k_{1})^{t}}<f,[G_{\mathbb{F},k_{1}},g_{k_{2}}]_{\mathbf{p}}>=\frac{\Gamma(k-1)^{t}\Gamma(k_{1}+p)^{t}}{(4\pi)^{t(k-1)}\Gamma(k_{1})^{t}\Gamma(p+1)^{t}}\sum_{\nu\succ
0}\frac{a_{f}(\nu)\overline{a_{g}(\nu)}}{\mathcal{N}(\nu)^{k-p-1}}$

Taking $g_{k_{2}}(\tau)=G_{\mathbb{F},k}(\tau)$ in Theorem 4.1 we get

$\displaystyle<f,[G_{\mathbb{F},k_{1}},G_{\mathbb{F},k_{2}}]_{\mathbf{p}}>=(-1)^{\frac{tk_{1}}{2}}\frac{D^{k_{1}-\frac{1}{2}}\Gamma(k-1)^{2}\Gamma(k_{1}+p)^{t}}{2^{t(k-1)}(2\pi)^{t(k+k_{1}-1)}\Gamma(p+1)^{t}}\sum_{\nu>0}\frac{a_{f}(\nu)\sigma_{k_{2}-1}(\nu)}{\mathcal{N}(\nu)^{k-p-1}}$
$\displaystyle=(-1)^{\frac{tk_{1}}{2}}\frac{D^{k_{1}-\frac{1}{2}}\Gamma(k-1)^{t}\Gamma(k_{1}+p)^{t}}{2^{t(k-1)}(2\pi)^{t(k+k_{1}-1)}\Gamma(p+1)^{t}}L(f,k-p-1)L(f,k_{1}+p)$
$\displaystyle(\mbox{since
$R_{n}(f)=i^{t(n+1)}D^{n+1}(2\pi)^{-t(n+1)}\Gamma(n+1)^{t}L(f,n+1)$})$
$\displaystyle=\frac{(-1)^{\frac{t(k-1)}{2}}\Gamma(k-1)^{t}}{D^{k-\frac{1}{2}}2^{t(k-1)}\Gamma(p+1)^{t}\Gamma(k-p-1)^{t}}R_{k-p-2}(f)R_{k_{1}+p-1}(f).$

And so we have that

$\displaystyle\sum_{f\in\mathcal{B}^{0}_{k}}\frac{R_{k-p-2}(f)R_{k_{1}+p-1}(f)f(\tau)}{i^{t(k-1)}D^{k-\frac{1}{2}}2^{t(k-3)}<f,f>}=\frac{2^{2t}\Gamma(p+1)^{t}\Gamma(k-p-1)^{t}}{\Gamma(k-1)^{t}}[G_{\mathbb{F},k_{1}},G_{\mathbb{F},k_{2}}]^{Hil}_{\mathbf{p}}.$

Now take $k_{1}=k-q-p-1,k_{2}=q+1-p$ for $q+p\equiv 1\pmod{2},p,q>0,$ to get

$\displaystyle\sum_{f\in\mathcal{B}^{0}_{k}}(-1)^{t}\frac{R_{{k-q-2}}(f)R_{{k-p-2}}(f)}{D^{k-\frac{1}{2}}(2i)^{t(k-3)}<f,\,f>}f(\tau)=\frac{2^{2t}\Gamma({k-p-1})^{t}\Gamma({p+1})^{t}}{\Gamma({k-1})^{t}}[G_{\mathbb{F},{k-1-q-p}},\,G_{\mathbb{F},{q+1-p}}]^{Hil}_{\mathbf{p}}.$

Since

$\displaystyle\sum_{f\in\mathcal{B}^{0}_{k}}R_{f}(X,Y)f(\tau)=\sum_{f\in\mathcal{B}^{0}_{k}}(-1)^{t}\frac{R_{f}^{ev}(X)R_{f}^{ev}(Y)+R_{f}^{ev}(Y)R_{f}^{ev}(X)}{D^{k-\frac{1}{2}}(2i)^{t(k-3)}<f,\,f>}f(\tau)$

we get

$\frac{\bigl{<}\,b_{\mathbf{k}}(X,Y;\cdot),\,f(\cdot)\,\bigr{>}}{\bigl{<}\,f,\,f\bigl{>}}=R_{f}(X,Y).$

Combining all together with Proposition 3.5 we conclude that

$\displaystyle C(X,Y;\tau;T)=F_{\tau}(T,-XYT)\,F_{\tau}(XT,YT)$

which completes a proof. ∎

##  5\. Conclusion

One of the main importance of modular forms in number theory is that spaces of
modular forms are generated by those with rational Fourier coefficients. The
”period theory” gives another natural rational structure of modular forms. A
striking result by Zagier [24] states that this rational information of
modular forms can be written as a single product of Kronecker series
$F_{\tau}(u,v)$ which is a Jacobi form. The recent results in [2, 21] show
that Eisenstein-Kronecker numbers have a rich arithmetic nature, such as a
connection with the special Hecke $L$-function over imaginary quadratic fields
and Katz’ two-variable p-adic Eisenstein measure.
In this paper, we identified the Kronecker series as a ”Kuznetsov lifting” of
holomorphic Hilbert Eisenstein series over totally real number fields with
strict class number 1. This is the first case to connect Kronecker series to
the critical values of Hilbert modular $L$-functions over a totally real
number field and it seems worthwhile to explore the hidden arithmetic
relations more.

On the other hand, in terms of geometric interpretation, a modular form can be
regarded as a section of a certain sheaf of differential forms on the open
modular curve on a congruence subgroup $\Gamma\subset SL_{2}(\mathbb{Z}).$ By
noting that the singular cohomology of the open modular curve is given by the
group cohomology $H^{*}(\Gamma,W)$ the comparison of de Rham and singular
cohomology can give an Eichler-isomorphism. Matsushima and Murakami [15]
extended the results to show that the space of automorphic forms on a
symmetric space $M$ is isomorphic to $H^{*}(M,S)$ for a certain locally
constant sheaf $S$ over $M.$ The cohomology of Hilbert surfaces in terms of
Hilbert cusp forms has been studied by many researchers including [11, 22].
Relating the critical values of $L$-functions of Hilbert cusp forms and
cohomology was first studied by Yoshida [22]. Following the work by Knopp [13]
and Kohnen-Zagier [14], which provide us the considerable new lights on
Eichler-Shimura isomorphisms such as rational structures of elliptic modular
forms, we are able to associate the space of Hilbert modular forms over the
totally real number fields to the parabolic cohomology group in terms of the
period polynomial by taking anti-derivative of Hilbert modular form [5].

## References

  * [1] A. Babei, L. Rolen and I.Wagner, The Riemann hypothesis for period polynomials of Hilbert modular forms, Journal of Number Theory, Vol. 218, January 2021, page 44 - 61.
  * [2] K. Bannai and S. Kobayashi, Algebraic theta functions and the p-adic interpolation of Eisenstein-Kronecker numbers, Duke Mathematical Journal, 153 (2010), no. 2, 229 - 295.
  * [3] G. Bol, Invarianten linearer Differentialgleichungen, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 16 (1949), no. 3-4, 1 - 28.
  * [4] J. Brunier, G. van der Geer, G. Harder and D. Zagier, The $1-2-3$ of modular forms, Lectures from the Summer School on Modular Forms and their Applications held in Nordfjordeid, June 2004, Edited by Kristian Ranestad, Universitext. Springer-Verlag, Berlin, 2008.
  * [5] Y. Choie,  Parabolic cohomology and Hilbert modular forms, Preprint (2021) .
  * [6] Y. Choie, H. Kim and O. Richter, Differential operators on Hilbert modular forms, Journal of Number Theory 122 (2007), no. 1, 25-36.
  * [7] Y. Choie and M. Lee, Jacobi-like forms, Pseudodifferential operators, and Quasimodular forms, Springer Monographs in Mathematics (2019).
  * [8] Y. Choie, Y. Park and D. Zagier, Periods of modular forms on $\Gamma_{0}(N)$ and Products of Jacobi Theta functions, Journal of the European Mathematical Society, Vol. 21, Issue 5, 1379 - 1410 (2019).
  * [9] M. Eichler, Eine Verallgemeinerung der Abelschen Integrale, Mathematische Zeitschrift, vol. 67 (1957) pp 267-298.
  * [10] M. Eichler and D. Zagier, The theory of Jacobi forms, Progress in Mathematics, 55, Birkhäuser Boston, Inc., Boston, MA, 1985.
  * [11] G. van der Geer, Hilbert modular surfaces, Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)], 16, Springer-Verlag, Berlin, 1988.
  * [12] M. Harris, $L$-functions of $2\times 2$ unitary groups and factorization of periods of Hilbert modulr forms, Journal of American Mathematical Society 6 (1993), no. 3, 637 - 719.
  * [13] M. Knopp, Some new results on the Eichler cohomology of automorphic forms, Bulletin of the American Mathematical Society, 80 (1974), 607 - 632.
  * [14] W. Kohnen and D. Zagier, Modular forms with rational periods, Modular forms (Durham, 1983), 197-249, Ellis Horwood Series in Mathematics and its Applications: Statistics and Operational Research, Horwood, Chichester, 1984.
  * [15] Y. Matsushima and S. Murakami, On vector bundle valued harmonic forms and automorphic forms on symmetric Riemannian manifolds, Annals of Mathematics (2) 78 (1963), 365 - 416.
  * [16] Y. Manin, Periods of cusp forms, and $p$-adic Hecke series, (Russian) Matematicheskii Sbornik (N.S.) 92 (134) (1973), 378 - 401, 503.
  * [17] $\underline{\makebox[40.00006pt]{}}$, On the critical values of certain Dirichlet series and the periods of automorphic forms, Inventiones Mathematicae 94 (1988), no. 2, 245 - 305.
  * [18] $\underline{\makebox[40.00006pt]{}}$, Algebraic relations between critical values of zeta functions and inner products, American Journal of Mathematics, 105 (1983), no. 1, 253 - 285.
  * [19] $\underline{\makebox[40.00006pt]{}}$, The special values of the zeta functions associated with Hilbert modular forms, Duke Mathematical Journal, 45 (1978), no. 3, 637- 679.
  * [20] $\underline{\makebox[40.00006pt]{}}$, Sur les intégrales attachées aux forms automorphes, Journal of the Mathematical Society of Japan 11 (1959), 291 - 311.
  * [21] J. Sprang,  Eisenstein-Kronecker series via the Poincaré bundle, Forum of Mathematics, Sigma 7 (2019), No. e34, 59 pp.
  * [22] H. Yoshida, Absolute CM-periods. Mathematical Surveys and Monographs, 106, American Mathematical Society, Providence, RI, 2003.
  * [23] A. Weil, Elliptic functions according to Eisenstein and Kronecker, Reprint of the 1976 original, Classics in Mathematics, Springer-Verlag, Berlin, 1999.
  * [24] D. Zagier, Periods of modular forms and Jacobi theta functions, Inventiones Mathematicae, 104, 449-265 (1991).
  * [25] $\underline{\makebox[40.00006pt]{}}$, The Rankin-Selberg method for automorphic functions which are not of rapid decay, Journal of the Faculty of Science, University of Tokyo, Section IA, Mathematics 28 (1981), no. 3, 415 - 437 (1982).
  * [26] $\underline{\makebox[40.00006pt]{}}$, On the values at negative integers of the Zeta-functions of a real quadratic fields, L’Enseignement Mathématique, 22 (1976), no. 1 - 2, 55 - 95.