# Global Analysis of LISA Data with Galactic Binaries and Massive Black Hole
Binaries

Stefan H. Strub<EMAIL_ADDRESS>Luigi Ferraioli  Cédric Schmelzbach
Simon C. Stähler  Domenico Giardini  Institute of Geophysics, ETH Zurich
Sonneggstrasse 5, 8092 Zurich, Switzerland

###### Abstract

The Laser Interferometer Space Antenna (LISA) is a planned space-based
observatory to measure gravitational waves in the millihertz frequency band.
This frequency band is expected to be dominated by signals from millions of
Galactic binaries and tens of merging massive black hole binaries. The LISA
Data Challenge 2a is focused on robust signal extraction from a blend of these
two types of gravitational wave signals. Here, we introduce a novel high
performance and cost-effective global fit pipeline extracting and
characterizing galactic binary and massive black hole binary signals and
estimate the noise of the residual. We perform the pipeline in a time-evolving
weekly analysis starting with an observation time of
$1\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$ until we reach a full
year. As expected we detect more galactic binaries and massive black hole
binaries bringing the noise estimate of the residual closer to the instrument
noise with each week of additional observation time. Furthermore, we present a
novel maximum likelihood estimate-based algorithm for extracting multiple
massive black hole binaries. Additionally, we demonstrate a massive black hole
binary signal extraction with a more accurate LISA response, considering
higher harmonic modes, in a noisy data set.

Gravitational Waves, LISA, Galactic Binaries, Massive Black Hole Binaries,
LISA Data Challenge, Global Fit

##  I Introduction

The detection of gravitational waves (GWs) by the LIGO detector in 2015 marked
a significant breakthrough in astrophysics [1]. This breakthrough and the
successful technology demonstration with LISA Pathfinder paved the way for
LISA, a space-borne observatory designed to observe low-frequency GWs in the
range from $0.1\text{\,}\mathrm{mHz}$ to $1\text{\,}\mathrm{Hz}$ [2].

In the LISA frequency band, the most common sources are tens of millions of
Galactic binaries (GBs), primarily white dwarf pairs, emitting quasi-
monochromatic GWs. Unlike merging black holes, GBs remain separate throughout
their lifetimes, allowing us to observe GBs with LISA continuously. It is
estimated that the analysis of the LISA data will reveal tens of thousands of
these interfering signals, with the remainder forming a galactic foreground
noise [3]. By precisely measuring the parameters of the GBs, valuable insights
into the evolution of binary systems will be gained [4, 5, 6, 7, 8].

In addition to signals from GBs, LISA will also detect GWs from massive black
hole binaries (MBHBs), with masses ranging from
$10^{5}\text{\,}\textup{M}_{\odot}$ to $10^{7}\text{\,}\textup{M}_{\odot}$
[9]. Unlike the GBs’ quasi-monochromatic GW signals, MBHBs mergers are
recorded as a chirp with a high signal-to-noise ratio (SNR) making them
especially prominent in the data collected [10]. Estimates suggest LISA could
detect 1 to 20 MBHB mergers annually [11, 12], offering a window into
fundamental physics, cosmology, and astrophysics [13, 14, 15, 9].

Analyzing the LISA data presents a difficult task due to the overlap of
numerous signals in the same data stream. To address this challenge, the LISA
Data Challenge (LDC) has emerged as a crucial platform [16]. Building on the
legacy of the Mock LISA Data Challenges (MLDC) [17, 18, 19], the LDC provides
a comprehensive suite of simulated datasets mimicking real-world LISA
observations. These datasets serve as a testing ground for researchers to
develop and refine their data analysis methods.

The LDC is structured into various challenges, each designed to explore a
specific aspect of LISA data analysis. LDC1, the first set of challenges,
presents researchers with six distinct scenarios, each focusing on a single
type of GW signal. For instance, LDC1-1 contains a single massive black hole
binary (MBHB) signal along with simulated instrument noise. LDC2a presents a
more complex challenge, tasking researchers with a dataset containing a
combination of 15 MBHBs, 30 million GBs, and instrument noise where the signal
extraction and parameter estimation have to be a global fit of all parameters.

Several methods have been proposed for extracting GB signals, including
maximum likelihood estimate (MLE) and Bayesian approaches: MLE methods are
used to find the best matching simulated signal to the data [20, 21, 22],
while Bayesian methods provide a posterior distribution that describes the
uncertainty of the source parameters. Successful approaches to Bayesian
analysis are Markov chain Monte Carlo (MCMC) based methods. Examples for MCMC
based algorithms are blocked annealed Metropolis-Hastings (BAM) [23, 24, 25],
an MCMC algorithm with simulated annealing, or the reversible jump Markov
chain Monte Carlo (RJMCMC) [26, 27] method, which allows for varying parameter
dimensions and thus variable numbers of GBs to construct the posterior
distribution. Recently, a hybrid method has been introduced that leverages
both MLE and MCMC [28, 29]. This approach is used to first find the MLE and
then to compute the posterior distribution around the identified MLE with a
parallelized MCMC method.

Methods to analyze the MBHB signals have historically used both MCMC-based or
MLE-based algorithms, as detailed in various studies [30, 31, 32, 33, 34, 35,
36, 37, 38]. In recent years, the focus has shifted from MLE algorithms
towards MCMC algorithms, offering robust posterior distribution computations
[39, 40, 41, 42, 10]. The analysis with a more accurate LISA response
including higher harmonic modes was recently explored in a noiseless data set
[10]. Additionally, the impact of data gaps on the analysis of MBHBs,
specifically in scenarios without noise, has been investigated [43].

Recently, a prototype global analysis of MBHBs and GBs has been published
[27]. The pipeline is a blocked MCMC sampler to obtain a simultaneous global
fit of the LDC2a data set. The method to extract the MBHBs is a parallel
tempering MCMC algorithm [42], and the GB extraction is based on a RJMCMC
algorithm [26].

In this article, we present a novel pipeline to analyze together MBHB and GBs.
The performances of the pipeline are evaluated with the LDC2a dataset.
Additionally, we propose a novel pipeline for extracting MBHBs within LISA
data. We revisit MLE approaches by introducing a robust genetic algorithm-
driven pipeline. This pipeline effectively extracts MBHBs from noisy data,
incorporating scenarios with and without higher modes. The global fit pipeline
combines the newly developed MBHB extraction method with the established GB
algorithm from our previous work [29] and a generalized noise estimate
algorithm.

In Sec. II we introduce the likelihood function, while in Sec. III we present
the MBHB algorithm, changes made on the previous GBs algorithm, and describe
the global fit pipeline. Our results, presented in Section IV, include testing
the MBHB algorithm on LDC1-1 and a self-injected higher mode signal, along
with the global fit pipeline on the LDC2a dataset. Additionally, Section V
showcases the results of robustness tests of the pipeline. Finally, in Section
VI, we engage in a discussion regarding the performance and future extensions
of the pipeline.

##  II The Likelihood function

LISA data $d(t)$ will consist of the sum of multiple GW signals $s(t,\theta)$
and noise $n$ which can be written as $d(t)=\sum s(t,\theta)+n(t)$. For
simplicity, we shall henceforth forego the time dependency notation, although
$d$, $s(\theta)$, and $n$ remain time series.

The likelihood function $p(d|\theta)$ serves to quantify the probability of
observing the data stream $d$ given the source parameters $\theta$
characterizing the GW signal. Employing the logarithm of the likelihood,
denoted as $\log p(d|\theta)$, is prevalent owing to its computational
tractability and is expressed as:

$\log p(d|\theta)=-\frac{1}{2}\langle d-s(\theta)|d-s(\theta)\rangle,$ (1)

where $\langle x(t)|y(t)\rangle$ signifies the scalar product between two time
domain signals $x(t)$ and $y(t)$. It is defined as:

$\langle
x(t)|y(t)\rangle=4\mathcal{R}\left(\int_{0}^{\infty}\frac{\tilde{x}(f)\tilde{y}^{*}(f)}{S(f)}\,df\right).$
(2)

Here, $\tilde{x}(f)$ denotes the Fourier transform of $x(t)$, and $S(f)$
represents the one-sided power spectral density of the noise. The noise
estimation is continuously refined during the analysis process, with the
specifics regarding noise estimation discussed in Sec.III.1.

To mitigate laser noise inherent in the laser measurements of LISA’s arms, the
technique of time-delay-interferometry (TDI) is employed, converting the six
Doppler measurements into three observables: X, Y, and Z [44, 45, 46, 47, 48].
Consequently, the data $d$ and the signal $s(\theta)$ entail TDI responses
across multiple channels. The scalar product here reads as:

$\langle
d-s\left(\theta\right)|d-s\left(\theta\right)\rangle=\sum_{\alpha\in\mathcal{M}}\langle
d_{\alpha}-s_{\alpha}\left(\theta\right)|d_{\alpha}-s_{\alpha}\left(\theta\right)\rangle$
(3)

where $\mathcal{M}={X,Y,Z}$ represents the default TDI combination.
Alternatively, $\mathcal{M}={A,E,T}$ where:

$\displaystyle A$ $\displaystyle=\frac{1}{\sqrt{2}}\left(Z-X\right)$ (4)
$\displaystyle E$ $\displaystyle=\frac{1}{\sqrt{6}}\left(X-2Y+Z\right)$
$\displaystyle T$ $\displaystyle=\frac{1}{\sqrt{3}}\left(X+Y+Z\right)$

These channels are specifically chosen for their lack of correlation with
instrument noise [49, 50]. While we utilize $A$, $E$, and $T$, we only
consider $A$ and $E$ for signals with frequencies $f<f_{\ast}/2=1/(4\pi
L)\approx$9.55\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$$, given the suppression
of the GW response for $T$ where $f_{\ast}$ is the transfer frequency and $L$
is the distance between two spacecraft of LISA [26].

##  III Time evolving global fit

In our pipeline to analyze data containing signals of multiple source types,
we first discuss the noise estimation and the signal extraction methods
individually. In the last subsection III.4, we present the global fit pipeline
combining the individual algorithms.

###  III.1 Noise Estimation

The noise estimation process is performed across the full frequency band and
undergoes a modified version as outlined in [29]. Here, we provide a concise
overview of the pipeline along the modifications.

Initially, if any GW signals have been previously detected, we subtract them
from the data. Specifically, in our scenario, we subtract the recovered GBs
denoted as $\tilde{\theta}_{\text{GB,recovered}}$ and recovered MBHBs denoted
as $\tilde{\theta}_{\text{MBHB,recovered}}$ from the data, and obtain the
residual

$\begin{split}d_{\text{residual}}&=d-\sum\limits_{\theta\in\tilde{\theta}_{\text{GB,recovered}}}s(\theta)-\sum\limits_{\theta\in\tilde{\theta}_{\text{MBHB,recovered}}}s(\theta)\\\
&=d-d_{\text{GB}}-d_{\text{MBHB}}\end{split}$ (5)

Here, $s(\theta)$ represents the signal corresponding to the MLE of each
recovered parameter set $\theta$, and $d_{\text{GB}}$ and $d_{\text{MBHB}}$
represent the sum of all GB and MBHB signals, respectively.

Subsequently, the residual data $d_{\text{residual}}$ is utilized to estimate
a noise curve across the frequency domain employing Welch’s method. To address
outlier peaks, particularly those associated with unresolved signals, a
smoothing procedure is implemented. A frequency window comprising 30 bins is
defined, and any values exceeding double the window’s median are adjusted to
the median value. This process is reiterated by shifting the window by 15
frequency bins until the entire Power Spectral Density (PSD) is smoothed.
Further smoothing is conducted using the Savitzky-Golay filter. Spline
interpolation is then applied to derive a PSD estimate, leading to the
residual noise curve.

Different from the algorithm described in [29], while implementing Welch’s
method, we transition from a fixed number of 500 windows to a fixed number of
$15\,000$ samples per window. This change aims to generalize the noise
estimate for varying observation times, considering that the number of
frequency bins changes with $T_{\text{obs}}$. Consequently, the window length
for the Savitzky-Golay filter is fixed at 10 frequency bins for frequencies
below $0.8\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ and at 70 bins for
frequencies above, facilitating a consistent approach across different
observation duration $T_{\text{obs}}$.

Figure 1: The noise estimates and PSD of the TDI A channel within the
$1\text{\,}\mathrm{y}\mathrm{r}$ LDC2a dataset are shown.
$S_{A\text{,instrument}}$ represents the noise PSD utilized in generating the
dataset. The discrepancy observed between $S_{A\text{,residual}}$ and the
instrument PSD within the $0.2\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ to
$5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ range is attributed to unresolved
Galactic Backgrounds (GBs). It is anticipated that the majority of GBs within
this frequency span remain unresolved, leading them to blend into the Galactic
foreground noise. The blue spike is an unsolvable GB signal at
$26\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ due to an issue with the dataset
instead of the analysis. At [51] we provide videos and figures of the noise
estimate for each week.

In Fig. 1 we present the noise estimate for the LDC2a dataset. The noise
estimates generally exhibit alignment with instrument noise, except for the
$0.2\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ to
$4\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ range, where deviations occur due
to unresolved background signals. Similar computation procedures are employed
for noise estimation in other TDI channels ($E$ and $T$).

###  III.2 Extracting Massive Black Hole Binary signals

To simulate the LISA response of a GW emitted by a MBHB we use the BBHx [41,
10] algorithm, accessible at [52]. This simulator utilizes 11 parameters
$\theta_{\text{MBHB}}=\left\\{\beta,\lambda,M,q,a_{1},a_{2},D_{L},\iota,\phi_{0},\psi,t_{\text{ref}}\right\\}$
to model the GW signal. Here, $\lambda$ and $\beta$ correspond to the sky
coordinates in terms of ecliptic longitude and ecliptic latitude,
respectively. The total mass $M=m_{1}+m_{2}$ represents the sum of the two
masses of the binary, and $q$ denotes the mass ratio $q=m_{2}/m_{1}$. The
spins $a_{1}$ and $a_{2}$ pertain to the objects with masses $m_{1}$ and
$m_{2}$, respectively, and are aligned with the orbital angular momentum of
the system; a negative spin indicates anti-parallel alignment. $D_{L}$ denotes
the luminosity distance, $\iota$ represents the inclination angle, $\phi_{0}$
signifies the initial phase, and $\psi$ corresponds to the polarization angle.
The parameter set $\theta_{\text{MBHB}}$ is defined within the prior
$\Theta_{\text{MBHB}}$ as shown in Table 1. The prior is linearly uniform
except for $\beta$, $\iota$, and $M$, which are sine-, cosine-, and log-
uniform, respectively.

Table 1: Boundaries of the prior distribution $\Theta_{\text{MBHB}}$ for MBHBs. Parameter | Lower Bound | Upper Bound
---|---|---
$\sin\beta$ | $-1$ | 1
$\lambda$ | $0$ | $2\pi$
$\log M$ ($\textup{M}_{\odot}$) | $\log 10^{6}$ | $\log 10^{7}$
$q$ | $1$ | $10$
$a_{1}$ | $-1$ | $1$
$a_{2}$ | $-1$ | $1$
$D_{L}$ (Gpc) | 0.5 | 1000
$\cos\iota$ | $-1$ | 1
$\phi_{0}$ | $-\pi$ | $\pi$
$\psi$ | 0 | $\pi$
$t_{\text{ref}}$ | $\min(t_{\text{segment}})$ | $\max(t_{\text{segment}})$

To obtain the MLE we can maximize the signal-to-noise ratio (S/N) defined as

$\rho=\frac{\langle d|s\left(\theta^{\prime}\right)\rangle}{\sqrt{\langle
s\left(\theta^{\prime}\right)|s\left(\theta^{\prime}\right)\rangle}}=\frac{\langle
d|s\left(\theta\right)\rangle}{\sqrt{\langle
s\left(\theta\right)|s\left(\theta\right)\rangle}}.$ (6)

which is independent of $D_{L}$ with
$\theta^{\prime}=\theta\setminus\\{D_{L}\\}$ and obtain analytically the
luminosity distance which maximizes the likelihood with

$D_{L}=\frac{\langle
s\left(\theta^{\prime}\right)|s\left(\theta^{\prime}\right)\rangle}{\langle
d|s\left(\theta^{\prime}\right)\rangle}.$ (7)

The MBHB signal extraction uses a time-evolving algorithm that extracts MBHBs
within LISA data time segments. This approach reflects a more realistic
scenario where the inspiral phase of MBHBs may be present within the data but
remain unextracted since the merger has not been recorded yet.

In Algorithm 1, we provide pseudo-code for the search algorithm, using
differential evolution (DE) [53] to maximize eq. 6, aiming to obtain the MLE
of MBHBs within a time series segment $d_{\text{segment}}$ spanning a set time
$t_{\text{segment}}$ of the data. To mitigate the potential issue of a merger
occurring at the border of $d_{\text{segment}}$, we extend the segment to
analyze $d_{\text{analyze}}$, incorporating padded regions before and after
the segment. Both the segment and padding lengths are adjustable, but for this
work, we use a one-week segment with one week of padding before and one day of
padding after the segment, resulting in a total analysis window length of 15
days. The padding before the merger can be extended as desired until the start
of measurement since MBHBs that merged before $d_{\text{segment}}$ have
already been extracted. Conversely, the padding duration after
$d_{\text{segment}}$, being in the future, is limited by available data. We
deem a padding duration of $1\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}$ after
$d_{\text{segment}}$ to be sufficient.

$\textbf{Function}\ \textbf{{MBHB\\_search}}(d_{\text{analyze}})$

$\tilde{\theta}_{\textrm{found}}\leftarrow\\{\,\\}$

$d_{\textrm{residual}}\leftarrow d_{\textrm{analyze}}$

$\rho_{\textrm{max}}\leftarrow\rho_{\textrm{threshold}}+1$

while $\rho_{\textrm{max}}>\rho_{\textrm{threshold}}$ do

$\tilde{\theta}^{\prime}_{\textrm{MLEs}}\leftarrow\\{\,\\}$

for $j$ in $\\{1,2,...,n_{\textrm{searches}}\\}$ do

$\theta^{\prime}_{\textrm{init}}$ randomly drawn from prior

$\theta^{\prime}_{\textrm{MLE}}\leftarrow\operatorname*{\arg\\!\max}\limits_{\theta^{\prime}}\rho(\theta^{\prime},d_{\textrm{residual}})$
using DE with $\theta^{\prime}_{\textrm{init}}$

if
$\rho(\theta^{\prime}_{\textrm{MLE}},d_{\textrm{analyze}})\leq\rho_{\textrm{threshold}}\textbf{
and }j>1\textbf{ do}$

$\textbf{ return }\tilde{\theta}_{\textrm{found}}$

$\tilde{\theta}^{\prime}_{\textrm{MLEs}}\leftarrow\tilde{\theta}^{\prime}_{\textrm{MLEs}}\cup\\{\theta^{\prime}_{\textrm{MLE}}\\}$

end for

$\theta^{\prime}_{\textrm{MLE}}\leftarrow\operatorname*{\arg\\!\max}\limits_{\theta^{\prime}\in\tilde{\theta}^{\prime}_{\textrm{MLEs}}}\rho(\theta^{\prime},d_{\textrm{residual}})$

$\rho_{max}\leftarrow\rho(\theta^{\prime}_{\textrm{MLE}},d_{\textrm{analyze}})$

if $\rho_{max}\leq\rho_{\textrm{threshold}}\textbf{ do}$

$\textbf{return }\tilde{\theta}_{\textrm{found}}$

Compute $D_{L}$ according to eq. (7) with $\theta^{\prime}_{\textrm{MLE}}$

$\theta_{\textrm{MLE}}\leftarrow\theta^{\prime}_{\textrm{MLE}}\cup\\{D_{L}\\}$

$\tilde{\theta}_{\textrm{found}}\leftarrow\tilde{\theta}_{\textrm{found}}\cup\\{\theta_{\textrm{MLE}}\\}$

$d_{\textrm{residual}}\leftarrow
d_{\textrm{residual}}-s(\theta_{\textrm{MLE}})$

end while

$\textbf{return }\tilde{\theta}_{\textrm{found}}$

Algorithm 1 The MBHB search algorithm to analyze a data segment
$d_{\text{analyze}}$.

To obtain a good fit with a high likelihood we set the number of searches for
each signal to $n_{\textrm{searches}}=5$. This repeats each search 5 times
with different random initial parameters $\theta^{\prime}_{\textrm{init}}$ and
then selects the result with the highest likelihood.

The MLE can then be used to subtract the MBHBs from the data containing
signals of multiple source types to analyze the other signals like GBs without
the presence of the MBHBs. Furthermore, the MLE can be used to obtain an
uncertainty estimate through the Fisher Information Matrix or any MCMC
algorithm with the MLE as its initial guess to skip the burn-in phase.

The uncertainty estimate with FIM is generally not valid by wrongfully
assuming Gaussian distribution of errors [38]. For completeness, we still use
the FIM to provide a fast and simple uncertainty estimate even though an
uncertainty estimate is not necessarily needed for the global fit. We compute
the uncertainties by taking the square root of the inverse of the Fisher
Information Matrix (FIM), defined as:

$F_{ij}=\langle\partial_{i}s(\theta)|\partial_{j}s(\theta)\rangle,$ (8)

where $\partial_{i}$ represents the partial derivative concerning the $i$th
component of the parameter vector $\theta$. The derivatives of the FIM are
computed using the second-order center finite difference method with a step
size of $10^{-9}$ times the exploration space determined by the prior
distribution $\Theta_{\text{MBHB}}$. The estimated uncertainty vector is
denoted as $\sigma=\sqrt{\text{diag}(F^{-1})}$.

###  III.3 Extracting Galactic Binary signals

The GB extraction algorithm, to construct $d_{\text{GB}}$, is described in an
earlier paper [29]. The algorithm splits the frequency domain data into
multiple segments $B_{\text{search}}$ to be analyzed. Where one segment is
double the size of the maximum expected bandwidth of a signal within the
segment. The segments are split into two groups of non-neighboring segments
called $B_{\text{even}}$ and $B_{\text{odd}}$. All even segments are first
analyzed in parallel and then the $B_{\text{odd}}$ segments are analyzed with
the signals found within the even segments subtracted. Lastly, the even
segments are reanalyzed with the signals found within the odd segments
subtracted.

Simulating a GB signal involves eight parameters
$\theta_{\text{GB}}=\left\\{\mathcal{A},\lambda,\beta,f,\dot{f},\iota,\phi_{0},\psi\right\\}$
[54]. Here $\mathcal{A}$ represents the amplitude, $\lambda$ and $\beta$
correspond to the sky coordinates in terms of ecliptic longitude and ecliptic
latitude, respectively. The parameter $f$ signifies the frequency of the GW,
while $\dot{f}$ denotes the first-order frequency derivative. The inclination
angle is denoted by $\iota$, $\phi_{0}$ signifies the initial phase and $\psi$
corresponds to the polarization angle. The prior distribution is listed in
Table 2.

Table 2: Boundaries of the prior distribution $\Theta$. Parameter | Lower Bound | Upper Bound
---|---|---
$\sin\beta$ | $-1$ | 1
$\lambda$ | $-\pi$ | $\pi$
$f$ | $0.3\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ | $f_{\text{Nyquist}}$
$\dot{f}$ | $-5\cdot 10^{-6}f^{13/3}$ | $1.02\cdot 10^{-6}f^{11/3}$
$\log\mathcal{A}$ | $\log\mathcal{A}(\rho=7)$ | $\log\mathcal{A}(\rho=1000)$
$\cos\iota$ | $-1$ | 1
$\phi_{0}$ | 0 | 2$\pi$
$\psi$ | 0 | $\pi$

The amplitude boundary is determined by a lower and upper bound signal to
noise ratio $\rho$ and is related to the amplitude by [26]

$\mathcal{A}\left(\rho\right)=2\rho\left(\frac{S\left(f\right)}{T_{obs}\,\sin^{2}\left(f/f_{\ast}\right)}\right)^{1/2}.$
(9)

###  III.4 Pipeline of the time evolving global fit

Combining individual algorithms to analyze data sets containing signals of
multiple source types, we introduce a time-evolving pipeline that dynamically
updates signal extraction as new data arrives.

In Fig. 2 we illustrate the pipeline of one iteration through a flow diagram.
For a likelihood-based signal extraction we need a noise estimate. Initially,
we estimate the noise using the original data since at the beginning no
signals are extracted and therefore $d_{\mathrm{MBHB}}=d_{\mathrm{GB}}=0$ for
the first iteration.

After estimating the noise we first extract the MBHBs since a new MBHB merger
could have been recorded, and these mergers are so loud that they would tamper
significantly with the GB extraction. Moreover, since the GBs are quasi
monochromatic, their signal does not change much from one time step to the
next one. Therefore the GB extraction from the previous time set is already a
good fit to the updated data and does not tamper much with the MBHB
extraction.

In each iteration, the complete dataset available up to that moment is
reanalyzed. With the updated $d_{\text{GB}}$ incorporating an increased number
of subtracted GB signals, there is reduced interference in the MBHB signal
extraction process. Consequently, it is recommended to reassess all MBHBs as
well to ensure optimal analysis outcomes.

Subsequently, we subtract the identified MBHBs $d_{\text{MBHB}}$ from the data
and proceed to extract the GBs. The GB algorithm autonomously generates an on-
the-fly noise estimate for each frequency segment, obviating the need for the
noise estimate $n$ across the full frequency spectrum [29].

Once new data of one week arrives, we update the dataset $d$ and start the
pipeline shown in Fig. 2 for another iteration.

Figure 2:  Flow chart of one iteration of the time-evolving extraction of
MBHBs and GBs and noise estimate. The parallelograms represent the inputs for
the processes, which are represented by rectangles. For the first iteration
$d_{\mathrm{MBHB}}=d_{\mathrm{GB}}=0$ which are then updated by the
pipeline.Start$d-d_{\mathrm{MBHB}}-d_{\mathrm{GB}}$ Estimate Noise $n$
$d-d_{\mathrm{GB}}$, $n$ Extract MBHBs $d_{\mathrm{MBHB}}$
$d-d_{\mathrm{MBHB}}$ Extract GBs $d_{\mathrm{GB}}$ Stop

##  IV Results

Since the algorithms for the noise estimate and GBs have already been
discussed in [28, 29] we now first show the results of the MBHB part
individually. Later we show the results for the LDC2a containing overlapping
GB and MBHB signals.

###  IV.1 MBHBs

Table 3: The injected and recovered parameters of the MBHB of the LDC1-1 data challenge. The uncertainties are computed with the inverse of the Fisher Information Matrix. Parameter | Injected | Recovered (2,2) mode | Recovered higher modes
---|---|---|---
$m_{1}$ ($\textup{M}_{\odot}$) | $2.599137\times 10^{6}$ | $2.60\times 10^{6}\pm 2\times 10^{4}$ | $2.523\times 10^{6}\pm 5\times 10^{3}$
$m_{2}$ ($\textup{M}_{\odot}$) | $1.242861\times 10^{6}$ | $1.239\times 10^{6}\pm 8\times 10^{3}$ | $1.275\times 10^{6}\pm 2\times 10^{3}$
$a_{1}$ | $7.534822\times 10^{-1}$ | $7.6\times 10^{-1}\pm 2\times 10^{-2}$ | $6.66\times 10^{-1}\pm 3\times 10^{-3}$
$a_{1}$ | $6.215875\times 10^{-1}$ | $6.0\times 10^{-1}\pm 5\times 10^{-2}$ | $8.33\times 10^{-1}\pm 4\times 10^{-3}$
$D_{L}$ (Gpc) | $5.600578\times 10^{1}$ | $5\times 10^{1}\pm 4\times 10^{1}$ | $3.5\times 10^{1}\pm 4\times 10^{0}$
$\iota$ | $1.224532\times 10^{0}$ | $1.3\times 10^{0}\pm 3\times 10^{-1}$ | $1.35\times 10^{0}\pm 3\times 10^{-2}$
$\lambda$ | $3.509100\times 10^{0}$ | $3.49\times 10^{0}\pm 7\times 10^{-2}$ | $3.526\times 10^{0}\pm 4\times 10^{-3}$
$\beta$ | $2.926963\times 10^{-1}$ | $2.7\times 10^{-1}\pm 6\times 10^{-2}$ | $2.86\times 10^{-1}\pm 3\times 10^{-3}$
$t_{c}$ (s) | $2.496000\times 10^{7}$ | $2.495999\times 10^{7}\pm 3\times 10^{1}$ | $2.4960004\times 10^{7}\pm 2\times 10^{0}$

Further, we test the pipeline on the signal simulated considering higher modes
$(l,m)=\\{(2,2),(3,3),(4,4),(2,1),(3,2),(4,3)\\}$ where the MBHB simulation of
the LDC1-1 and LDC2a data sets only considers the primary $(l,m)=(2,2)$ mode
[10].

In the following, we present the results of the MBHB algorithm described in
III.2 for the LDC1-1 data set, which consists of instrument noise and a single
MBHB signal. As mentioned in [55] the signal simulation also used in this work
sets $t_{\text{ref}}$ at $f_{\text{max}}$ where $f_{\text{max}}$ is set
internally which is different from the signal in LDC1-1 is simulated where
$f^{2}A_{22}$ is maximized instead. Therefore to avoid discrepancies due to
using different simulators we subtract the MBHB signal from the LDC data and
self-inject the signal with the same parameters using BBHX for the signal
simulations.

Figure 3: LDC1-1 TDI A component at the merger of the MBHB signal. The black
line is the data. The thicker orange line is the self injected signal
according to the parameters given by the LDC catalog and the thin green line
is the recovered signal.

Since the observation time $T_{\text{obs}}=1yr$ is rather long we split the
data into segments of $1\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$
duration.

In Fig. 3 we present the found signal alongside the injection and the data.
All analyses have been conducted using noisy data. The padding is usually set
to 7 days before and 1 day after the segment but for the analysis considering
higher harmonic modes, we obtained better results by extending the padding
before the segment to 28 days.

In Table 3 we present the recovered parameters and their uncertainties. Since
we compute the uncertainties by taking the inverse of the FIM we wrongfully
assume a Gaussian posterior distribution for all parameters. Nonetheless, the
obtained uncertainties are consistent with uncertainties computed with an MCMC
sampler as presented in [10].

###  IV.2 Multiple type signal sources

Here we show the results of our analysis of the LDC2a containing overlapping
GB and MBHB signals and instrument noise. We analyzed for GBs across the
entire frequency spectrum
$f\in[$0.3\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$,$100\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$]$.
This range spans from regions where GB detection is improbable due to
instrument noise to the Nyquist frequency. The extraction of MBHBs was carried
out in time domain segments of the entire one-year dataset.

Since the pipeline as described in section III.4 is a time-evolving algorithm
where the pipeline in Fig. 2 runs once for each time step. We decide to
analyze the data in weekly steps. Therefore we can set a realistic premise to
analyze data as new data is coming in as in a realistic setup with LISA being
online.

For both the MBHB and GB search algorithms we set $\rho_{\text{threshold}}=9$
as the threshold to accept a found signal within the set of recovered signals.
To simulate the LISA response of MBHBs we use the BBHX implementation [52] and
for the GBs we use fastGB provided by the LDC [16].

Figure 4: The weekly results of the mixed signal LDC2a data set. In black we
show the TDI A component of the LDC2a data set to be analyzed and in orange
and green we show the injected and resolved MBHBs $d_{\text{MBHB}}$. The blue
dots and circles are the recovered and matched GBs respectively. The purple
crosses are the match ratio of matched GBs divided by recovered GBs. Figure 5:
Histogram of the difference of resolved and injected parameters of matched
signals for different observation times. Note that for $f$ and $\mathcal{A}$
we present the relative difference. Figure 6:  The top plot illustrates the
cumulative distribution function of scaled error $\delta$ of all extracted
GBs, while the bottom plot depicts the survival function. Each color
represents a different number of weeks of observation time.

To analyze the accuracy of the recovered GB signals $\theta_{\text{rec}}$ we
can match them with the injected parameter sets $\theta_{inj}$. This is only
possible with simulated data. In other publications, the overlap $\mathcal{O}$
is used to define a metric of similarity between two signals

$\mathcal{O}(s(\theta_{\text{rec}}),s(\theta_{\text{inj}}))=\frac{\langle
s(\theta_{\text{rec}})),s(\theta_{\text{inj}})\rangle}{\sqrt{\langle
s(\theta_{\text{rec}}),s(\theta_{\text{rec}})\rangle\langle
s(\theta_{\text{inj}}),s(\theta_{\text{inj}})\rangle}}.$ (10)

which is independent of the amplitude of the signals. Therefore in this work,
we use the scaled error

$\delta(s(\theta_{\text{rec}}),s(\theta_{\text{inj}}))=\frac{\langle
s(\theta_{\text{rec}})-s(\theta_{\text{inj}}),s(\theta_{\text{rec}})-s(\theta_{\text{inj}})\rangle}{\sqrt{\langle
s(\theta_{\text{rec}}),s(\theta_{\text{rec}})\rangle\langle
s(\theta_{\text{inj}}),s(\theta_{\text{inj}})\rangle}}$ (11)

which is dependent on the amplitude of the signals and does not run the risk
of matching an injected parameter set with low amplitude which is part of the
foreground noise. Particularly at lower frequencies, where there are fewer
frequency bins and low amplitude signals are prevalent, there exists a
heightened risk of substantial overlap with a recovered signal, despite the
scaled error not indicating a match.

In Fig. 6 we observe a rapid change of the cumulative distribution function
until
$T_{\text{obs}}=$10\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
where the scaled error rapidly declines and for
$T_{\text{obs}}>$26\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
onwards the cumulative distribution function stabilized. The steady increase
in total counts with longer observation times is as expected. We define a
recovered GB as matched if it forms a pair with an injected parameter set
$\theta_{\text{inj}}$ where
$\mathcal{O}(s(\theta_{\text{rec}}),s(\theta_{\text{inj}}))<0.3$.

The results of the weekly analysis of the LDC2a data set are shown in Fig. 4.
As expected there is a steady increase in the amount of recovered signals as
the observation time increases.

It is noteworthy that from week 16 onwards, the match ratio (ratio of matched
signals to extracted signals) stabilizes within the range of $85\%$ and
$88\%$. Before this, the match ratio is lower, indicating that the recovered
signals do not correspond well to the injected ones. This observation is also
evident in Fig. 5, where analyses with
$T_{\text{obs}}<$16\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
exhibit a high number of signals with low accuracy in sky location. Also from
week 16 onwards, the number of recovered GBs and matched GBs increases by
roughly 129 and 115 for each week.

Even though there is a MBHB merging right at
$t_{c}=$19.06\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$$, the analysis
of the data set with an observation time of 19 weeks shows no change in the
trend of the amount of recovered signals nor on the match rate. This shows
that the pipeline successfully extracts GBs without interference from MBHB
which merge soon afterwards. This is because the signal of this MBHB until
$T_{\text{obs}}=$19\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
is only in the inspiral phase and its frequency is
$0.4\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ or less. This is well visible in
Fig.9 which we will discuss later.

Also observable from Fig. 5 is that for all parameters there is a clear trend
of decreasing errors with longer observation times, as expected. The mostly
flat histograms for $\phi_{0}$ and $\psi$ for
$T_{\text{obs}}<$10\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
stem from their intrinsic degeneracy [29]. Nonetheless, as $T_{\text{obs}}$
increases, it becomes clear that this degeneracy progressively diminishes.

Figure 7: Scatter plot of recovered GBs of the LDC2a data set of different
observation times. Each of the shown recovered signals matched with an
injected signal with $\delta<0.3$. The bottom plot shows an overlap of the
three analyses. Figure 8: Histogram of the recovered GBs of the LDC2a data set
of all weekly analyses.

In Fig. 7 we present a scatter plot of $\mathcal{A}$ and $f$ for three
different observation times. It is well noticeable that the number of
recovered GBs increases overall and that also GBs of lower $\mathcal{A}$ are
recovered. An increase of $T_{\text{obs}}$ increases the SNR which allows for
the detection of fainter signals, as well as signals of lower frequency $f$.
Furthermore, Figure 8 depicts a noticeable trend where the distribution of
recovered GBs shifts toward lower frequencies as observation time lengthens.
Specifically, at
$T_{\text{obs}}=$30\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$,
the majority of resolved GBs are centered around a frequency of approximately
$5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$, whereas by
$T_{\text{obs}}=$52\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$,
the frequency peak shifts to around
$4\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$. This shift occurs as the number of
recovered GBs converges towards the number of injected GBs with
$f>$5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$$, underscoring no further GBs at
this or higher frequencies are left to be resolved. This pattern highlights
the critical role of mission duration in the effective recovery of GBs,
especially at frequencies around $4\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$.

###  IV.3 Noise

In Fig. 9 we show the noise estimate obtained as described in Sec. III.1 for
the weekly analysis of the LDC2a data set. The GBs and MBHBs found in the
previous week are subtracted from the data. The noise estimate then consists
of an estimate for each channel
$n=\\{S_{A,residual},S_{E,residual},S_{T,residual}\\}$. For
$T_{\text{obs}}<$5\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
the noise curve is approximately linear for
$f\in[$0.3\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$,$10\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$]$
where the detectable GBs are expected. Since the GBs are quai-monochromatic
their signal is constantly recorded by LISA. On one hand, the galactic
confusion foreground is not as prominent yet and becomes louder by an increase
of $T_{\text{obs}}$ for
$f\in[$0.3\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$,$2.5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$]$.
On the other hand the signals with
$f\in[$2.5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$,$10\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$]$
are also getting louder, but in this frequency range there are fewer signals
which are overlapping. Therefore the GB extraction algorithm can analyze and
subtract the GBs with increasing accuracy with the increase of
$T_{\text{obs}}$.

It is also interesting to mention the increased noise at
$f\in[$0.2\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$,$0.5\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$]$
only for
$T_{\text{obs}}=$19\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$$. This
increase is due to the MBHB merging at
$T_{\text{obs}}\approx$19.06\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$$
which is not extracted yet but the inspiral phase is already well detectable.

Figure 9: Weekly noise estimate of the residual of the analysis of the LDC2a
data set. The outlier between 0.3 and 0.4 mHz is due to a MBHB merger in week
20 whose inspiral phase is already observable in week 19 (see text). The noise
estimate of each week can be accessed at [51].

###  IV.4 Computational Cost

In Fig. 10 we show the individual time and costs for each analysis as well as
the cumulative sum. As shown, our pipeline can keep the computational cost at
a neglectable level, which allows us to perform a weekly analysis with
simulated LISA data without spending much resources.

The jump from 7 to 8 weeks is due to a change in the search repetition for
each GB signal from 3 to 2, which lowers the computational cost. The jump from
week 13 to week 14 is due to the change in signal generation of FastGB. At
$T_{\text{obs}}<$3\text{\,}\text{months}$$ the implemented determination of
the number of frequency bins of a signal fails. Therefore we set the number of
frequency bins to 128 for all signals with
$T_{\text{obs}}<$3\text{\,}\text{months}$$, which appears to be more cost-
effective.

The CPU core times shown are for cost visualization only. The computation is
highly parallelizable since all frequency segments can be divided into two
groups of non-neighboring segments, each segment within a group can be
analyzed in parallel. This brings the time for the analysis down to minutes if
$\approx 3000$ CPU cores are available. For the interested reader the
parallelization is discussed in more detail in [29].

Figure 10: Computational CPU core time and the correlated cost in USD if run
on commercial hardware [56].

It is very unlikely that a new signal with
$f>$30\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$$ becomes detectable if the data
has already been analyzed for
$T_{\text{obs}}>$10\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$.
Thus, to further reduce the costs one could limit the search by an upper
frequency bound other than the maximum determined by $f_{\text{Nyquist}}$.
Furthermore, the number of searches for each GB can also be reduced from two
to only one to save computational costs.

##  V Robustness

Since the pipeline of this article has a low computational cost it is feasible
to conduct some robustness tests. First, we test the pipeline for one
iteration for the LDC2a data set with
$T_{\text{obs}}=$26\text{\,}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}\mathrm{s}$$
where we repeat the analysis ten times with different seeds. In Table 4 the
number of recovered and matched GBs for different seeds is presented. The
number of recovered and matched GBs have a standard deviation of 29 and 25 GBs
respectively. This is a rather low deviation compared to the mean of 5116 and
4342 GBs.

Furthermore, we matched the recovered signals of an analysis with the signals
of the other nine analyses. In Table 5 we present the number of times $N$ a
recovered signal is matched to the other analyses. For example, out of all the
$5086$ recovered signals of the analysis with seed $=1$ did $4184$ of them
also get extracted from all nine $N=9$ other analyses and only $N=47$ signals
were recovered from only one $N=1$ other analysis. One would expect that the
GBs with a higher SNR are recovered by all analyses, which is observed in Fig.
11 where most of the GBs with high amplitude are recovered by all or most of
the analyses. The signal strength depends on the amplitude and the inclination
angle. Therefore the signal strength does not perfectly correlate with the
amplitude but the correlation is still high enough for this observation.

Table 4: Results of recovering signals from the LDC2a dataset for $T_{\mathrm{obs}}=$0.5\text{\,}\mathrm{y}\mathrm{r}$$ with different seeds. First, the MBHBs are recovered, and then the GBs. A recovered signal is matched if it is paired with an injected signal and the scaled error is $\delta<0.3$. Seed | Recovered | Matched | Match rate
---|---|---|---
1 | $5086$ | $4311$ | $85\%$
2 | $5141$ | $4385$ | $85\%$
3 | $5156$ | $4361$ | $85\%$
4 | $5126$ | $4359$ | $85\%$
5 | $5120$ | $4327$ | $85\%$
6 | $5139$ | $4340$ | $84\%$
7 | $5117$ | $4364$ | $85\%$
8 | $5051$ | $4299$ | $85\%$
9 | $5095$ | $4327$ | $85\%$
10 | $5132$ | $4351$ | $85\%$
Mean | $5116.3$ | $4342.4$ | $85\%$
Std | $29$ | $25$ | $0.3\%$
Table 5: All signal for each seed gets potentially matched with a signal from each different seed. The table shows how many signals get matched with signals from all 9 other seeds. As an example for seed = 1, $4184$ signals got matched with a signal from each other 9 catalogs of recovered signals. Seed $N$ | 9 | 8 | 7 | 6 | 5 | 4 | 3 | 2 | 1
---|---|---|---|---|---|---|---|---|---
1 | 4184 | 334 | 146 | 82 | 54 | 60 | 57 | 44 | 47
2 | 4184 | 351 | 161 | 89 | 73 | 74 | 59 | 35 | 43
3 | 4187 | 346 | 152 | 83 | 71 | 81 | 60 | 44 | 46
4 | 4185 | 349 | 155 | 86 | 75 | 71 | 50 | 44 | 45
5 | 4184 | 351 | 164 | 75 | 75 | 69 | 52 | 37 | 35
6 | 4187 | 347 | 161 | 84 | 64 | 67 | 52 | 41 | 44
7 | 4187 | 352 | 153 | 75 | 60 | 65 | 52 | 36 | 50
8 | 4182 | 300 | 134 | 66 | 53 | 64 | 38 | 47 | 37
9 | 4186 | 341 | 150 | 86 | 59 | 53 | 52 | 28 | 54
10 | 4183 | 354 | 157 | 92 | 61 | 55 | 50 | 33 | 58
Mean | 4184.9 | 342.5 | 153.3 | 81.8 | 64.5 | 65.9 | 52.2 | 38.9 | 45.9
Std | 1.8 | 16.1 | 8.7 | 7.7 | 8.4 | 8.5 | 6.2 | 6.0 | 7.0
Table 6: The variables of interest include the count of detectable injected GB sources, the count of recovered sources, the count of matches with injected sources, and the match rate. The match rate is determined by dividing the number of matched signals by the total number of recovered signals. The overlap is included to get an evaluation comparable to other analyses [21, 22, 27, 57, 29]. Challenge |  $T_{\mathrm{obs}}$ [yr] | Recovered | $\delta<0.3$ | Match rate δ<0.3 | $\mathcal{O}>0.9$ | Match rate O>0.9
---|---|---|---|---|---|---
LDC2a | $0.5$ | $5086$ | $4311$ | $85\%$ | $4251$ | $84\%$
LDC2a without MBHBs | $0.5$ | $5099$ | $4326$ | $85\%$ | $4251$ | $83\%$
LDC2a | $1$ | $8313$ | $7192$ | $87\%$ | $7124$ | $86\%$
LDC2a without MBHBs | $1$ | $8323$ | $7217$ | $87\%$ | $7161$ | $86\%$
Figure 11:  Scatter plot of the recovered GBs of the analysis with seed $=1$.
The color represents the number of times $N$ this signal is recovered in the
other nine analyses.

In Table 6 we present the number of recovered and matched GBs of the LDC2a
data set where once the MBHBs are extracted by the pipeline and once without
the MBHBs in the data. As we can see the GB extraction does not change much,
which indicates a good quality MBHB extraction that is comparable with the
analysis as if there were no MBHBs in the first place.

##  VI Conclusion

This article introduces a novel method for extracting MBHBs and presents a
cost-effective global analysis of synthetic LISA data including MBHBs and GBs.
We expand upon the GB extraction algorithm [29] by incorporating the MBHB
extraction method, detailed in Sections III.2, and integrating both techniques
into a global fit pipeline shown in Fig. 2.

The algorithm undergoes testing on the LDC1-1 dataset, featuring a solitary
MBHB signal amidst instrument noise, and the LDC2a dataset, comprising 15
MBHBs, 30 million GBs, and instrument noise. The MBHB algorithm effectively
extracts the MBHB from the LDC1-1 dataset. Furthermore, the successful
extraction of self-injected signals, considering higher modes alongside
instrument noise, is presented in Sec. IV.1.

In the challenge of the LDC2a dataset, the MBHB algorithm extracts MBHBs such
that the GB extraction and noise estimation is without interference
effectively treating obtaining a residual as if the MBHBs were initially
absent from the dataset. Only if a strong MBHB signal merges right at the
boundary in the next week, the MBHB is visible in the noise estimate but still
not interfering with the GB extraction.

Thanks to the low cost of this pipeline, where analyzing
$T_{\text{obs}}=$1\text{\,}\mathrm{y}\mathrm{r}$$ costs only
$700\text{\,}\mathrm{C}\mathrm{P}\mathrm{U}\,\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\,\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}$,
corresponding to $20\text{\,}\mathrm{U}\mathrm{S}\mathrm{D}$, we can conduct
weekly analyses as new data is incorporated. Note that 52 weekly analyses
require a cumulative cost of only
$800\text{\,}\mathrm{U}\mathrm{S}\mathrm{D}$. Furthermore, the pipeline
operates on multiple frequency segments in parallel, thereby reducing the
analysis duration to mere minutes.

Starting from week 16, we observe a linear increase of 129 recovered GBs per
week with a constant match ratio of $85\%$ to $88\%$. Moreover, weekly
analyses reveal the correlation between $T_{\text{obs}}$ and the difference
between recovered and injected parameters.

Moreover, the noise estimation is effectively extended to arbitrary
$T_{\text{obs}}$. The weekly analysis unveils the evolving noise curve,
showcasing enhanced GB extraction for signals with
$f>$2\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$$ as $T_{\text{obs}}$ increases,
aligning the noise estimate of the residual closer to the instrument noise
level.

Furthermore, we leverage the low cost of the pipeline to conduct robustness
tests by repeating analyses of
$T_{\text{obs}}=$0.5\text{\,}\mathrm{y}\mathrm{r}$$ with different seeds.

The result is a robust, computationally efficient, high-quality GB and MBHB
signal extraction and noise estimation pipeline publicly available at [58,
59].

The next steps involve incorporating other GW signal types such as extreme
mass ratio inspirals, glitches, and data gaps. Additionally, the pipeline
could be extended to include pre-merger MBHB parameter estimation, as the
current pipeline focuses solely on post-merger MBHB signal extraction.

##  VII Acknowledgements

We are thankful to the LDC working group [16] for their efforts in
establishing and backing the LDC1-1 and LDC2a. Additionally, we appreciate the
GPU implementation of FASTLISARESPONSE [60] and the signal simulation of the
MBHBs BBHx [41, 10] which is available at [52]. We thank Franziska Riegger for
her helpful feedback. The computations to extract GBs were conducted on the
CPUs of the Euler cluster at ETH Zürich, and we extend our sincere
appreciation for this support. This project receives funding from the Swiss
National Science Foundation (SNF 200021_185051).

## References

  * Abbott  _et al._ [2019] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, _et al._ , Gwtc-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs, Physical Review X 9, 031040 (2019).
  * Colpi  _et al._ [2024] M. Colpi, K. Danzmann, M. Hewitson, K. Holley-Bockelmann, P. Jetzer, G. Nelemans, A. Petiteau, D. Shoemaker, C. Sopuerta, R. Stebbins, _et al._ , Lisa definition study report, arXiv preprint arXiv:2402.07571 (2024).
  * Digman and Cornish [2022] M. C. Digman and N. J. Cornish, Lisa gravitational wave sources in a time-varying galactic stochastic background, The Astrophysical Journal 940, 10 (2022).
  * Taam  _et al._ [1980] R. Taam, B. Flannery, and J. Faulkner, Gravitational radiation and the evolution of cataclysmic binaries, The Astrophysical Journal 239, 1017 (1980).
  * Willems  _et al._ [2008] B. Willems, A. Vecchio, and V. Kalogera, Probing white dwarf interiors with lisa: periastron precession in eccentric double white dwarfs, Physical review letters 100, 041102 (2008).
  * Nelemans  _et al._ [2010] G. Nelemans, L. Yungelson, M. v. d. Sluys, and C. A. Tout, The chemical composition of donors in am cvn stars and ultracompact x-ray binaries: observational tests of their formation, Monthly Notices of the Royal Astronomical Society 401, 1347 (2010).
  * Littenberg and Yunes [2019] T. B. Littenberg and N. Yunes, Binary white dwarfs as laboratories for extreme gravity with lisa, Classical and Quantum Gravity 36, 095017 (2019).
  * Piro [2019] A. L. Piro, Inferring the presence of tides in detached white dwarf binaries, The Astrophysical Journal Letters 885, L2 (2019).
  * Klein  _et al._ [2016] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder, F. Ohme, _et al._ , Science with the space-based interferometer elisa: Supermassive black hole binaries, Physical Review D 93, 024003 (2016).
  * Katz [2022] M. L. Katz, Fully automated end-to-end pipeline for massive black hole binary signal extraction from lisa data, Physical Review D 105, 044055 (2022).
  * Berti  _et al._ [2016] E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Spectroscopy of kerr black holes with earth-and space-based interferometers, Physical review letters 117, 101102 (2016).
  * Katz  _et al._ [2020a] M. L. Katz, L. Z. Kelley, F. Dosopoulou, S. Berry, L. Blecha, and S. L. Larson, Probing massive black hole binary populations with lisa, Monthly Notices of the Royal Astronomical Society 491, 2301 (2020a).
  * Petiteau  _et al._ [2011] A. Petiteau, S. Babak, and A. Sesana, Constraining the dark energy equation of state using lisa observations of spinning massive black hole binaries, The Astrophysical Journal 732, 82 (2011).
  * Gair  _et al._ [2013] J. R. Gair, M. Vallisneri, S. L. Larson, and J. G. Baker, Testing general relativity with low-frequency, space-based gravitational-wave detectors, Living Reviews in Relativity 16, 1 (2013).
  * Bogdanović [2014] T. Bogdanović, Supermassive black hole binaries: the search continues, in  _Gravitational Wave Astrophysics: Proceedings of the Third Session of the Sant Cugat Forum on Astrophysics_ (Springer, 2014) pp. 103–119.
  * [16] S. Babak and A. Petiteau, https://lisa-ldc.lal.in2p3.fr/.
  * Arnaud  _et al._ [2006] K. A. Arnaud, S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, C. Cutler, S. L. Larson, B. Sathyaprakash, M. Vallisneri, A. Vecchio, _et al._ , The mock lisa data challenges: an overview, arXiv preprint gr-qc/0609105 (2006).
  * Babak  _et al._ [2008] S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, J. Crowder, S. L. Larson, E. Plagnol, E. K. Porter, M. Vallisneri, A. Vecchio, _et al._ , The mock lisa data challenges: from challenge 1b to challenge 3, Classical and Quantum Gravity 25, 184026 (2008).
  * Babak  _et al._ [2010] S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, S. L. Larson, I. Mandel, S. T. McWilliams, A. Petiteau, E. K. Porter, E. L. Robinson, _et al._ , The mock lisa data challenges: from challenge 3 to challenge 4, Classical and Quantum Gravity 27, 084009 (2010).
  * Bouffanais and Porter [2016] Y. Bouffanais and E. K. Porter, Detecting compact galactic binaries using a hybrid swarm-based algorithm, Physical Review D 93, 064020 (2016).
  * Zhang  _et al._ [2021] X. Zhang, S. D. Mohanty, X. Zou, and Y. Liu, Resolving galactic binaries in lisa data using particle swarm optimization and cross-validation, arXiv preprint arXiv:2103.09391 (2021).
  * Gao  _et al._ [2023] P. Gao, X.-L. Fan, Z.-J. Cao, and X.-H. Zhang, Fast resolution of galactic binaries in lisa data, Physical Review D 107, 123029 (2023).
  * Crowder and Cornish [2007a] J. Crowder and N. J. Cornish, Solution to the galactic foreground problem for lisa, Phys. Rev. D 75, 043008 (2007a).
  * Crowder and Cornish [2007b] J. Crowder and N. J. Cornish, Extracting galactic binary signals from the first round of mock lisa data challenges, Classical and Quantum Gravity 24, S575 (2007b).
  * Littenberg [2011] T. B. Littenberg, Detection pipeline for galactic binaries in lisa data, Physical Review D 84, 063009 (2011).
  * Littenberg  _et al._ [2020] T. B. Littenberg, N. J. Cornish, K. Lackeos, and T. Robson, Global analysis of the gravitational wave signal from galactic binaries, Phys. Rev. D 101, 123021 (2020).
  * Littenberg and Cornish [2023] T. B. Littenberg and N. J. Cornish, Prototype global analysis of lisa data with multiple source types, arXiv preprint arXiv:2301.03673 (2023).
  * Strub  _et al._ [2022] S. H. Strub, L. Ferraioli, C. Schmelzbach, S. C. Stähler, and D. Giardini, Bayesian parameter estimation of galactic binaries in lisa data with gaussian process regression, Phys. Rev. D 106, 062003 (2022).
  * Strub  _et al._ [2023] S. H. Strub, L. Ferraioli, C. Schmelzbach, S. C. Stähler, and D. Giardini, Accelerating global parameter estimation of gravitational waves from galactic binaries using a genetic algorithm and gpus, Physical Review D 108, 103018 (2023).
  * Vecchio [2004] A. Vecchio, Lisa observations of rapidly spinning massive black hole binary systems, Physical Review D 70, 042001 (2004).
  * Berti  _et al._ [2005] E. Berti, A. Buonanno, and C. M. Will, Estimating spinning binary parameters and testing alternative theories of gravity with lisa, Physical Review D 71, 084025 (2005).
  * Arun [2006] K. Arun, Parameter estimation of coalescing supermassive black hole binaries with lisa, Physical Review D 74, 024025 (2006).
  * Cornish and Porter [2006] N. J. Cornish and E. K. Porter, Mcmc exploration of supermassive black hole binary inspirals, Classical and Quantum Gravity 23, S761 (2006).
  * Röver  _et al._ [2007] C. Röver, A. Stroeer, E. Bloomer, N. Christensen, J. Clark, M. Hendry, C. Messenger, R. Meyer, M. Pitkin, J. Toher, _et al._ , Inference on inspiral signals using lisa mldc data, Classical and Quantum Gravity 24, S521 (2007).
  * Feroz  _et al._ [2009] F. Feroz, J. R. Gair, M. P. Hobson, and E. K. Porter, Use of the multinest algorithm for gravitational wave data analysis, Classical and Quantum Gravity 26, 215003 (2009).
  * Gair and Porter [2009] J. R. Gair and E. K. Porter, Cosmic swarms: a search for supermassive black holes in the lisa data stream with a hybrid evolutionary algorithm, Classical and Quantum Gravity 26, 225004 (2009).
  * Petiteau  _et al._ [2010] A. Petiteau, Y. Shang, S. Babak, and F. Feroz, Search for spinning black hole binaries in mock lisa data using a genetic algorithm, Physical Review D 81, 104016 (2010).
  * Porter and Cornish [2015] E. K. Porter and N. J. Cornish, Fisher versus bayes: A comparison of parameter estimation techniques for massive black hole binaries to high redshifts with elisa, Physical Review D 91, 104001 (2015).
  * Cornish and Shuman [2020] N. J. Cornish and K. Shuman, Black hole hunting with lisa, Physical Review D 101, 124008 (2020).
  * Marsat  _et al._ [2021] S. Marsat, J. G. Baker, and T. D. Canton, Exploring the bayesian parameter estimation of binary black holes with lisa, Physical Review D 103, 083011 (2021).
  * Katz  _et al._ [2020b] M. L. Katz, S. Marsat, A. J. K. Chua, S. Babak, and S. L. Larson, Gpu-accelerated massive black hole binary parameter estimation with lisa, Physical Review D 102, 023033 (2020b).
  * Cornish [2022] N. J. Cornish, Low latency detection of massive black hole binaries, Physical Review D 105, 044007 (2022).
  * Dey  _et al._ [2021] K. Dey, N. Karnesis, A. Toubiana, E. Barausse, N. Korsakova, Q. Baghi, and S. Basak, Effect of data gaps on the detectability and parameter estimation of massive black hole binaries with lisa, Physical Review D 104, 044035 (2021).
  * Tinto and Armstrong [1999] M. Tinto and J. W. Armstrong, Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation, Physical Review D 59, 102003 (1999).
  * Armstrong  _et al._ [1999] J. W. Armstrong, F. B. Estabrook, and M. Tinto, Time-delay interferometry for space-based gravitational wave searches, The Astrophysical Journal 527, 814 (1999).
  * Estabrook  _et al._ [2000] F. B. Estabrook, M. Tinto, and J. W. Armstrong, Time-delay analysis of lisa gravitational wave data: Elimination of spacecraft motion effects, Physical Review D 62, 042002 (2000).
  * Dhurandhar  _et al._ [2002] S. V. Dhurandhar, K. R. Nayak, and J.-Y. Vinet, Algebraic approach to time-delay data analysis for lisa, Physical Review D 65, 102002 (2002).
  * Tinto and Dhurandhar [2014] M. Tinto and S. V. Dhurandhar, Time-delay interferometry, Living Reviews in Relativity 17, 1 (2014).
  * Prince  _et al._ [2002] T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, Lisa optimal sensitivity, Physical Review D 66, 122002 (2002).
  * Vallisneri [2005] M. Vallisneri, Synthetic lisa: Simulating time delay interferometry in a model lisa, Physical Review D 71, 022001 (2005).
  * Strub [2024a] S. Strub, https://doi.org/10.5281/zenodo.10852821 (2024a).
  * Katz [2021] M. Katz, mikekatz04/BBHx: First official public release, 10.5281/zenodo.5730688 (2021).
  * Storn and Price [1997] R. Storn and K. Price, Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces, Journal of global optimization 11, 341 (1997).
  * Cornish and Littenberg [2007] N. J. Cornish and T. B. Littenberg, Tests of bayesian model selection techniques for gravitational wave astronomy, Physical Review D 76, 083006 (2007).
  * Katz  _et al._ [2022] M. L. Katz, J.-B. Bayle, A. J. K. Chua, and M. Vallisneri, Assessing the data-analysis impact of lisa orbit approximations using a gpu-accelerated response model, Physical Review D 106, 103001 (2022).
  * [56] Google-Cloud, https://cloud.google.com/pricing/, accessed on 2023-06-12.
  * Lackeos  _et al._ [2023] K. Lackeos, T. B. Littenberg, N. J. Cornish, and J. I. Thorpe, The lisa data challenge radler analysis and time-dependent ultra-compact binary catalogues, Astronomy & Astrophysics 678, A123 (2023).
  * Strub [2024b] S. Strub, https://doi.org/10.5281/zenodo.10853529 (2024b).
  * Strub [2024c] S. Strub, https://doi.org/10.5281/zenodo.10853495 (2024c).
  * Katz and Bayle [2022] M. Katz and J.-B. Bayle, mikekatz04/lisa-on-gpu: First official public release!, 10.5281/zenodo.5867731 (2022).


# Generalized hypergeometric Bernoulli numbers

Kalyan Chakraborty  Harish-Chandra Research Institute, Allahabad, 211019,
Uttar Pradesh, India
Kerala school of Mathematics, Kunnamangalam, Kozhikode, 673571, Kerala, India
<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Takao Komatsu  Department of
Mathematical Sciences, School of Science, Zhejiang Sci-Tech University,
Hangzhou 310018, China<EMAIL_ADDRESS>
###### Abstract.

We introduce generalized hypergeometric Bernoulli numbers for Dirichlet
characters. We study their properties, including relations, expressions and
determinants. At the end in Appendix we derive first few expressions of these
numbers.

###### Key words and phrases:

Bernoulli numbers; hypergeometric Bernoulli numbers; hypergeometric functions;
Dirichlet characters

###### 1991 Mathematics Subject Classification:

Primary 11B68; Secondary 11B37, 33C15.

##  1\. Introduction

Let ${}_{1}F_{1}(a;b;z)$ denote the confluent hypergeometric function defined
by

${}_{1}F_{1}(a;b;z)=\sum_{n=0}^{\infty}\frac{(a)^{(n)}}{(b)^{(n)}}\frac{z^{n}}{n!}$

with the rising factorial $(x)^{(n)}=x(x+1)\dots(x+n-1)$ ($n\geq 1$) and
$(x)^{(0)}=1$. Also for $N\geq 1$, hypergeometric Bernoulli numbers $B_{N,n}$
([2, 6, 7, 8, 9, 13, 18]) are given by

(1.1)
$\frac{1}{{}_{1}F_{1}(1;N+1;t)}=\frac{t^{N}/N!}{e^{t}-\sum_{n=0}^{N-1}t^{n}/n!}=\sum_{n=0}^{\infty}B_{N,n}\frac{t^{n}}{n!}$

with $B_{1,n}=B_{n}$, the classical Bernoulli numbers, defined by

$\frac{t}{e^{t}-1}=\sum_{n=0}^{\infty}B_{n}\frac{t^{n}}{n!}\,.$

In [10] a result of Kamano on a multiple binomial sum is generalized, higher
order hypergeometric Bernoulli polynomials are defined, and its differential
equation and recursion formulas are studied among other properties. It is
well-known that some of the known properties of the Bernoulli numbers can be
obtained from fundamental relationships between complete and elementary
symmetric functions. In [11] two new closed form representations for the
Apostol-Bernoulli polynomials are derived as natural generalizations of those
of the Bernoulli polynomials. In [17] an infinite family of relationships
between complete and elementary symmetric functions is introduced. From these
relationships, the author deduces connections between some lacunary recurrence
relations with gaps of length $2r$ for the Bernoulli numbers and the integer
partitions into at most $r-1$ parts. In [20] new bounds for the ratio of two
adjacent even-indexed Bernoulli numbers are obtained, solving Qi’s conjecture
on the related topic, getting a tighter lower bound for this ratio is
obtained, and moreover, some conjectures on related topics are proposed.

In [14], hypergeometric Cauchy numbers are defined and studied as analog of
the hypergeometric Bernoulli numbers. This concept is further generalized in
[15] by introducing the general hypergeometric Cauchy polynomials. In [16],
the classical Euler numbers are generalized to define the higher order
hypergeometric Euler numbers. Recently, in [12], by using the Hasse-
Teichmüller derivatives, two explicit expressions for the related numbers of
higher order Appell polynomials are obtained. One of them presents a
determinant expression for the related numbers of higher order Appell
polynomials, which involves several determinant expressions of special
numbers, such as the higher order generalized hypergeometric Bernoulli and
Cauchy numbers

In this paper, we introduce generalized hypergeometric Bernoulli numbers for
Dirichlet characters. Let $\chi$ be a Dirichlet character with conductor $f$
and $N$ be a positive integer. We define the generalized hypergeometric
Bernoulli numbers $B_{N,n,\chi}$ by

(1.2) $\sum_{a=1}^{f}\frac{\chi(a)\
t^{N}e^{at}/N!}{e^{ft}-\sum_{n=0}^{N-1}\frac{(ft)^{n}}{n!}}=\sum_{n=0}^{\infty}B_{N,n,\chi}\frac{t^{n}}{n!}\quad(|t|<2\pi/f)\,.$

One recovers the classical generalized Bernoulli number (see, e.g., [1, 3, 4,
19]) $B_{n,\chi}:=B_{1,n,\chi}$, which are simply Bernoulli numbers twisted by
a Dirichlet character. As seen in [3, Chapter 4], generalized Bernoulli
numbers and Bernoulli polynomials are related. A closed form expression for
sums of products of any number of generalized Bernoulli numbers is given in
[4] and shortened recurrence relations for generalized Bernoulli numbers and
polynomials attached to a primitive Dirichlet character is being studied in
[1].

##  2\. Basic properties

We begin with some basic expressions of generalized hypergeometric Bernoulli
numbers.

######  Proposition 1.

Let $\chi={\bf 1}$ be the trivial character of conductor $1$. For any $N\geq
1$ and $n\geq 0$, put $N(n):=$min$\\{N-1,n\\}$. Then we have

$B_{N,n,{\bf
1}}=\begin{cases}\displaystyle{\sum_{m=0}^{N(n)}\binom{n}{m}B_{N,n-m}}&\text{if}\
n\neq N,\\\ \\\
\displaystyle{1+\sum_{m=0}^{N-1}\binom{N}{m}B_{N,N-m}}&\text{if}\
n=N.\end{cases}$

Especially, for $N=1$, we have

$B_{1,n,{\bf 1}}=\begin{cases}B_{n}&\text{if}\ n\neq 1,\\\ \\\ 1/2&\text{if}\
n=1.\end{cases}$

###### Proof.

By the definition of generalized hypergeometric Bernoulli numbers, we have

$\displaystyle\sum_{n=0}^{\infty}B_{N,n,{\bf 1}}\frac{t^{n}}{n!}$
$\displaystyle=\frac{t^{N}e^{t}/N!}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}$
$\displaystyle=\frac{t^{N}/N!(e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!})+t^{N}/N!\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}$
$\displaystyle=\frac{t^{N}}{N!}+\sum_{n=0}^{N-1}\frac{t^{n}}{n!}\times\sum_{n=0}^{\infty}B_{N,n}\frac{t^{n}}{n!}$
$\displaystyle=\frac{t^{N}}{N!}+\sum_{m=0}^{N-1}\sum_{n=0}^{\infty}B_{N,n}\frac{t^{m+n}}{m!n!}$
$\displaystyle=\frac{t^{N}}{N!}+\sum_{\ell=0}^{N-1}\sum_{m=0}^{\ell}B_{N,\ell-m}\frac{t^{\ell}}{m!(\ell-m)!}+\sum_{\ell=N}^{\infty}\sum_{m=0}^{N-1}B_{N,\ell-m}\frac{t^{\ell}}{m!(\ell-m)!}$
$\displaystyle=\sum_{\ell=0}^{N-1}\sum_{m=0}^{\ell}\binom{\ell}{m}B_{N,\ell-m}\frac{t^{\ell}}{\ell!}+\left(1+\sum_{m=0}^{N-1}\binom{N}{m}B_{N,N-m}\right)\frac{t^{N}}{N!}$
$\displaystyle+\sum_{\ell=N+1}^{\infty}\sum_{m=0}^{N-1}\binom{\ell}{m}B_{N,\ell-m}\frac{t^{\ell}}{\ell!},$

and we get the assertion. ∎

Howard [8, 9] defined the hypergeometric Bernoulli polynomials
$B_{N,n}(x)\in{\mathbb{Q}}[x]$ by

(2.1)
$\frac{t^{N}e^{xt}/N!}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}=\sum_{n=0}^{\infty}B_{N,n}(x)\frac{t^{n}}{n!}.$

Note that $B_{n}(x):=B_{1,n}(x)$ be the classical Bernoulli polynomial and
$B_{N,n}(0)=B_{N,n}$. For classical case, the following relation between the
generalized Bernoulli numbers and the Bernoulli polynomials is well-known [19,
Proposition 4.1]:

$B_{n,\chi}=f^{n-1}\sum_{a=1}^{f}\chi(a)B_{n}\left(\frac{a}{f}\right).$

The following proposition is a generalization of the above classical case.

######  Proposition 2.

Let $\chi$ be a Dirichlet character of conductor $f$. For any $N\geq 1$ and
$n\geq 0$,

$B_{N,n,\chi}=f^{n-N}\sum_{a=1}^{f}\chi(a)B_{N,n}\left(\frac{a}{f}\right).$

###### Proof.

By the definition (of hypergeometric Bernoulli polynomial), we have,

$\sum_{n=0}^{\infty}B_{N,n}\left(\frac{a}{f}\right)\frac{(ft)^{n}}{n!}=\frac{(ft)^{N}e^{a/f\times
ft}/N!}{e^{ft}-\sum_{n=0}^{N-1}\frac{(ft)^{n}}{n!}}.$

Now multiplying by $\chi(a)$ and taking a sum on $a$ from $1$ to $f$, we have

$\displaystyle
f^{n}\sum_{n=0}^{\infty}\sum_{a=1}^{f}\chi(a)B_{N,n}\left(\frac{a}{f}\right)\frac{t^{n}}{n!}$
$\displaystyle=f^{N}\sum_{a=1}^{f}\frac{\chi(a)t^{N}e^{at}/N!}{e^{ft}-\sum_{n=0}^{N-1}\frac{(ft)^{n}}{n!}}$
$\displaystyle=f^{N}\sum_{n=0}^{\infty}B_{N,n,\chi}\frac{t^{n}}{n!}.$

Therefore,

$B_{N,n,\chi}=f^{n-N}\sum_{a=1}^{f}\chi(a)B_{N,n}\left(\frac{a}{f}\right)$

for any $n\geq 0$. ∎

Consequently, we can express the generalized hypergeometric Bernoulli numbers
in terms of hypergeometric Bernoulli numbers

######  Corollary 1.

For any integers $N\geq 1$ and $n\geq 0$,

$B_{N,n,\chi}=\sum_{a=1}^{f}\chi(a)\sum_{k=0}^{n}\binom{n}{k}B_{N,k}a^{n-k}f^{k-N}.$

###### Proof.

Since

$\displaystyle\sum_{n=0}^{\infty}B_{N,n}(x)\frac{t^{n}}{n!}$
$\displaystyle=\frac{t^{N}e^{xt}/N!}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}$
$\displaystyle=\left(\sum_{n=0}^{\infty}B_{N,n}\frac{t^{n}}{n!}\right)\left(\sum_{m=0}^{\infty}x^{m}\frac{t^{m}}{m!}\right)$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\binom{n}{k}B_{N,k}x^{n-k}\frac{t^{n}}{n!},$

comparing the coefficients on both sides, we get

$B_{N,n}(x)=\sum_{k=0}^{n}\binom{n}{k}B_{N,k}x^{n-k}.$

By Proposition 2 with $x=a/f$, we obtain the desired result. ∎

##  3\. Fundamental properties of the generalized hypergeometric Bernoulli
polynomials

We also define the generalized hypergeometric Bernoulli polynomials
$B_{N,n,\chi}$ by

(3.1) $\sum_{a=1}^{f}\frac{\chi(a)\
t^{N}e^{(x+a)t}/N!}{e^{ft}-\sum_{n=0}^{N-1}\frac{(ft)^{n}}{n!}}=\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}\quad(|t|<2\pi/f)\,.$

When $N=1$ in (3.1), $B_{n,\chi}(x):=B_{1,n,\chi}(x)$ is the classical
generalized Bernoulli polynomial. When $\chi={\bf 1}$,
$B_{N,n}(x):=B_{N,n,{\bf 1}}(x)$ is the hypergeometric Bernoulli polynomial in
(2.1).

######  Proposition 3.

Let $\chi={\bf 1}$ be the trivial character of conductor $1$. For any $N\geq
1$ and $n\geq 0$, we have

$B_{N,n,{\bf
1}}(x)=\begin{cases}\displaystyle{\sum_{i=0}^{n}\binom{n}{i}B_{N,n-i}}(x)&\text{if}\
0\leq n\leq N-1,\\\ \\\
\displaystyle{\binom{n}{N}x^{n-N}+\sum_{i=0}^{N-1}\binom{n}{i}B_{N,n-i}}(x)&\text{if}\
n\geq N.\end{cases}$

Especially, for $N=1$, we have

$B_{1,n,{\bf 1}}(x)=\begin{cases}1&\text{if}\ n=0,\\\
nx^{n-1}+B_{n}(x)&\text{if}\ n\geq 1.\end{cases}$

###### Proof.

By the definition in (3.1), we have

$\displaystyle\sum_{n=0}^{\infty}B_{N,n,{\bf
1}}(x)\frac{t^{n}}{n!}=\frac{t^{N}e^{(x+1)t}/N!}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}$
$\displaystyle=\frac{t^{N}/N!\bigl{(}e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}\bigr{)}e^{xt}+t^{N}/N!e^{xt}\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}{e^{t}-\sum_{n=0}^{N-1}\frac{t^{n}}{n!}}$
$\displaystyle=\frac{t^{N}}{N!}e^{xt}+\left(\sum_{n=0}^{N-1}\frac{t^{n}}{n!}\right)\left(\sum_{n=0}^{\infty}B_{N,n}(x)\frac{t^{n}}{n!}\right)$
$\displaystyle=\sum_{n=0}^{\infty}\frac{x^{n}}{N!n!}t^{N+n}+\sum_{n=0}^{\infty}\sum_{i=0}^{\min\\{n,N-1\\}}\binom{n}{i}B_{N,n-i}(x)\frac{t^{n}}{n!}$
$\displaystyle=\sum_{n=N}^{\infty}\binom{n}{N}x^{n-N}\frac{t^{n}}{n!}+\sum_{n=0}^{N-1}\sum_{i=0}^{n}\binom{n}{i}B_{N,n-i}(x)\frac{t^{n}}{n!}$
$\displaystyle\quad+\sum_{n=N}^{\infty}\sum_{i=0}^{N-1}\binom{n}{i}B_{N,n-i}(x)\frac{t^{n}}{n!}\,.$

Comparing the coefficients on both sides, we get the result. ∎

######  Proposition 4.

For $n\geq 0$,

$B_{N,n,\chi}(x+y)=\sum_{k=0}^{n}\binom{n}{k}B_{N,k,\chi}(x)y^{n-k}\,.$

###### Proof.

By the definition in (3.1), we have

$\displaystyle\sum_{n=0}^{\infty}B_{N,n,\chi}(x+y)\frac{t^{n}}{n!}$
$\displaystyle=\left(\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}\right)e^{yt}$
$\displaystyle=\left(\sum_{k=0}^{\infty}B_{N,k,\chi}(x)\frac{t^{k}}{k!}\right)\left(\sum_{l=0}^{\infty}\frac{y^{l}t^{l}}{l!}\right)$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\binom{n}{k}B_{N,k,\chi}(x)y^{n-k}\frac{t^{n}}{n!}\,.$

Comparing the coefficients on both sides, we get the result. ∎

When $x=0$ and $y$ is replaced by $x$ in Proposition 2, we have a relation
between generalized hypergeometric Bernoulli polynomials and numbers.

######  Corollary 2.

For $n\geq 0$,

$B_{N,n,\chi}(x)=\sum_{k=0}^{n}\binom{n}{k}B_{N,k,\chi}x^{n-k}\,.$

There is another relation of generalized hypergeometric Bernoulli polynomials
with their numbers in terms of Stirling numbers of the second kind
$\left\\{n\atop k\right\\}$, whose generating function is given by

$\frac{(e^{t}-1)^{\ell}}{\ell!}=\sum_{n=\ell}^{\infty}\left\\{n\atop\ell\right\\}\frac{t^{n}}{n!}\,.$

######  Proposition 5.

For $n\geq 0$,

$B_{N,n,\chi}(x)=\sum_{k=0}^{n}\sum_{l=0}^{k}\binom{n}{k}\left\\{k\atop
l\right\\}(x)_{l}B_{N,n-k,\chi}\,,$

where $(x)_{l}=x(x-1)\cdots(x-l+1)$ ($l\geq 1$) denotes the falling factorial
with $(x)_{0}=1$.

###### Proof.

$\displaystyle\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}$
$\displaystyle=\left(\sum_{n=0}^{\infty}B_{N,n,\chi}\frac{t^{n}}{n!}\right)\bigl{(}(e^{t}-1)+1\bigr{)}^{x}$
$\displaystyle=\left(\sum_{n=0}^{\infty}B_{N,n,\chi}\frac{t^{n}}{n!}\right)\left(\sum_{l=0}^{\infty}\binom{x}{l}l!\sum_{n=l}^{\infty}\left\\{n\atop
l\right\\}\frac{t^{n}}{n!}\right)$
$\displaystyle=\left(\sum_{n=0}^{\infty}B_{N,n,\chi}\frac{t^{n}}{n!}\right)\left(\sum_{n=0}^{\infty}\sum_{l=0}^{n}\left\\{n\atop
l\right\\}(x)_{l}\frac{t^{n}}{n!}\right)$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\sum_{l=0}^{k}\binom{n}{k}\left\\{k\atop
l\right\\}(x)_{l}B_{N,n-k,\chi}\frac{t^{n}}{n!}\,.$

Comparing the coefficients on both sides, we get the result. ∎

We can see that the generalized hypergeometric Bernoulli numbers are Appell.

######  Proposition 6.

For $n\geq 1$,

$\frac{d}{dx}B_{N,n,\chi}(x)=nB_{N,n-1,\chi}(x)\,.$

###### Proof.

From the definition in (3.1),

$\displaystyle\frac{d}{dx}\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}$
$\displaystyle=t\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}$
$\displaystyle=\sum_{n=0}^{\infty}(n+1)B_{N,n,\chi}(x)\frac{t^{n+1}}{(n+1)!}$
$\displaystyle=\sum_{n=1}^{\infty}nB_{N,n-1,\chi}(x)\frac{t^{n}}{n!}\,.$

Comparing the coefficients on both sides, we get the result. ∎

##  4\. Recurrence relations

From now on, for convenience sake let us put

$S_{n}:=\sum_{a=1}^{f}\chi(a)a^{n}\quad\hbox{and}\quad
S_{n}(x):=\sum_{a=1}^{f}\chi(a)(x+a)^{n}\,.$

The following are recurrence relations of the generalized hypergeometric
Bernoulli numbers and polynomials with Dirichlet character.

######  Proposition 7.

For $n\geq 1$,

(4.1)
$\displaystyle\binom{N+n}{n}^{-1}\sum_{i=0}^{n}\binom{N+n}{i}B_{N,i,\chi}f^{N+n-i}=S_{n}\,,$
(4.2)
$\displaystyle\binom{N+n}{n}^{-1}\sum_{i=0}^{n}\binom{N+n}{i}B_{N,i,\chi}(x)f^{N+n-i}=S_{n}(x)\,.$

Remark. When $N=1$ in Proposition 7, we have

$\displaystyle\sum_{i=0}^{n-1}\binom{n}{i}\frac{B_{i+1,\chi}}{i+1}f^{n-i}+\frac{B_{0,\chi}}{n+1}f^{n+1}=S_{n}\,,$
$\displaystyle\sum_{i=0}^{n-1}\binom{n}{i}\frac{B_{i+1,\chi}(x)}{i+1}f^{n-i}+\frac{B_{0,\chi}(x)}{n+1}f^{n+1}=S_{n}(x)\,.$

which are [1, (1.12),(1.13)]. When $N=f=1$ in Proposition 7, we have the well-
known basic recurrence relations

$\displaystyle\sum_{i=0}^{n-1}\binom{n}{i}\frac{B_{i+1}}{i+1}+\frac{1}{n+1}=0\,;$
$\displaystyle\sum_{i=0}^{n-1}\binom{n}{i}\frac{B_{i+1}(x)}{i+1}+\frac{1}{n+1}=x^{n}\,.$

When $\chi=\mathbf{1}$ is the trivial character of conductor $1$, the relation
(4.1) is reduced to

$\sum_{i=0}^{n}\binom{N+n}{i}B_{N,i,{\mathbf{1}}}=0\,.$

which is [2, Proposition 1 (3)].

###### Proof of Proposition 7..

We shall prove the relation (4.2). The relation (4.1) is the special case when
$x=0$. From the definition in (3.1),

$\sum_{a=1}^{f}\chi(a)\frac{t^{N}}{N!}e^{(x+a)t}=\left(e^{ft}-\sum_{n=0}^{N-1}\frac{(ft)^{n}}{n!}\right)\left(\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}\right)\,.$

The right-hand side is equal to

$\left(\sum_{n=0}^{\infty}\frac{(ft)^{N+n}}{(N+n)!}\right)\left(\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}\right)=\sum_{n=0}^{\infty}\sum_{i=0}^{n}\frac{B_{N,i,\chi}(x)}{i!}\frac{f^{N+n-i}}{(N+n-i)!}t^{N+n}\,.$

The left-hand side is equal to

$\sum_{n=0}^{\infty}\sum_{a=1}^{f}\chi(a)\frac{t^{N}}{N!}(x+a)^{n}\frac{t^{n}}{n!}=\sum_{n=0}^{\infty}\sum_{a=1}^{f}\chi(a)(x+a)^{n}\frac{t^{N+n}}{N!n!}\,.$

Comparing the coefficients on both sides, we get

$\sum_{i=0}^{n}\frac{B_{N,i,\chi}(x)}{i!}\frac{f^{N+n-i}}{(N+n-i)!}=\sum_{a=1}^{f}\frac{\chi(a)(x+a)^{n}}{N!n!}\,,$

yielding the relation (4.2). ∎

##  5\. Expressions

We give expressions of the generalized hypergeometric Bernoulli polynomials
and numbers.

######  Theorem 1.

For $n\geq 1$,

(5.1) $\displaystyle B_{N,n,\chi}$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}T_{r}(k)S_{n-k}\,,$
(5.2) $\displaystyle B_{N,n,\chi}(x)$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}T_{r}(k)S_{n-k}(x)\,.$

where

$T_{r}(k):=(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=k\atop i_{1},\dots,i_{r}\geq
1}\frac{1}{(N+i_{1})!\cdots(N+i_{r})!}$

with $T_{0}(0)=1$ and $T_{0}(k)=0$ ($k\geq 1$).

Remark. When $\chi=\mathbf{1}$ is the trivial character of conductor $1$,
expression (5.1) reduces to

$B_{N,n,{\mathbf{1}}}=n!\sum_{r=1}^{n}(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=n\atop
i_{1},\dots,i_{r}\geq 1}\frac{1}{(N+i_{1})!\cdots(N+i_{r})!}\,,$

which is [2, Proposition 1 (4)].

###### Proof of Theorem 1..

We shall prove (5.2) as (5.1) can be easily deduced. From the definition in
(3.1),

$\displaystyle\sum_{n=0}^{\infty}B_{N,n,\chi}(x)\frac{t^{n}}{n!}$
$\displaystyle=\left(\sum_{a=1}^{f}\chi(a)\sum_{n=0}^{\infty}(x+a)^{n}\frac{t^{n}}{n!}\right)\frac{1}{f^{N}}\left(\sum_{n=0}^{\infty}\frac{N!f^{n}}{(N+n)!}t^{n}\right)^{-1}$
$\displaystyle=\frac{1}{f^{N}}\left(\sum_{n=0}^{\infty}S_{n}(x)\frac{t^{n}}{n!}\right)\sum_{r=0}^{\infty}\left(-\sum_{n=1}^{\infty}\frac{N!f^{n}}{(N+n)!}t^{n}\right)^{r}$
$\displaystyle=\frac{1}{f^{N}}\left(\sum_{n=0}^{\infty}S_{n}(x)\frac{t^{n}}{n!}\right)\sum_{r=0}^{\infty}(-1)^{r}\sum_{n=r}^{\infty}\sum_{i_{1}+\cdots+i_{r}=n\atop
i_{1},\dots,i_{r}\geq 1}\frac{(N!)^{r}f^{n}}{(N+i_{1})!\cdots(N+i_{r})!}t^{n}$
$\displaystyle=\frac{1}{f^{N}}\left(\sum_{n=0}^{\infty}S_{n}(x)\frac{t^{n}}{n!}\right)\sum_{n=0}^{\infty}\sum_{r=0}^{n}(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=n\atop
i_{1},\dots,i_{r}\geq
1}\frac{n!f^{n}}{(N+i_{1})!\cdots(N+i_{r})!}\frac{t^{n}}{n!}$
$\displaystyle=\frac{1}{f^{N}}\sum_{n=0}^{\infty}\sum_{k=0}^{n}\binom{n}{k}S_{n-k}(x)\sum_{r=0}^{k}(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=k\atop
i_{1},\dots,i_{r}\geq
1}\frac{k!f^{k}}{(N+i_{1})!\cdots(N+i_{r})!}\frac{t^{n}}{n!}\,.$

Comparing the coefficients on both sides, we get (5.2). Note that the
summation is empty when $k\geq 1$ and $r=0$. ∎

There are alternative expressions with additional binomial coefficients.

######  Theorem 2.

For $n\geq 1$,

(5.3) $\displaystyle B_{N,n,\chi}$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}\binom{k+1}{r+1}\widetilde{T}_{r}(k)S_{n-k}\,,$
(5.4) $\displaystyle B_{N,n,\chi}(x)$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}\binom{k+1}{r+1}\widetilde{T}_{r}(k)S_{n-k}(x)\,.$

where

$\widetilde{T}_{r}(k):=(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=k\atop
i_{1},\dots,i_{r}\geq 0}\frac{1}{(N+i_{1})!\cdots(N+i_{r})!}$

with $\widetilde{T}_{0}(0)=1$ and $\widetilde{T}_{0}(k)=0$ ($k\geq 1$).

Remark. When $\chi=\mathbf{1}$ is the trivial character of conductor $1$,
(5.3) is reduced to

$B_{N,n,{\mathbf{1}}}=n!\sum_{r=1}^{n}\binom{n+1}{r+1}(-N!)^{r}\sum_{i_{1}+\cdots+i_{r}=n\atop
i_{1},\dots,i_{r}\geq 0}\frac{1}{(N+i_{1})!\cdots(N+i_{r})!}\,,$

which is [2, Proposition 1 (5)].

###### Proof of Theorem 2..

We shall prove the expression (5.4). Put

$w:=\sum_{n=1}^{\infty}\frac{N!f^{n}}{(N+n)!}t^{n}\,.$

Then, we have

$\displaystyle B_{N,n,\chi}(x)$
$\displaystyle:=\left.\frac{d^{n}}{dt^{n}}\frac{1}{f^{N}}(1+w)^{-1}\sum_{a=1}^{f}\chi(a)e^{(x+a)t}\right|_{t=0}$
$\displaystyle=\frac{1}{f^{N}}\sum_{k=0}^{n}\binom{n}{k}\left.\frac{d^{k}}{dt^{k}}\sum_{l=0}^{\infty}(-w)^{l}\right|_{t=0}\left.\frac{d^{n-k}}{dt^{n-k}}\sum_{\ell=0}^{\infty}S_{\ell}(x)\frac{t^{\ell}}{\ell!}\right|_{t=0}$
$\displaystyle=\frac{1}{f^{N}}\sum_{k=0}^{n}\binom{n}{k}\sum_{l=0}^{k}\sum_{r=0}^{l}(-1)^{r}\binom{l}{r}$
$\displaystyle\qquad\times\left.\frac{d^{k}}{dt^{k}}\left(\sum_{m=0}^{\infty}\frac{N!f^{m}}{(N+m)!}t^{m}\right)^{r}\right|_{t=0}\left.\sum_{\ell=0}^{\infty}S_{n-k+\ell}(x)\frac{t^{\ell}}{\ell!}\right|_{t=0}$
$\displaystyle=\frac{1}{f^{N}}\sum_{k=0}^{n}\binom{n}{k}\sum_{l=0}^{k}\sum_{r=0}^{l}(-1)^{r}\binom{l}{r}$
$\displaystyle\qquad\times\left.\frac{d^{k}}{dt^{k}}\sum_{m=0}^{\infty}\sum_{i_{1}+\cdots+i_{r}=m\atop
i_{1},\dots,i_{r}\geq
0}\frac{(N!)^{r}f^{m}}{(N+i_{1})!\cdots(N+i_{r})!}t^{m}\right|_{t=0}S_{n-k}(x)$
$\displaystyle=\frac{1}{f^{N}}\sum_{k=0}^{n}\binom{n}{k}\sum_{l=0}^{k}\sum_{r=0}^{l}\binom{l}{r}\left.\sum_{m=k}^{\infty}\frac{m!f^{m}}{(m-k)!}\widetilde{T}_{r}(m)t^{m-k}\right|_{t=0}S_{n-k}(x)$
$\displaystyle=\frac{1}{f^{N}}\sum_{k=0}^{n}\binom{n}{k}\sum_{l=0}^{k}\sum_{r=0}^{l}\binom{l}{r}k!f^{k}\widetilde{T}_{r}(k)S_{n-k}(x)$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}\left(\sum_{l=r}^{k}\binom{l}{r}\right)\widetilde{T}_{r}(k)S_{n-k}(x)$
$\displaystyle=\sum_{k=0}^{n}\frac{k!}{f^{N-k}}\binom{n}{k}\sum_{r=0}^{k}\binom{k+1}{r+1}\widetilde{T}_{r}(k)S_{n-k}(x)\,.$

∎

##  6\. Determinants

Recurrence formula from Proposition 7

(6.1)
$B_{N,n,\chi}=\frac{S_{n}}{f^{N}}-\sum_{i=0}^{n-1}\frac{N!n!f^{n-i}}{(N+n-i)!i!}B_{N,i,\chi}$

assists us to derive determinant expressions of generalized hypergeometric
Bernoulli numbers and polynomials. Such determinant expressions were studied
by Glaisher ([5]) extensively for Bernoulli, Cauchy, Euler and more numbers.

######  Theorem 3.

For $n\geq 1$,

(6.7) $\displaystyle B_{N,n,\chi}$
$\displaystyle=(-1)^{n}n!\left|\begin{array}[]{ccccc}\frac{N!}{(N+1)!}f&1&0&&\\\
\frac{N!}{(N+2)!}f^{2}&\frac{N!}{(N+1)!}f&&&\\\ \vdots&&&&0\\\
\frac{N!}{(N+n-1)!}f^{n-1}&\cdots&\frac{N!}{(N+2)!}f^{2}&\frac{N!}{(N+1)!}f&1\\\
\widehat{S}_{n}&\widehat{S}_{n-1}&\cdots&\widehat{S}_{2}&\widehat{S}_{1}\end{array}\right|\,$

where

$\widehat{S}_{n}=\frac{1}{f^{N}}\sum_{a=1}^{f}\chi(a)\left(\frac{N!}{(N+n)!}f^{n}-\frac{a^{n}}{n!}\right).$

Also

(6.13) $\displaystyle B_{N,n,\chi}(x)$
$\displaystyle=(-1)^{n}n!\left|\begin{array}[]{ccccc}\frac{N!}{(N+1)!}f&1&0&&\\\
\frac{N!}{(N+2)!}f^{2}&\frac{N!}{(N+1)!}f&&&\\\ \vdots&&&&0\\\
\frac{N!}{(N+n-1)!}f^{n-1}&\cdots&\frac{N!}{(N+2)!}f^{2}&\frac{N!}{(N+1)!}f&1\\\
\widehat{S}_{n}(x)&\widehat{S}_{n-1}(x)&\cdots&\widehat{S}_{2}(x)&\widehat{S}_{1}(x)\end{array}\right|\,$

where

$\widehat{S}_{n}(x)=\frac{1}{f^{N}}\sum_{a=1}^{f}\chi(a)\left(\frac{N!}{(N+n)!}f^{n}-\frac{(x+a)^{n}}{n!}\right)\,.$

Remark. When $\chi=\mathbf{1}$ is the trivial character of conductor $1$, by
$\widehat{S}_{n}=N!/(N+n)!$, the determinant (6.7) is reduced to

$B_{N,n,{\mathbf{1}}}=(-1)^{n}n!\left|\begin{array}[]{ccccc}\frac{N!}{(N+1)!}&1&0&&\\\
\frac{N!}{(N+2)!}&\frac{N!}{(N+1)!}&&&\\\ \vdots&&&&0\\\
\frac{N!}{(N+n-1)!}&\cdots&\frac{N!}{(N+2)!}&\frac{N!}{(N+1)!}&1\\\
\frac{N!}{(N+n)!}&\frac{N!}{(N+n-1)!}&\cdots&\frac{N!}{(N+2)!}&\frac{N!}{(N+1)!}\end{array}\right|\,$

which is [2, Theorem 1].

###### Proof of Theorem 3..

For simplicity, let

$\alpha_{i}=\frac{N!f^{i}}{(N+i)!}\quad(i\geq 1)\,,$

then

$\widehat{S}_{n}=\alpha_{n}\frac{S_{0}}{f^{N}}-\frac{S_{n}}{f^{N}n!}=\frac{1}{f_{N}}\left(\alpha_{n}S_{0}-\frac{S_{n}}{n!}\right)\,.$

Furthermore, by putting $\beta_{n}=(-1)^{n}B_{N,n,\chi}/n!$, we shall show
that

(6.19) $\displaystyle\beta_{n}=\left|\begin{array}[]{ccccc}\alpha_{1}&1&0&&\\\
\alpha_{2}&\alpha_{1}&&&\\\ \vdots&&&&0\\\
\alpha_{n-1}&\cdots&\alpha_{2}&\alpha_{1}&1\\\
\widehat{S}_{n}&\widehat{S}_{n-1}&\cdots&\widehat{S}_{2}&\widehat{S}_{1}\end{array}\right|$

by using the recurrence relation

(6.20)
$\beta_{n}=\frac{(-1)^{n}S_{n}}{f^{N}n!}-\sum_{i=0}^{n-1}(-1)^{n-i}\alpha_{n-i}\beta_{i}\quad(n\geq
1)\,.$

By definition,

$\beta_{0}=B_{N,0,\chi}=\frac{S_{0}}{f^{N}}\,.$

For $n=1$, by (6.20),

$\beta_{1}=\alpha_{1}\beta_{0}-\frac{S_{1}}{f^{N}}=\widehat{S}_{1}\,.$

Hence (6.19) is valid. Assume that (6.19) is valid up to $n-1$. Then by
expanding the determinant along the first row repeatedly, the right-hand side
of (6.19) is equal to

$\displaystyle\alpha_{1}\beta_{n-1}-\left|\begin{array}[]{ccccc}\alpha_{2}&1&0&&\\\
\alpha_{3}&\alpha_{1}&&&\\\ \vdots&\vdots&&&0\\\
\alpha_{n-1}&\alpha_{n-3}&\cdots&\alpha_{1}&1\\\
\widehat{S}_{n}&\widehat{S}_{n-2}&\cdots&\widehat{S}_{2}&\widehat{S}_{1}\end{array}\right|$
$\displaystyle=\alpha_{1}\beta_{n-1}-\alpha_{2}\beta_{n-2}+\cdots+(-1)^{n-1}\alpha_{n-2}\beta_{2}+(-1)^{n}\left|\begin{array}[]{cc}\alpha_{n-1}&1\\\
\widehat{S}_{n}&\widehat{S}_{1}\end{array}\right|$
$\displaystyle=\alpha_{1}\beta_{n-1}-\alpha_{2}\beta_{n-2}+\cdots+(-1)^{n-1}\alpha_{n-2}\beta_{2}+(-1)^{n}\alpha_{n-1}\beta_{1}+(-1)^{n+1}\widehat{S}_{n}$
$\displaystyle=\beta_{n}\,.$

The last equality is due to the recurrence relation (6.20). Then the first
relation (6.7) is proved. The second equation (6.13) can be proved similarly,
by using the recurrence relation

$\frac{B_{N,n,\chi}(x)}{n!}=\frac{S_{n}(x)}{f^{N}n!}-\sum_{i=0}^{n-1}\frac{N!f^{n-i}}{(N+n-i)!}\frac{B_{N,i,\chi}(x)}{i!}$

or

$\beta_{n}(x)=\frac{(-1)^{n}S_{n}(x)}{f^{N}n!}-\sum_{i=0}^{n-1}(-1)^{n-i}\alpha_{n-i}\beta_{i}(x)\quad(n\geq
1)\,.$

∎

## Acknowledgements

The authors are grateful to three anonymous referees for their precious
comments and advices. This work has been partly done when the second author
(T.K.) visited Harish-Chandra Research Institute in February 2020, where the
first author (K. C.) was working before moving to his current institute. T. K.
thanks K. C. for the warm hospitality of the institute. This work was
initiated from a hint by Professor Miho Aoki of Shimane University. Both
authors thank her.

## References

  * [1] T. Agoh, Shortened recurrence relations for generalized Bernoulli numbers and polynomials, J. Number Theory 176 (2017), 149–173.
  * [2] M. Aoki, T. Komatsu and G. K. Panda, Several properties of hypergeometric Bernoulli numbers, J. Inequal. Appl. 2019 (2019), Paper No.113. 24p.
  * [3] T. Arakawa, T. Ibukiyama and M. Kaneko, Bernoulli numbers and zeta functions. With an appendix by Don Zagier, Springer Monographs in Mathematics, Springer, Tokyo, 2014.
  * [4] K.-W. Chen and M. Eie, A note on generalized Bernoulli numbers, Pacific J. Math. 199 (2001), 41–59.
  * [5] J. W. L. Glaisher, Expressions for Laplace’s coefficients, Bernoullian and Eulerian numbers etc. as determinants, Messenger (2) 6 (1875), 49–63.
  * [6] A. Hassen and H. D. Nguyen, Hypergeometric Bernoulli polynomials and Appell sequences, Int. J. Number Theory 4 (2008), 767–774.
  * [7] A. Hassen and H. D. Nguyen, Hypergeometric zeta functions, Int. J. Number Theory 6 (2010), 99–126.
  * [8] F. T. Howard, A sequence of numbers related to the exponential function, Duke Math. J. 34 (1967), 599–615.
  * [9] F. T. Howard, Some sequences of rational numbers related to the exponential function, Duke Math. J. 34 (1967), 701–716.
  * [10] S. Hu and M.-S. Kim, On hypergeometric Bernoulli numbers and polynomials, Acta Math. Hungar. 154(1) (2018), 134–146.
  * [11] S. Hu and M.-S. Kim, Two closed forms for the Apostol-Bernoulli polynomials, Ramanujan J. 46 (2018), 103–117.
  * [12] S. Hu and T. Komatsu, Explicit expressions for the related numbers of higher order Appell polynomials, Quaest. Math. 43(8) (2020), 1019–1029.
  * [13] K. Kamano, Sums of products of hypergeometric Bernoulli numbers, J. Number Theory 130 (2010), 2259–2271.
  * [14] T. Komatsu, Hypergeometric Cauchy numbers, Int. J. Number Theory 9 (2013), 545–560.
  * [15] T. Komatsu and P. Yuan, Hypergeometric Cauchy numbers and polynomials, Acta Math. Hungar. 153 (2017), 382–400.
  * [16] T. Komatsu and W. Zhang, Several properties of multiple hypergeometric Euler numbers, Tokyo J. Math. 42(2) (2019), 551–570.
  * [17] M. Merca, Bernoulli numbers and symmetric functions, Rev. R. Acad. Cienc. Exactas Fis. Nat. Ser. A Mat. RACSAM 114(1) (2020), Paper No.20, 16 pp.
  * [18] H. D. Nguyen, Sums of products of hypergeometric Bernoulli polynomials, MAA-NJ Section - Spring Meeting, Middle Country College, NJ, April 10, 2010.
  * [19] L. C. Washington, Introduction to cyclotomic fields, second edition, GTM 83, Springer-Verlag, New York, 1997.
  * [20] L. Zhu, New bounds for the ratio of two adjacent even-indexed Bernoulli numbers, Rev. R. Acad. Cienc. Exactas Fis. Nat. Ser. A Mat. RACSAM 114(2) (2020), Paper No. 83, 13 pp.

## Appendix

We derive first few values of $B_{N,n,\chi}$:

$\displaystyle B_{N,0,\chi}$ $\displaystyle=\frac{S_{0}}{f^{N}}\,,$
$\displaystyle B_{N,1,\chi}$
$\displaystyle=-\frac{S_{0}}{(N+1)f^{N-1}}+\frac{S_{1}}{f^{N}}\,,$
$\displaystyle B_{N,2,\chi}$
$\displaystyle=\frac{2S_{0}}{(N+1)^{2}(N+2)f^{N-2}}-\frac{2S_{1}}{(N+1)f^{N-1}}+\frac{S_{2}}{f^{N}}\,,$
$\displaystyle B_{N,3,\chi}$
$\displaystyle=\frac{3!(N-1)S_{0}}{(N+1)^{3}(N+2)(N+3)f^{N-3}}+\frac{3!S_{1}}{(N+1)^{2}(N+2)f^{N-2}}-\frac{3S_{2}}{(N+1)f^{N-1}}+\frac{S_{3}}{f^{N}}\,,$
$\displaystyle B_{N,4,\chi}$
$\displaystyle=\frac{4!(N^{3}-N^{2}-6N+2)S_{0}}{(N+1)^{4}(N+2)^{2}(N+3)(N+4)f^{N-4}}+\frac{4!(N-1)S_{1}}{(N+1)^{3}(N+2)(N+3)f^{N-3}}$
$\displaystyle\qquad+\frac{12S_{2}}{(N+1)^{2}(N+2)f^{N-2}}-\frac{4S_{3}}{(N+1)f^{N-1}}+\frac{S_{4}}{f^{N}}\,.$

Now

$\displaystyle B_{N,0}$ $\displaystyle=1\,,$ $\displaystyle B_{N,1}$
$\displaystyle=-\frac{1}{N+1}\,,$ $\displaystyle B_{N,2}$
$\displaystyle=\frac{2}{(N+1)^{2}(N+2)}\,,$ $\displaystyle B_{N,3}$
$\displaystyle=\frac{3!(N-1)}{(N+1)^{3}(N+2)(N+3)}\,,$ $\displaystyle B_{N,4}$
$\displaystyle=\frac{4!(N^{3}-N^{2}-6N+2)}{(N+1)^{4}(N+2)^{2}(N+3)(N+4)}\,.$

From ([2]) we can see that,

$B_{N,n,\chi}=\sum_{k=0}^{n}\binom{n}{k}\frac{S_{n-k}B_{N,k}}{f^{N-k}}\,,$

as shown in Corollary 1.

When $N=1$ and $\chi$ is not the trivial character, as $S_{0}=0$, we find that

$\displaystyle B_{1,0,\chi}$ $\displaystyle=0\,,$ $\displaystyle B_{1,1,\chi}$
$\displaystyle=\frac{S_{1}}{f}\,,$ $\displaystyle B_{1,2,\chi}$
$\displaystyle=-S_{1}+\frac{S_{2}}{f}\,,$ $\displaystyle B_{1,3,\chi}$
$\displaystyle=\frac{fS_{1}}{2}-\frac{3S_{2}}{2}+\frac{S_{3}}{f}\,,$
$\displaystyle B_{1,4,\chi}$ $\displaystyle=fS_{2}-2S_{3}+\frac{S_{4}}{f}\,.$