# Determination of QPO properties in the presence of strong broad-band noise:
a case study on the data of MAXI J1820+070

Deng-Ke Zhou1,2, Shuang-Nan Zhang1,2, Li-Ming Song1,2, Jin-Lu Qu1,2, Liang
Zhang1,3, Xiang Ma1, You-Li Tuo1, Ming-Yu Ge1, Yanan Wang3, Shu Zhang1 and
Lian Tao1
1Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing, China
2University of Chinese Academy of Sciences, Chinese Academy of Sciences,
Beijing, China.
3Physics & Astronomy, University of Southampton, Southampton, Hampshire SO17
1BJ, UK Corresponding author. E-mail<EMAIL_ADDRESS>author.
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Accurate calculation of the phase lags of quasi-periodic oscillations (QPOs)
will provide insight into their origin. In this paper we investigate the phase
lag correction method which has been applied to calculate the intrinsic phase
lags of the QPOs in MAXI J1820+070. We find that the traditional additive
model between BBN and QPOs in the time domain is rejected, but the convolution
model is accepted. By introducing a convolution mechanism in the time domain,
the Fourier cross-spectrum analysis shows that the phase lags between QPOs
components in different energy bands will have a simple linear relationship
with the phase lags between the total signals, so that the intrinsic phase
lags of the QPOs can be obtained by linear correction. The power density
spectrum (PDS) thus requires a multiplicative model to interpret the data. We
briefly discuss a physical scenario for interpreting the convolution. In this
scenario, the corona acts as a low-pass filter, the Green’s function
containing the noise is convolved with the QPOs to form the low-frequency part
of the PDS, while the high-frequency part requires an additive component. We
use a multiplicative PDS model to fit the data observed by _Insight_ -HXMT.
The overall fitting results are similar compared to the traditional additive
PDS model. Neither the width nor the centroid frequency of the QPOs obtained
from each of the two PDS models were significantly different, except for the
r.m.s. of the QPOs. Our work thus provides a new perspective on the coupling
of noise and QPOs.

###### keywords:

X-rays: binaries – methods: analytical – methods: data analysis

††pubyear: 2021††pagerange: Determination of QPO properties in the presence of
strong broad-band noise: a case study on the data of MAXI
J1820+070–Determination of QPO properties in the presence of strong broad-band
noise: a case study on the data of MAXI J1820+070

##  1 Introduction

Decades of research on black hole binaries (BHBs) show that their X-ray
emission is variable on different time scales, including the low-frequency
(mHz to 30Hz) quasi-periodic oscillations (QPOs) and the broad-band noise
(BBN) (Psaltis et al. 1999; Ingram et al. 2009; Ingram & Done 2011; Motta
2016). The study of the timing signals can effectively diagnose the geometric
characteristics of the disk and the corona near the black hole (Belloni &
Hasinger 1990; Belloni et al. 2002; Ingram 2016; Ingram & Motta 2019). The
disk and the corona near the black hole continuously radiate X-ray photons
outward due to various radiation mechanisms (thermal radiation, Compton
Radiation and so on). Photons with different energy arrive at the observer at
different times because they may come from different radiation regions (Lin et
al. 2000; Rapisarda et al. 2016), or undergo different scattering processes
(Cui 1999; Poutanen 2001), or have very complex mechanisms that cause delays
(Morgan et al. 1997; Wijnands et al. 1999; Qu et al. 2010). Therefore,
analyzing the phase/time lags of photons between different energy bands helps
us better understand the geometric or radiometric characteristics of X-ray
BHBs. A common analysis method is based on Fourier cross-spectrum, which
measures the frequency-dependent phase lag spectrum (FDPLS) between the
signals in two different energy bands (van der Klis et al. 1987). This method
allows to study the phase lags of two signals as a function of Fourier
frequencies. Thus the phase lags between different components of timing
signals, which usually originate from different physical processes (Narayan &
Yi 1995; Done et al. 2007; Ingram & Done 2011), can be studied separately. For
example, Zhang et al. (2020) conducted a systematic study on the phase lag of
the type-C QPO and found that the phase lag behaviour of the sub-harmonic of
the QPO is very similar to that of the QPO fundamental component but the
second harmonic of the QPO shows a quit different phase lag behaviour. Uttley
et al. (2011) investigated the phase lag of the BBN components of GX-339 and
found that the large lags can be explained by viscous propagation of mass
accretion fluctuations in the disk.

The traditional way to obtain the phase lag of the QPO components is to assume
that the other components contribute weakly to the lag in the QPO frequency
range, and then directly treat the values in the QPO frequency range as the
phase lag of the QPO components (e.g., Morgan et al. 1997; Wijnands et al.
1999; Kara et al. 2019; Zhang et al. 2020). However, the coexistence of
various components makes it difficult to calculate any of the individual
component. In particular, when the BBN is sufficiently strong in the QPO
frequency range, there is no reason to ignore the effect of the BBN on the
measured QPO phase lag. Despite attempts by some authors to ameliorate this
dilemma by fitting different components of the cross-spectrum (e.g., Qu et al.
2004), there is no broad consensus on how to obtain the intrinsic phase lag of
the QPO in the presence of strong BBN. Therefore, it is difficult to determine
the intrinsic properties (including phase lag) of the QPO in the presence of
strong BBN.

In a recent work (Ma et al. 2021, hereafter Ma21) the authors attempt to
correct for the original phase lags, which gives a clear physical picture
using the corrected phase lags. Ma21 investigated the behavior of the QPO
phase lags in MAXI J1820+070 using _Insight_ -HXMT observations and proposed a
method to obtain the intrinsic phase lag of the QPO. In their analysis of the
phase lags, they find that by subtracting the phase lags below the QPO
frequency range they can obtain consistent QPO phase lags as functions of
photon energy for all observations and can explain the lag behavior through
the precession of a compact jet above the black hole. On the data they used,
the PDS shows that the BBN components are too strong compared to the QPO
component to ignore the contribution to the phase lags in the QPO frequency
range (see panel c of figure 1 in Ma21). If they do not correct the phase lags
for the QPO, the phase lags obtained from the original FDPLS will be affected
by the BBN and thus are not intrinsic phase lags of the QPO. Although Ma21
applied this method to obtain consistent results of the phase lags, the
rationale for doing so was not explained in detail, so the plausibility of
this correction method needs to be tested. On the data they analyzed, some of
the observations obtained phase lags with little difference before and after
the correction, but some of the phase lags changed significantly (even the
sign is totally reversed) before and after the correction. Therefore, we
believe it is necessary to investigate under what conditions the correction is
effective and how the QPO component is related to the BBN component. The
motivation of this paper is to explore the mechanism behind the correction
method used by Ma21 and to investigate the QPO properties in the presence of
strong BBN in conjunction with the results obtained by Ma21.

Since we want to obtain the properties of a certain component (in our case,
the QPO), and what we observe is some kind of superposition of all components,
we have to face the problem of how these components contribute to the total
signal. In this paper, when we refer to the term signal, we are referring to
the light curve or the underlying time series. The total signal is defined as
the time series that we directly observe and the sub-signals are the sub-
components such as QPO and BBN that make up the total signal. Traditionally,
it is believed that the BBN and the QPO are additive in the time domain and
that they are incoherent at any frequency, which is why the PDS is fitted by
the sum of several Lorentzian functions. Ingram & van der Klis (2013) proposed
a possible relationship between QPO and BBN, where the QPO component and the
BBN component are multiplied in the time domain. In this case, a convolution
model is required for the fitting of the PDS in the frequency domain. Another
way in which the QPO component and the BBN component are combined into a total
signal in the time domain is convolution, which is usually caused by the
response of the QPO signal in the region where the BBN component is generated
(a model similar to this mechanism can be found in Cabanac et al. 2010). The
calculation of the FDPLS involves the Fourier orthogonal decomposition of the
signal, so it can be expected that if the sub-signals form the total signal in
different ways, then the relationship between the FDPLS of the total signal
and the FDPLS of the sub-signals must be different.

This paper is structured as follows: Section 2 analyzes the relationships
between the FDPLS as well as the PDS of total signals and sub-signals. An
algorithm to generate two signals satisfied specific PDS and FDPLS
simultaneously is also proposed. Besides, one possible way of coupling the QPO
component and the BBN component in the time domain is discussed. In Section 3,
based on the results of Ma21’s analysis of MAXI J1820+070 on phase lags, we
argue that the QPO component and the BBN component constitute the total signal
by convolution in the time domain. Using the data of MAXI J1820+070, we fit
the PDS in different energy bands using the multiplicative PDS model and the
traditional additive PDS model, and compare their differences. In addition, we
also performed some simulations to rule out the possibility that the total
signal appears to be the sum of the sub-signals in the time domain. Section 4
discusses and summarizes the whole paper.

##  2 theory and simulation

###  2.1 phase lag relationship

Suppose that the expressions of non-zero mean signals $r_{1}(t)$, $r_{2}(t)$,
$q_{1}(t)$, $q_{2}(t)$ at frequency $f_{0}$ can be written as:

$\displaystyle r_{1}(t)$ $\displaystyle=R_{1}\sin(2\pi
f_{0}t+\phi_{r_{1}})+c_{r_{1}},$ (1) $\displaystyle r_{2}(t)$
$\displaystyle=R_{2}\sin(2\pi f_{0}t+\phi_{r_{2}})+c_{r_{2}},$ $\displaystyle
q_{1}(t)$ $\displaystyle=Q_{1}\sin(2\pi f_{0}t+\phi_{q_{1}})+c_{q_{1}},$
$\displaystyle q_{2}(t)$ $\displaystyle=Q_{2}\sin(2\pi
f_{0}t+\phi_{q_{2}})+c_{q_{2}},$

where $c_{r_{1}}$, $c_{r_{2}}$, $c_{q_{1}}$ and $c_{q_{2}}$ are the mean
values of the corresponding signals; $R_{1}$, $R_{2}$, $Q_{1}$, $Q_{2}$ and
$\phi_{r_{1}}$, $\phi_{r_{2}}$, $\phi_{q_{1}}$, $\phi_{q_{2}}$ are the
amplitudes and the initial phases of the corresponding signals, respectively.
The frequency $f_{0}$ can take any non-negative value including $0$. When $0$
is taken, it indicates that this is a constant signal. If the total signal is
the sum of the sub-signals in the time domain (hereafter this kind of total
signal is called the additive signal), i.e. :

$\displaystyle s_{1}(t)$ $\displaystyle=r_{1}(t)+q_{1}(t),$ (2) $\displaystyle
s_{2}(t)$ $\displaystyle=r_{2}(t)+q_{2}(t),$

then the phase difference (i.e. phase lag) between $s_{1}(t)$ and $s_{2}(t)$
can be written as:

$\displaystyle\Delta\phi_{\rm add}(s_{2},s_{1};f_{0})$ (3)
$\displaystyle=\phi_{s_{2}}-\phi_{s_{1}}$ $\displaystyle={\rm
Arg}[R_{2}\cos(\phi_{r_{2}})+Q_{2}\cos(\phi_{q_{2}}),R_{2}\sin(\phi_{r_{2}})+Q_{2}\sin(\phi_{q_{2}})]$
$\displaystyle-{\rm
Arg}[R_{1}\cos(\phi_{r_{1}})+Q_{1}\cos(\phi_{q_{1}}),R_{1}\sin(\phi_{r_{1}})+Q_{1}\sin(\phi_{q_{1}})].$

Here we use $\rm Arg$ $[a,b]$ to denote the argument of the complex $a+ib$,
where $i$ is the imaginary unit. It can be seen from equation (3) that if the
total signal is the additive signal, the phase lag between the total signals
depends on the amplitude and initial phase of each sub-signal.

If the total signal is convoluted by the sub-signals (hereafter this kind of
total signal is called the convolved signal), then the phase lag between the
total signal and the phase lag between the sub-signals satisfies a linear
relationship, the proof of which will be given below. Still assume that the
sub-signals satisfy equation (1), but at this time the total signals are equal
to the convolution of the sub-signals:

$\displaystyle s_{1}(t)$ $\displaystyle=r_{1}(t)\otimes q_{1}(t),$ (4)
$\displaystyle s_{2}(t)$ $\displaystyle=r_{2}(t)\otimes q_{2}(t),$

where the sign $\otimes$ represents the convolution operation. The Fourier
transform of the convolution of two signals is equal to the multiplication of
their respective Fourier transforms. We can obtain the cross-correlation
function (CCF) of $s_{1}(t)$ and $s_{2}(t)$ in the frequency domain:

$\displaystyle{\rm CCF}(f)$
$\displaystyle=\frac{R_{1}R_{2}Q_{1}Q_{2}}{64\pi^{2}}e^{-i[\Delta\phi(r_{2},r_{1})+\Delta\phi(q_{2},q_{1})]}\delta^{4}(f-f_{0}),$
(5)

where $\Delta\phi(r_{2},r_{1})=\phi_{r_{2}}-\phi_{r_{1}}$ and
$\Delta\phi(q_{2},q_{1})=\phi_{q_{2}}-\phi_{q_{1}}$ are the phase lag of the
sub-signals. The phase lag of the two total signals $s_{1}$ and $s_{2}$ can be
obtained by taking the argument of their CCF:

$\displaystyle\Delta\phi_{\rm con}(s_{2},s_{1};f_{0})$ $\displaystyle={\rm
Arg}[{\rm CCF}(f)]=\Delta\phi(r_{2},r_{1})+\Delta\phi(q_{2},q_{1}).$ (6)

That is to say, if the total signal is the convolved signal, the phase lag of
the total signals is equal to the sum of the phase lag of the sub-signals.

We also note that Rapisarda et al. (2014) argued that the QPO component is
multiplied together with the broad component to form the observed signal. We
now consider the phase lag relationship between the total signal composed of
single frequency sub-signals by multiplying them together (hereafter, this
kind of total signal is called the multiplicative signal). $s_{1}(t)$ and
$s_{2}(t)$ now are written as

$\displaystyle s_{1}(t)$ $\displaystyle=r_{1}(t)\times q_{1}(t)$ (7)
$\displaystyle=c_{q_{1}}R_{1}\sin(2\pi
f_{0}t+\phi_{r_{1}})+c_{r_{1}}Q_{1}\sin(2\phi f_{0}t+\phi_{q_{1}})$
$\displaystyle+\frac{1}{2}R_{1}Q_{1}\cos(2\pi\times
2f_{0}t+\phi_{r_{1}}+\phi_{q_{1}})+c_{s_{1}},$ $\displaystyle s_{2}(t)$
$\displaystyle=r_{2}(t)\times q_{2}(t)$ $\displaystyle=c_{q_{2}}R_{2}\sin(2\pi
f_{0}t+\phi_{r_{2}})+c_{r_{2}}Q_{2}\sin(2\phi f_{0}t+\phi_{q_{2}})$
$\displaystyle+\frac{1}{2}R_{2}Q_{2}\cos(2\pi\times
2f_{0}t+\phi_{r_{2}}+\phi_{q_{2}})+c_{s_{2}},$

where $c_{s_{1}}$ and $c_{s_{2}}$ are constants. Thus, the two total signals
$s_{1}(t)$ and $s_{2}(t)$ contain two non-zero frequency components, one at
$f_{0}$ and the other at $2f_{0}$. We can see that the first two terms of
$s_{1}(t)$ and $s_{2}(t)$ are in fact additive signals and thus the results on
additive signals can be used. Thus, the phase lags of them can be written as:

$\displaystyle\Delta\phi_{\rm mul}(s_{2},s_{1};f_{0})$
$\displaystyle=\Delta\phi_{\rm add}(s^{\prime}_{2},s^{\prime}_{1};f_{0}),$ (8)
$\displaystyle\Delta\phi_{\rm mul}({s_{2}},{s_{1}};2f_{0})$
$\displaystyle=\Delta\phi(r_{2},r_{1})+\Delta\phi(q_{2},q_{1}),$

where $s^{\prime}_{1}=c_{q_{1}}R_{1}\sin(2\pi
f_{0}t+\phi_{r_{1}})+c_{r_{1}}Q_{1}\sin(2\pi f_{0}t+\phi_{q_{1}})$ and
$s^{\prime}_{2}=c_{q_{2}}R_{2}\sin(2\pi
f_{0}t+\phi_{r_{2}})+c_{r_{2}}Q_{2}\sin(2\pi f_{0}t+\phi_{q_{2}})$. This is
very interesting because the multiplicative signal seems to contain properties
of both additive and convolved signals: on one hand the phase lag at frequency
$f_{0}$ follows the pattern of the additive signal and on the other hand the
phase lag at frequency $2f_{0}$ follows the pattern of the convolved signal.
However, in general the mean value of the actual signal is larger than its
amplitude, so it is expected that the total FDPLS of the multiplicative signal
should be closer to the pattern of the additive signal, as we will see in the
simulation section.

For the general signals $r_{n}$, $q_{n}$, $s_{n}$ ($n=0,1,....N-1$), their
discrete-time Fourier series are

$\begin{split}r_{n}=\frac{1}{N}\sum_{k=0}^{N-1}R_{k}e^{i2\pi\frac{k}{N}n},\\\
q_{n}=\frac{1}{N}\sum_{k=0}^{N-1}Q_{k}e^{i2\pi\frac{k}{N}n},\\\ \
s_{n}=\frac{1}{N}\sum_{k=0}^{N-1}S_{k}e^{i2\pi\frac{k}{N}n},\end{split}$ (9)

where $R_{k}$, $Q_{k}$, and $S_{k}$ are the discrete Fourier transforms of
$r_{n}$, $q_{n}$, and $s_{n}$, respectively. Thus $r_{n}$, $q_{n}$, $s_{n}$
can be treated as a superposition of many trigonometric functions with
different amplitude, different frequencies and different initial phases. For
the additive and convolved signals discuss above, at the specified frequency,
these signals have the same properties as the corresponding single frequency
signals. So the phase lags of the additive signal follows equation (3) at each
frequency, and the phase lags of the convolved signal follows equation (6) at
each frequency. In summary, for the additive/convolved signals, their FDPLS
follow equation (3) or (6) at each specific frequency, respectively. In this
case, each frequency corresponds to a set of parameters (amplitude and initial
phase) for calculating the phase lag. Unfortunately, it is clear from equation
(7) that additional frequency components appear in the multiplicative signal
that are not identical to the sub-signals, so the conclusion for the single-
frequency signal cannot be generalized to the general signal. Nevertheless, we
can use simulation (see sec 2.4) to explore the phase lag relationship between
the multiplicative signals.

###  2.2 PDS relationship

Assuming $s(t)$ is the additive signal, i.e. $s(t)=r(t)+q(t)$. Let $P_{\rm
s}(f)$, $P_{\rm r}(f)$ and $P_{\rm q}(f)$ be the PDS of the signals $s(t)$,
$r(t)$ and $q(t)$, respectively. Considering the Fourier transform is linear,
one obtains:

$\displaystyle P_{s}(f)$ $\displaystyle=|\mathscr{F}(r+q)|^{2}$ (10)
$\displaystyle=|\mathscr{F}(r)+\mathscr{F}(q)|^{2}$
$\displaystyle=P_{r}(f)+P_{q}(f)+\mathscr{F}(r)^{*}\mathscr{F}(q)+\mathscr{F}(r)\mathscr{F}(q)^{*}.$

The last two terms are actually the cross-spectrum of the signals $r(t)$ and
$q(t)$. If $r(t)$ and $q(t)$ are incoherent at all frequencies, then their
cross-spectrum will converge to zero after averaging many signal realizations.
Therefore the above equation is simplified to

$<P_{\rm s}(f)>=<P_{r}(f)>+<P_{q}(f)>,$ (11)

where the <> sign indicates the average of many realizations of the signals.
This indicates that the PDS of the sum of two incoherent signals is equal to
the sum of their respective PDS.

Assuming $s(t)$ is the convolved signal, i.e. $s(t)=r(t)\otimes q(t)$. The
Fourier transform of the total signal $s(t)$ is equal to the multiplication of
the Fourier transforms of the sub-signals $r(t)$, $q(t)$, i.e.

$\displaystyle\mathscr{F}(s)=\mathscr{F}(r)\mathscr{F}(q),$ (12)
$\displaystyle<P_{s}(f)>=<P_{r}(f)><P_{q}(f)>.$

That is, the PDS of the convolved signal is equal to the multiplication of the
PDS of the corresponding sub-signals. It can be easily generalized that if
$s(t)=r(t)\otimes q(t)+p(t)$, and $r(t)\otimes q(t)$ is incoherent with
$p(t)$, their PDS will satisfy

$<P_{\rm s}(f)>=<P_{\rm r}(f)><P_{\rm q}(f)>+<P_{\rm p}(f)>.$ (13)

The PDS properties for convolved signals can be generalized to multiplicative
signals simply based on the symmetry of the Fourier transform, that is, the
PDS of the multiplicative signal is the convolution of the PDS of the sub-
signals.

###  2.3 An algorithm for simultaneously simulating signals with specified PDS
and FDPLS

In order to verify the correctness of the above theoretical analysis as well
as to facilitate the analysis below, some simulations need to be done. An
algorithm is thus needed to generate two signals with specified PDS and
specified FDPLS simultaneously. The algorithm steps are as follows:

  1. 1)

Use Timmer & Koenig (1995) (TK95 in the following) algorithm to generate two
signals $s(t)$ and $s^{\prime}(t)$ that satisfy the specified PDS. Because the
phase given to the signal by TK95 algorithm is random, the phase lag between
these two signals is now on average zero. Denote their Fourier transforms as
$S(f)$, $S^{\prime}(f)$, respectively.

  2. 2)

Given the FDPLS $\phi(f)$, then calculate ${\rm CCF}(f)$ according to the
following equation:

${\rm
CCF}(f)=\begin{cases}|S(f)||S^{\prime}(f)|\\{\cos[\phi(f)]+i\sin[\phi(f)]\\}&\text{$f\neq
0$},\\\ |S(f)||S^{\prime}(f)|&\text{$f=0$},\end{cases}$ (14)

where $i$ is the imaginary unit.

  3. 3)

The complex array ${\rm CCF}(f)$ obtained in step 2 is divided by the complex
conjugate of $S(f)$ to obtain a new complex array. Then performing inverse
Fourier transform to it to obtain the signal $s^{\prime\prime}(t)$. Expressed
in mathematical notation, it is

$s^{\prime\prime}(t)=\mathscr{F}^{-1}[\frac{{\rm CCF}(f)}{S^{*}(f)}].$ (15)

The underlying PDS of $s^{\prime\prime}(t)$ is the same as the PDS of $s(t)$,
but the FDPLS between $s^{\prime\prime}(t)$ and $s(t)$ will satisfy the given
FDPLS. In summary, $s(t)$ and $s^{\prime\prime}(t)$ satisfy both the given PDS
and the given FDPLS. In this paper, all PDS and FDPLS are extracted using the
X-ray astronomy python package stingray (Huppenkothen et al. 2019, version
0.3), and all PDS and FDPLS fitting are done by XSPEC (Arnaud 1996, version
12.11.1) or lmfit (Newville et al. 2014, version 1.0.2).

###  2.4 simulation

Table 1: Timing properties of $r_{1}(t)$, $r_{2}(t)$, $q_{1}(t)$, $q_{2}(t)$ (see section 2.4 for the definition of them). signals | bin size (s) | $\nu_{c}$ | $\omega$ | mean rate (cts/s) | exposure (s) | fractional r.m.s | PDS type | FDPLS type
---|---|---|---|---|---|---|---|---
$r_{1}(t)$ | 0.01 | 0 | 3 | 2000 | 2000 | 30% | BBN | constant
$r_{2}(t)$ | 0.01 | 0 | 4 | 2000 | 2000 | 20% | BBN
$q_{1}(t)$ | 0.01 | 1 | 0.1 | 2000 | 2000 | 15% | QPO | dip
$q_{2}(t)$ | 0.01 | 1 | 0.2 | 2000 | 2000 | 10% | QPO
Figure 1: Simulation results of the FDPLS and PDS. Panel a: PDS of signals
$r_{1}(t)$, $r_{2}(t)$, $q_{1}(t)$, $q_{2}(t)$. Panel b: simulated FDPLS (blue
dots) between $r_{1}(t)$ and $r_{2}(t)$ and simulated FDPLS (black dots)
between $q_{1}(t)$ and $q_{2}(t)$. Panel c: simulated FDPLS (green dots)
between $r_{2}(t)+q_{2}(t)$ and $r_{1}(t)+q_{1}(t)$ and simulated FDPLS (blue
dots) between $r_{2}(t)\times q_{2}(t)$ and $r_{1}(t)\times r_{2}(t)$. Panel
d: simulated FDPLS between $r_{2}(t)\otimes q_{2}(t)$ and $r_{1}(t)\otimes
q_{1}(t)$. The red curves in panels b, c, and d are theoretically calculated
curves. Error bars correspond to $1\sigma$ confidence intervals.

Four signals $r_{1}(t)$, $r_{2}(t)$, $q_{1}(t)$, and $q_{2}(t)$ with time
resolution of 0.01 s are simulated according to the algorithm proposed in
subsection 2.3. The PDS of all these signals are characterized by the
Lorentzian function, which takes the form of

$L(f)=\frac{K(\omega/(2\pi))}{(\omega/2)^{2}+(f-f_{c})^{2}},$ (16)

where $K$, $\omega$ and $f_{c}$ denote the normalization factor, the full
width at half maximum (FWHM) and the centroid frequency, respectively. The PDS
of $r_{1}(t)$ and $r_{2}(t)$ are modeled by setting the centroid frequency of
the Lorentzian function to zero and taking a large $\omega$, which simulates
BBN, while $q_{1}(t)$ and $q_{2}(t)$ are modeled by taking the appropriate
non-zero centroid frequency and $\omega$, which simulates QPO. In addition,
the theoretical FDPLS $\phi(f)$ is also set. The FDPLS between BBNs is set to
be constant, while the FDPLS between QPOs is set to have a dip-like feature
near the centroid frequency (as seen in MAXI J1820+070). That is,

$\phi(f)=\begin{cases}0.5&\text{for BBN},\\\
-0.5e^{\frac{(f-1)^{2}}{0.04}}&\text{for QPO}.\end{cases}$ (17)

The timing properties of these four signals are summarized in table 1. We then
split each signal into multiple 20-sec segments and calculated the PDS of each
segment with Leahy normalization (Leahy et al. 1983). The PDS is rebined by a
logarithmic factor of 0.03 and we finally obtain the averaged PDS with the
frequency range of 0.05-50.53 Hz. The FDPLS is obtained using cross-spectrum
analysis.

The results of the simulated PDS and the FDPLS are shown in panel a and panel
b of Fig 1, respectively. When the total signal is assumed to be additive or
multiplicative signal, the FDPLS between the total signals is shown in panel c
of Fig. 1. When the total signal is assumed to be the convolved signal, the
FDPLS between the total signals is shown in panel d of Fig. 1. As stated in
the theoretical analysis section, the FDPLS of the multiplicative signal is
very close to the FDPLS of the additive signal (see the green data points and
the blue data points in panel c of Fig. 1). Due to the symmetry of the Fourier
transform to convolution and multiplication, the PDS section will only compare
the differences between convolved and additive signals. In panels b, c, d of
Fig. 1, the data points are obtained by simulation and the red dashed lines
are obtained from our theoretical calculation (the theoretical curve drawn in
panel c of Fig. 1 is for the additive signal, and we did not draw the
theoretical curve for the multiplicative signal because of the analytical
difficulties). The theoretical curve shown in panel c of Fig. 1 is calculated
by using the value of the simulated data (i.e., amplitude, initial phase), and
it appears to fluctuate around the data points, which is due to the randomness
deliberately introduced by the simulation algorithm (see TK95 for detail). The
difference between panel c and panel d of Fig. 1 is mainly due to the
different dependence of the FDPLS on the different kinds of signals (additive,
multiplicative, convolved signals) on each sub-signal. The FDPLS between the
convolved signals depends only on the FDPLS between sub-signals, independent
of the other properties of the sub-signals. This is not the case for the
additive and multiplicative signals. So it can be seen from panel c that the
FDPLS depends on the relative power of sub-signals, while in panel d the FDPLS
does not depend on the shape of the PDS of the sub-signals. In conclusion, the
simulation results of the FDPLS are consistent with the theoretical analysis.
The results of the PDS simulation results are shown in Fig. 2.

The PDS of $r_{1}(t)$ and $q_{1}(t)$ are shown in the left panel of Fig. 2,
and the PDS of the additive and convolved signals are shown in the middle and
right panels of Fig. 2, respectively. We can see that the PDS of the additive
signal is the sum of the PDS of the sub-signals, while the PDS of the
convolved signal is the multiplication of the PDS of the sub-signals. The
solid lines running through the data points in the PDS are the best-fit using
the additive and multiplicative Lorentzian models for the additive signal and
the convolved signal, respectively. Overall, the simulation results are in
good agreement with those of the theoretical analysis.

Figure 2: PDS of simulation results. Left panel: the cyan and blue data points
are the PDS of the signals $r_{1}(t)$ and $q_{1}(t)$, respectively. The solid
lines are the best-fit using Lorentzian model (considering the contribution of
Poisson noise requires adding a constant to the Lorentzian model). Middle
panel: the blue points are the PDS of the sum of the signals $r_{1}(t)$ and
$q_{1}(t)$. The red solid line is the best-fit using two summed Lorentzian
functions. Right panel: the blue points are the PDS of the convolution of the
signals $r_{1}(t)$ and $q_{1}(t)$. The red solid line is the best-fit using
two multiplicative Lorentzian functions (considering the contribution of
Poisson noise requires adding a constant to each Lorentzian model). Error bars
correspond to $1\sigma$ confidence intervals.

###  2.5 A possible mechanism for introducing a convolution mechanism in the
time domain

Assuming that the orbit of matter around a black hole is circular and
Keplerian. The equation can be derived based on the conservation of mass and
angular momentum (e.g. Ingram 2016), i.e.

$\frac{\partial\Sigma}{\partial t}=\frac{3}{R}\frac{\partial}{\partial
R}[R^{\frac{1}{2}}\frac{\partial}{\partial R}(\nu\Sigma R)],$ (18)

where $\Sigma=\rho H$ is the surface density of the corona or disk, and $\nu$
is the kinematic viscosity. Assuming that the surface density at $t=0$ is
$\Sigma(t=0,R)=\delta(R-R_{0})$ and $\nu$ is a constant, we will get

$g(R,t)=\frac{m}{12\pi\nu
t}(\frac{R}{R_{0}})^{-\frac{1}{4}}I_{\frac{1}{4}}(\frac{RR_{0}}{6\nu
t})e^{-\frac{R_{0}^{2}+R^{2}}{12\nu t}},$ (19)

where $I_{\frac{1}{4}}$ is the modified Bessel function and $g(R,t)$ is called
the Green’s function of the system. Under the condition that the system is
linear, the surface density of any initial fluctuation $q(t)$ at position
$R=R_{0}$ is the convolution of that fluctuation with the Green’s function,
i.e., $\Sigma(R,t)=q(t)\otimes g(R,t)$ (Ingram 2016). Denoting the mass
accretion rate as $\dot{M}(R,t)$, then the luminosity corresponding to such a
accretion rate is $L(R,t)\propto\dot{M}(R,t)\propto\Sigma(R,t)\propto
q(t)\otimes g(R,t)$. If we check the region $R<<R_{0}$, then we will get
$g(R,t)\propto t^{\frac{5}{4}}e^{-\frac{R_{0}^{2}}{12\nu t}}$. The PDS of such
a damped exponential signal is a zero-centred Lorentzian function (Ingram
2016). By introducing two types of white noise, one associated with the
Green’s function and the other superimposed on the QPO signal, we assume that
the observed signal is expressed in the time domain as $s(t)=g(t)\otimes
wn_{1}\otimes[q(t)+wn_{2}]$. We have assumed that q(t) has the form of QPO.
Considering that the white noise and QPO signals are incoherent, a PDS of the
combined signal will has the form

$P(R,f)\propto P_{b1}(R,f)P_{q}(R,f)+P_{b2}(R,f),$ (20)

where $P_{b1}$ denotes the first zero-centred Lorentzian function (i.e., the
BBN1 component), $P_{q}$ denotes the non-zero centred Lorentzian function
(i.e., the QPO component) and $P_{b2}$ denotes the second zero-centred
Lorentzian function (i.e., the BBN2 component). Note that the former term of
the above summation is due to the fluctuation propagation in the form of QPO
and the latter term is due to the fluctuation propagation in the form of white
noise, which dominates different frequency ranges (we will see this in section
3).

Furthermore, it is worth noting that the above result is valid only when
$R<<R_{0}$ and the assumptions about the white noise and QPO fluctuations are
satisfied. The total observed luminosity is the integral of the differential
luminosity over the entire corona after considering the emissivity (Ingram &
Done 2011), but the form is very complicated. Nonetheless, it is still
worthwhile to start with a simple model to explain the data. For this reason,
when fitting the PDS of the real data with the multiplicative PDS model in
section 3, only a form similar to equation (20) will be considered.

Figure 3: QPO phase lag correction of three typical _Insight_ -HXMT
observations (reproduced from the data in Ma21). Each row represents the
results of one observation. The FDPLS between 1-2.6 keV and 100-150 keV energy
bands are shown in the left panels. Middle panels are the original QPO energy-
dependent phase lags. Right panels are the intrinsic QPO energy-dependent
phase lags after correction. The intrinsic QPO phase lags are obtained by
subtracting the average of the phase lag of the BBN component (marked by the
red dots on the left panels) from the original QPO phase lag (the averaged
value marked by the cyan dots on left panels). The red arrows indicate the
high energy band data that we will model in Fig. 6. Error bars correspond to
$1\sigma$ confidence intervals. Figure 4: Several examples of the best-fit of
the observed PDS data (ObsID P0114661078) in different energy bands with
additive or multiplicative PDS models. The PDS model used in the left panels
is the additive PDS model (i.e., equation (21)). The PDS model used in the
right panels is the multiplicative PDS model (i.e., equation (22)). The
contribution of Poisson noise in all PDS has been subtracted. Error bars
correspond to $1\sigma$ confidence intervals. Figure 5: The parameters of the
QPO (fundamental and harmonic components) as functions of photon energy which
are obtained with the traditional additive and multiplicative PDS model,
respectively. Error bars correspond to $1\sigma$ confidence intervals.

##  3 The FDPLS and PDS of MAXI J1820+070

MAXI J1820+070 was discovered by the Monitor of All-sky X-ray Image (MAXI)
during the outburst on 11 March 2018 (Kawamuro et al. 2018). It was confirmed
to be a BHB (Torres et al. 2019). _Insight_ -HXMT carried out observations
three days after its discovery and obtained rich data with total exposure time
of over 2000 ks. Ma21 has carried out a detailed temporal analysis of these
data and, in particular, detailed calculations of the phase lags in different
energy bands. Fig. 3 shows three typical observations from top to bottom, with
a clear dip-like feature appearing near the QPO frequency range (the averaged
value of this frequency range shown by the cyan dots denotes the original QPO
phase lag). The phase lags of low-frequency BBN component are marked with red
dots, denoted as background phase lag. The intrinsic phase lag of the QPO is
obtained by subtracting the average of the background phase lag from the
original QPO phase lag. After the correction, the absolute value of phase lag
of the QPO becomes larger as the energy increases in all three observations.
For the sake of clarity, the detailed correction steps used in Ma21 are re-
summarized as follows:

  1. 1)

Calculate the FDPLS and identify the ccentroid frequency $f_{0}$ and the FWHM
$\omega$ of the QPO according to the PDS.

  2. 2)

The original phase lag of the QPO is defined as the average of the phase lags
in the frequency range $f_{0}\pm\omega/2$.

  3. 3)

The background phase lag is defined as the average of the phase lags below the
QPO frequency range.

  4. 4)

The intrinsic phase lag of the QPO is then defined as the original phase lag
minus the background phase lag.

Such a correction actually implies two assumptions. The first assumption is
that the phase lags of the BBN component at the QPO frequency range are the
same as the phase lags below the QPO frequency range, at least their averaged
values must be approximately equal. The second assumption is that the total
phase lags (i.e., the observed original phase lags) at the QPO frequency range
are equal to the sum of the BBN component phase lags and the QPO intrinsic
phase lags. The correction of the phase lags is valid only when these two
assumptions are satisfied simultaneously. The first assumption can be
considered to be approximately satisfied. This is because the FDPLS obtained
from MAXI J1820+070 shows that the phase lag does not vary significantly with
frequency below the QPO frequency range. Thus it is reasonable to assume that
at the QPO frequency range, the phase lags of the BBN component are
approximately equal to the phase lags below the QPO frequency range. As to
whether the second assumption can be satisfied, we need to first make an
assumption about how the BBN component and the QPO component synthesize the
observed signal. From the discussion of section 2 we know that if the observed
signal is considered to be the sum of the BBN component and the QPO component
(which is the default assumption in most of the literatures), the second
condition cannot be satisfied. The second condition can be satisfied only when
the observed signal is considered as a convolution of the BBN component and
the QPO component.

In the case that the observed signal is the convolved signal, the PDS of the
signal needs to be fitted by a multiplicative PDS model. In section 2.5 we
introduced a multiplicative PDS model, so we use the multiplicative PDS model
to fit the PDS in different energy bands and make a comparison with the
traditional additive PDS model. We are not going to explore the multiplicative
signal (i.e., the total signal is the multiplication of the sub-signals in the
time domain) further because we are mainly concerned with the additive and
convolved signals. For the additive PDS model, the form is

$\begin{split}P_{\rm
add}(f)=&\frac{K_{1}(\omega_{1}/(2\pi))}{(\omega_{1}/2)^{2}+(f-f_{c_{1}})^{2}}+\frac{K_{2}(\omega_{2}/(2\pi))}{(\omega_{2}/2)^{2}+(f-f_{c_{2}})^{2}}\\\
&+\frac{K_{3}(\omega_{3}/(2\pi))}{(\omega_{3}/2)^{2}+(f-f_{c_{3}})^{2}}+\frac{K_{4}(\omega_{4}/(2\pi))}{(\omega_{4}/2)^{2}+(f-f_{c_{4}})^{2}}+c.\end{split}$
(21)

For the multiplicative PDS model, the form is

$\begin{split}P_{\rm
mul}(f)=&\frac{K_{1}(\omega_{1}/(2\pi))}{(\omega_{1}/2)^{2}+(f-f_{c_{1}})^{2}}\times[\frac{K_{2}(\omega_{2}/(2\pi))}{(\omega_{2}/2)^{2}+(f-f_{c_{2}})^{2}}\\\
&+\frac{K_{3}(\omega_{3}/(2\pi))}{(\omega_{3}/2)^{2}+(f-f_{c_{3}})^{2}}]+\frac{K_{4}(\omega_{4}/(2\pi))^{2}}{(\omega_{4}/2)^{2}+(f-f_{c_{4}})^{2}}+c,\end{split}$
(22)

which is consistent with the multiplicative PDS model we discussed in section
2.5. We use these two models to fit the PDS in different energy bands in a
representative _Insight_ -HXMT observation (ObsID P0114661078). The data
reduction process in this paper is the same as Ma21. We extract the light
curves with the time resolution of 0.03125 s in each energy band (1-2.6 keV,
2.6-4.8 keV, 4.8-7.0 keV, 7-11 keV, 11-23 keV, 25-35 keV, 35-48 keV, 48-67
keV, 67-100 keV, 100-150 keV and 150-200 keV). We then split the light curves
into multiple 32-sec segments and calculate the PDS of each segment with
Miyamoto normalization (Miyamoto et al. 1991) for the convenience of
calculating fractional r.m.s later, and finally obtain the averaged PDS with
the frequency range of 1/32 to 16 Hz. After subtracting the contribution of
Poisson noise in the PDS, we fitted the PDS with the additive and
multiplicative PDS models. Some fitting examples in different energy bands are
shown in Fig. 4. From top to bottom, Fig. 4 shows the PDS fitting results for
the three energy bands. The left panels are fitted using the traditional
additive PDS model while the right panels are fitted using the multiplicative
PDS model. In the multiplicative PDS model, the low-frequency zero-centered
Lorentzian component is multiplied onto the QPO component instead of being
added, resulting in the left side of the QPO component being lifted up in the
right panels of Fig. 4. We can see that the total fitting results are similar
but the individual components have some differences. Basing on the best-fit,
we calculate the centroid frequency, the FWHM and the fractional r.m.s111The
fractional r.m.s is calculated in the same way as for Bu et al. (2015), but
ignoring the background correction, since the correction coefficients are the
same for both models and our aim is only to compare their differences.
Neglecting the background correction leads to a lower fractional r.m.s for the
energy band with a lower signal-to-noise ratio (in this paper it is the higher
energy band), but it does not change our conclusion. of the QPO and parameters
of the BBN on each energy band. And the results are listed in table 2, 3 and
4. As shown in Fig. 5, the centroid frequency and the FWHM of the QPO as
functions of photon energy calculated according to the two models are similar,
but the fractional r.m.s given by the two models are significantly different.
For the fundamental component of QPO, the fractional r.m.s of QPO given by the
traditional additive PDS model is about $2\sim 3$ times higher than that of
the multiplicative PDS mode, but the trend is the same for both results. For
the harmonic component of the QPO, the difference in the fractional r.m.s
given by the two PDS models is not very significant.

Table 2: Best-fit results for the fundamental frequency component of the QPO
obtained using additive and multiplicative PDS models (i.e. equations (21) and
(22)), respectively.

[b]  QPO frequency (Hz) QPO FWHM (Hz) energy band (keV) additive PDS model
multiplicative PDS model additive PDS model multiplicative PDS model 1.0-2.6
$0.44\pm 0.01$ $0.45\pm 0.01$ $0.13\pm 0.04$ $0.14\pm 0.03$ 2.6-4.8 $0.43\pm
0.01$ $0.45\pm 0.01$ $0.14\pm 0.03$ $0.14\pm 0.03$ 4.8-7.0 $0.43\pm 0.01$
$0.45\pm 0.01$ $0.13\pm 0.05$ $0.12\pm 0.04$ 7.0-11.0 $0.43\pm 0.01$ $0.45\pm
0.01$ $0.11\pm 0.06$ $0.11\pm 0.05$ 11.0-23.0 $0.43\pm 0.01$ $0.44\pm 0.01$
$0.12\pm 0.03$ $0.10\pm 0.03$ 25.0-35.0 $0.42\pm 0.01$ $0.43\pm 0.01$ $0.11\pm
0.03$ $0.09\pm 0.02$ 35.0-48.0 $0.43\pm 0.01$ $0.44\pm 0.01$ $0.10\pm 0.03$
$0.08\pm 0.02$ 48.0-67.0 $0.43\pm 0.01$ $0.44\pm 0.01$ $0.10\pm 0.03$ $0.08\pm
0.02$ 67.0-100.0 $0.43\pm 0.01$ $0.44\pm 0.01$ $0.11\pm 0.03$ $0.10\pm 0.02$
100.0-150.0 $0.42\pm 0.01$ $0.43\pm 0.01$ $0.06\pm 0.03$ $0.09\pm 0.03$
150.0-200.0 $0.44\pm 0.03$ $0.45\pm 0.03$ $0.10\pm 0.13$ $0.08\pm 0.11$ QPO
r.m.s %* reduced $\chi^{2}$ energy band (keV) additive PDS model
multiplicative PDS model additive PDS model multiplicative PDS model 1.0-2.6
$9.14\pm 0.31$ $3.83\pm 0.10$ 0.48 0.71 2.6-4.8 $10.26\pm 0.31$ $4.15\pm 0.07$
0.74 0.88 4.8-7.0 $8.99\pm 0.40$ $3.36\pm 0.09$ 0.61 0.69 7.0-11.0 $7.59\pm
0.56$ $3.06\pm 0.13$ 0.66 0.66 11.0-23.0 $7.86\pm 0.20$ $2.59\pm 0.05$ 0.48
0.57 25.0-35.0 $5.65\pm 0.11$ $1.72\pm 0.02$ 1.02 1.55 35.0-48.0 $5.54\pm
0.13$ $1.72\pm 0.03$ 0.89 1.14 48.0-67.0 $4.75\pm 0.09$ $1.44\pm 0.02$ 0.59
0.71 67.0-100.0 $5.24\pm 0.13$ $1.64\pm 0.02$ 0.45 0.58 100.0-150.0 $3.88\pm
0.14$ $1.88\pm 0.04$ 0.43 0.44 150.0-200.0 $2.05\pm 0.12$ $0.76\pm 0.06$ 0.44
0.41

  * *

No background correction is applied to the r.m.s because we are only
interested in the difference between the results of fitting using the additive
PDS model and the multiplicative PDS model, and the correction factors are the
same for both models.

Table 3: The same as table 2 but for the harmonic frequency component of the QPO. | QPO frequency(Hz) | QPO FWHM(Hz)
---|---|---
energy band (keV) | additive PDS model | multiplicative PDS model | additive PDS model | multiplicative PDS model
1.0-2.6 | $0.90\pm 0.02$ | $0.98\pm 0.02$ | $0.37\pm 0.08$ | $0.54\pm 0.09$
2.6-4.8 | $0.90\pm 0.02$ | $0.94\pm 0.02$ | $0.34\pm 0.09$ | $0.33\pm 0.08$
4.8-7.0 | $0.91\pm 0.02$ | $0.94\pm 0.02$ | $0.26\pm 0.09$ | $0.24\pm 0.10$
7.0-11.0 | $0.91\pm 0.06$ | $0.98\pm 0.06$ | $0.47\pm 0.21$ | $0.47\pm 0.18$
11.0-23.0 | $0.88\pm 0.02$ | $0.89\pm 0.02$ | $0.10\pm 0.10$ | $0.12\pm 0.12$
25.0-35.0 | $0.91\pm 0.03$ | $0.96\pm 0.03$ | $0.42\pm 0.12$ | $0.35\pm 0.11$
35.0-48.0 | $0.85\pm 0.09$ | $0.85\pm 0.23$ | $0.42\pm 0.34$ | $0.41\pm 0.34$
48.0-67.0 | $0.93\pm 0.04$ | $0.97\pm 0.02$ | $0.75\pm 0.09$ | $0.36\pm 0.08$
67.0-100.0 | $0.91\pm 0.02$ | $0.95\pm 0.02$ | $0.31\pm 0.07$ | $0.25\pm 0.06$
100.0-150.0 | $0.94\pm 0.03$ | $1.01\pm 0.03$ | $0.45\pm 0.14$ | $0.39\pm 0.13$
150.0-200.0 | $1.01\pm 0.04$ | $1.12\pm 0.04$ | $0.39\pm 0.14$ | $0.65\pm 0.12$
| QPO r.m.s % |  |
energy band (keV) | additive PDS model | multiplicative PDS model |  |
1.0-2.6 | $9.55\pm 1.22$ | $10.04\pm 0.58$ |  |
2.6-4.8 | $9.47\pm 1.28$ | $7.76\pm 0.79$ |  |
4.8-7.0 | $8.57\pm 1.65$ | $6.63\pm 1.30$ |  |
7.0-11.0 | $9.29\pm 1.89$ | $7.57\pm 0.98$ |  |
11.0-23.0 | $4.53\pm 1.19$ | $3.35\pm 2.41$ |  |
25.0-35.0 | $8.76\pm 1.44$ | $5.93\pm 0.92$ |  |
35.0-48.0 | $5.55\pm 2.10$ | $6.16\pm 5.68$ |  |
48.0-67.0 | $9.42\pm 0.69$ | $4.51\pm 0.55$ |  |
67.0-100.0 | $6.37\pm 0.88$ | $4.04\pm 0.77$ |  |
100.0-150.0 | $5.59\pm 1.07$ | $3.54\pm 0.71$ |  |
150.0-200.0 | $4.43\pm 0.71$ | $4.77\pm 0.35$ |  |
Table 4: The same as table 2 but for the BBN components. | BBN1 FWHM (Hz) | BBN2 FWHM (Hz)
---|---|---
energy band (keV) | additive PDS model | multiplicative PDS model | additive PDS model | multiplicative PDS model
1.0-2.6 | $0.48\pm 0.08$ | $0.27\pm 0.04$ | $3.17\pm 0.51$ | $3.01\pm 0.35$
2.6-4.8 | $0.62\pm 0.19$ | $0.28\pm 0.04$ | $3.89\pm 0.76$ | $3.02\pm 0.44$
4.8-7.0 | $0.87\pm 0.47$ | $0.30\pm 0.06$ | $3.70\pm 1.61$ | $2.54\pm 0.61$
7.0-11.0 | $0.54\pm 0.24$ | $0.32\pm 0.09$ | $10.27\pm 9.14$ | $5.56\pm 3.56$
11.0-23.0 | $0.92\pm 1.92$ | $0.42\pm 0.38$ | $2.56\pm 1.59$ | $2.25\pm 0.44$
25.0-35.0 | $0.69\pm 0.35$ | $0.38\pm 0.07$ | $4.70\pm 1.89$ | $3.26\pm 0.95$
35.0-48.0 | $0.58\pm 0.80$ | $\cdots$ | $6.24\pm 15.27$ | $2.06\pm 2.09$
48.0-67.0 | $0.65\pm 0.20$ | $0.41\pm 0.05$ | $10.83\pm 2.89$ | $2.77\pm 0.52$
67.0-100.0 | $0.69\pm 0.54$ | $0.43\pm 0.09$ | $2.74\pm 0.82$ | $2.44\pm 0.35$
100.0-150.0 | $0.71\pm 0.45$ | $0.50\pm 0.08$ | $3.76\pm 1.71$ | $3.10\pm 0.87$
150.0-200.0 | $1.86\pm 0.18$ | $0.61\pm 0.08$ | $18.88\pm 12.81$ | $6.55\pm 2.43$
| BBN1 r.m.s % | BBN2 r.m.s %
energy band (keV) | additive PDS model | multiplicative PDS model | additive PDS model | multiplicative PDS model
1.0-2.6 | $16.32\pm 1.30$ | $\cdots$ | $20.12\pm 1.18$ | $22.07\pm 0.67$
2.6-4.8 | $14.36\pm 2.21$ | $\cdots$ | $20.56\pm 1.19$ | $23.23\pm 0.73$
4.8-7.0 | $14.40\pm 5.28$ | $\cdots$ | $18.60\pm 3.27$ | $22.31\pm 1.29$
7.0-11.0 | $11.55\pm 2.11$ | $\cdots$ | $15.69\pm 5.61$ | $15.50\pm 2.27$
11.0-23.0 | $21.02\pm 0.75$ | $\cdots$ | $19.22\pm 7.50$ | $20.85\pm 1.03$
25.0-35.0 | $9.87\pm 2.65$ | $\cdots$ | $14.66\pm 1.34$ | $16.88\pm 0.92$
35.0-48.0 | $11.34\pm 1.83$ | $\cdots$ | $7.87\pm 3.83$ | $11.38\pm 4.97$
48.0-67.0 | $8.01\pm 1.10$ | $\cdots$ | $10.53\pm 0.93$ | $13.04\pm 0.59$
67.0-100.0 | $6.19\pm 4.02$ | $\cdots$ | $13.18\pm 1.96$ | $14.56\pm 0.56$
100.0-150.0 | $5.33\pm 2.21$ | $\cdots$ | $9.17\pm 1.30$ | $10.51\pm 0.66$
150.0-200.0 | $10.49\pm 0.56$ | $\cdots$ | $11.46\pm 4.92$ | $9.34\pm 0.77$
Figure 6: Modeling results for the FDPLS between signals $s_{\rm l}$ and $s_{\rm h}$ and their respective PDS (see section 3 for the definitions of $s_{\rm l}$ and $s_{\rm h}$). In panels a and b, the solid lines are the best-fit using the models proposed in section 3 and the best-fit parameters are listed in table 2. The yellow bands running through the two panels are the QPO FWHM frequency ranges (including fundamental and harmonic frequencies). Error bars correspond to $1\sigma$ confidence intervals. Table 5: Timing properties of $s_{\rm lq}$,$s_{\rm lnq}$,$s_{\rm hq}$ and $s_{\rm hnq}$ (see section 3 for the definition of them).  signals | bin size (s) | mean rate (cts/s) | exposure (s) | fractional r.m.s | PDS type | FDPLS type
---|---|---|---|---|---|---
$s_{\rm lq}$ | 0.03125 | 142.5 | 8000 | 27.00% | QPO | dip
$s_{\rm hq}$ | 0.03125 | 63.5 | 8000 | 9.43% | QPO
$s_{\rm lnq}$ | 0.03125 | 142.5 | 8000 | 51.00% | Non-QPO | constant
$s_{\rm hnq}$ | 0.03125 | 63.5 | 8000 | 15.70% | Non-QPO
Table 6: PDS and FDPLS fitting results for $s_{\rm l}$ and $s_{\rm h}$. $s_{\rm l}$ PDS model |  $s_{\rm h}$ PDS model | FDPLS model
---|---|---
parameter name | value | parameter name | value | parameter name | value
$K_{1}$ | $0.007\pm 0.002$ | $K_{1}$ | $0.001\pm 0.000$ | $A_{1}$ | $-0.109\pm 0.015$
$f_{c1}$ | $0.442\pm 0.008$ | $f_{c1}$ | $0.430\pm 0.010$ | $\mu_{1}$ | $0.405\pm 0.011$
$\omega_{1}$ | $0.119\pm 0.033$ | $\omega_{1}$ | $0.090\pm 0.030$ | $\sigma_{1}$ | $0.062\pm 0.008$
$K_{2}$ | $0.009\pm 0.002$ | $K_{2}$ | $0.001\pm 0.000$ | $A_{2}$ | $-0.204\pm 0.041$
$f_{c2}$ | $0.900\pm 0.020$ | $f_{c2}$ |  $0.860$ (frozen) | $\mu_{2}$ | $0.776\pm 0.037$
$\omega_{2}$ | $0.370\pm 0.083$ | $\omega_{2}$ |  $0.200$ (frozen) | $\sigma_{2}$ | $0.185\pm 0.024$
$K_{3}$ | $0.027\pm 0.004$ | $K_{3}$ |  $0.006$ (frozen) | c | $1.102\pm 0.040$
$f_{c3}$ |  $0.000$ (frozen) | $f_{c3}$ |  $0.000$ (frozen) | $\cdots$ | $\cdots$
$\omega_{3}$ | $0.477\pm 0.080$ | $\omega_{3}$ | $2.061\pm 0.425$ | $\cdots$ | $\cdots$
$K_{4}$ | $0.040\pm 0.005$ | c | $0.016\pm 0.000$ | $\cdots$ | $\cdots$
$f_{c4}$ |  $0.000$ (frozen) | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
$\omega_{4}$ | $3.175\pm 0.512$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
c | $0.007\pm 0.000$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$

In addition, we would like to know how the phase lags should look like in MAXI
J1820 for additive and convolved signals, so some simulations are done. We use
the same data (ObsID P0114661078) used above as an example to show how
effective this correction is. From panel h and panel i of Fig. 3 we can see
that the QPO phase lags between the high and low energy bands before and after
correction are completely flipped. We focus on two energy bands of these
signals: the reference energy band (1-2.6 keV) and the high energy band
(100-150 keV). The phase lag between the reference energy band and the high
energy band is indicated by the red arrows in panel h and panel i of Fig. 3.
After the data reduction in the same way as Ma21, we extract two light curves
in reference energy band (denote as $s_{\rm l}$) and high energy band (denote
as $s_{\rm h}$). The mean count rate of $s_{\rm l}$ and $s_{\rm h}$ are 285
counts/s and 127 counts/s, respectively, and both of the effective exposure
time are 8 ks. The FDPLS and PDS of $s_{\rm l}$ and $s_{\rm h}$ are first
fitted to obtain the best models, and after that, the best models are used for
simulations. For the FDPLS between $s_{\rm l}$ and $s_{\rm h}$, the model
takes the form

${\rm
lag}(f)=\frac{A_{1}}{\sigma_{1}\sqrt{2\pi}}e^{[{-{(f-\mu_{1})^{2}}/{{2\sigma_{1}}^{2}}}]}+\frac{A_{2}}{\sigma_{2}\sqrt{2\pi}}e^{[{-{(f-\mu_{2})^{2}}/{{2\sigma_{2}}^{2}}}]}+c,$
(23)

where the first and second terms represent the dip-like phase lags of the QPO
components (including fundamental and harmonic frequencies) and the last term
represent the phase lags of the non-QPO components (BBN components). For the
PDS of $s_{\rm l}$ and $s_{\rm h}$, a constant term and the sum of four
Lorentzian functions are used to fit the data, i.e., the PDS model is the same
as equation (21). We find that the above PDS model does not require the high-
frequency zero-centred Lorentzian component for $s_{\rm h}$ when fitting the
PDS of $s_{\rm h}$. Moreover, in fitting the PDS of $s_{\rm h}$ we found that
the harmonic frequency component of QPO is not well constrained due to the low
signal-to-noise ratio of the data. We thus fix the parameters of the QPO
harmonic frequency component, which does not affect the goodness of fit, but
is useful for our subsequent simulation of the QPO components. The fitting
results are shown in Fig. 6 and the best-fit parameters of the above models
are listed in table 6. It is worth pointing out that we are using an additive
PDS model to fit the PDS here, which is correct for additive signal, but not
for convolved signal, which should use a multiplicative PDS model. However, we
note that the FDPLS relationship of the convolved signal depends only on the
FDPLS of the sub-signals and is independent of the PDS of the sub-signals, so
the PDS model we use here has no effect on the FDPLS calculation of the
convolved signal.

We then simulate four sub-signals based on the best FDPLS and PDS models
obtained above. The PDS of the QPO component is modeled using the sum of the
non-zero centred Lorentzian components and the PDS of the BBN component is
modeled using the sum of the zero-centred Lorentzian components. We first
simulate four signals using TK95 algorithm, noted as $s_{\rm lq}$, $s_{\rm
lnq}$, $s_{\rm hq}$, $s_{\rm hnq}$, which stand for the QPO component and the
BBN component in the reference energy band, and the QPO component and the BBN
component in the high energy band, respectively. After that, we use the
algorithm proposed in section 2.3 to make the FDPLS of $s_{\rm lq}$ and
$s_{\rm hq}$ satisfy the Gaussian components of the best-fit model and make
the FDPLS of $s_{\rm lnq}$ and $s_{\rm hnq}$ satisfy the constant component of
the best-fit model. The timing properties of these four signal are listed in
table 5. Note that in calculating the FDPLS we split the signal into QPO and
non-QPO components, which in effect assumes that the contribution of the
additive component of equation (20) to the overall FDPLS can be neglected,
i.e. the destruction of this additive component to the additivity of the phase
lag between the convolution components can be neglected, as we will explain in
detail in the subsequent discussion section. The PDS of $s_{\rm lq}$, $s_{\rm
lnq}$, $s_{\rm hq}$, $s_{\rm hnq}$ are shown in top panel of Fig. 7. The FDPLS
between $s_{\rm lq}$ and $s_{\rm hq}$ are two dips, and the FDPLS between
$s_{\rm lnq}$ and $s_{\rm hnq}$ are constant, which are shown in the middle
panel of Fig. 7.

Figure 7: Simulation results based on the data modeling shown in Fig. 6. Upper
panel: the simulated PDS of the signals $s_{\rm lq}$,$s_{\rm lnq}$,$s_{\rm
hq}$ and $s_{\rm hnq}$. Middle panel: the simulated FDPLS (blue dots) of
$s_{\rm lq}$ and $s_{\rm hq}$ and the simulated FDPLS (green dots) of $s_{\rm
lnq}$ and $s_{\rm hnq}$. Bottom panel: the calculated FDPLS of the total
signal when the total signal is an/a additive/convolved signal (cyan/red data
points). The gray dots are the observed FDPLS. The yellow bands running
through panels are the QPO FWHM frequency ranges (including fundamental and
harmonic frequencies). Error bars correspond to $1\sigma$ confidence
intervals.

Then, $s_{\rm lq}$ and $s_{\rm lnq}$ are added/convolved to get the
additive/convolved signal. $s_{\rm hq}$ and $s_{\rm hnq}$ are added/convolved
to get the other additive/convolved signal. The FDPLS of the
additive/convolved signals is shown by the cyan/red dotted lines in the lower
panel of Fig. 7. The gray dots in the figure are the observed data. As can be
seen from Fig. 7, the simulated results for the additive signal are very
different from those given by the data. But the simulated results for the
convolved signal match the data perfectly. This simulation result indicates
that the observed data can distinguish between convolved and additive signals
for FDPLS.

##  4 discussion and summary

In this paper, we investigate the mechanism behind the phase lag correction
that was successfully applied for the first time by Ma21 for MAXI J1820+070,
where the strong BBN and the QPO coexist. After correcting the phase lag of
the QPO, the absolute value of QPO phase lag increases monotonically with
photon energy in all observations. Ma21 explained the phase lag behavior of
the QPO by employing a compact jet with precession. In this scenario, the
high-energy photons come from the part of the jet closer to the black hole,
and the precession of the compact jet causes the QPO phenomenon and allows the
high-energy photons to reach the observer first, resulting in a soft lag.
Because the phase lag behavior can have a large impact on physical
conclusions, it is necessary to investigate the rationality of this correction
method. Since we want to obtain the intrinsic properties of the QPO, and what
we observe is some kind of superposition of the QPO and the BBN components, we
have to face the question of how these components constitute the total signal.
We found that the correction method is effective only when the sub-signals are
synthesized into the total signal by convolution. If the total signal is the
convolved signal, the intrinsic phase lags of the QPO can be obtained by
subtracting the phase lags of the BBN component from the original phase lags
of the total signal, as successfully implemented in Ma21.

If the observed total signals are convolved signals, the corresponding PDS
cannot be fitted simply by summing a series of Lorentzian functions
(conventions in most of the literatures) but require a multiplicative PDS
model. We then try to introduce the convolution mechanism by assuming the
propagation of the QPO waves in the corona (may be due to the magneto-acoustic
wave propagating within the corona, e.g. Cabanac et al. 2010). The fluctuation
propagation in the form of Dirac delta function resulting in the Green’s
function. Any form of timing fluctuation will be the convolution of that
fluctuation and the Green’s function (Ingram 2016). We assume that the Green’s
function is first convolved with the white noise and then convolved with the
QPO signal to form the low-frequency part of the observed signal, while the
high-frequency part is the result of the convolution of the Green’s function
with the two white noise components. If the Green’s function and the QPO
signal are convoluted in the time domain, the total PDS will be the
multiplication of their respective PDS according to the convolution theorem.
Based on this, we introduce a multiplicative PDS model to fit the observed PDS
in a representative _Insight_ -HXMT observation in 11 different energy bands.
For comparison, we also fitted the same data using the traditional additive
PDS model. Overall, both additive and multiplicative PDS models fit the
observed data well, but the individual components have some differences. The
two models give little difference in the centroid frequency as well as in the
FWHM of the QPO . For the fundamental frequency component of the QPO, the
fraction of r.m.s of the QPO given by the traditional additive PDS model is
about $2\sim 3$ higher than that of the multiplicative PDS model, but the
trend is the same for both results. For the harmonic frequency component of
QPO, the fractional rms given by the two models are not significantly
different.

For the traditional additive PDS model, the low-frequency zero-centred
Lorentzian component can be considered as the variability due to the
propagation of the white noise fluctuation in the outer region of the corona
(Ingram & Done 2011), and the narrow Lorentzian components stand for the
fundamental and harmonic components of the QPO, and the high frequency zero-
centred Lorentzian component is responsible for the variability due to the
propagation of the white noise fluctuation in the inner region of the corona
(Ingram & Done 2011). All these terms are simply added together, which means
that there is no coherence between them.

For the multiplicative PDS model, we find that the low frequency component of
the PDS can be fitted by multiplying the fundamental and harmonic components
of the QPO with a zero-centred Lorentzian function, in addition to an
additional additive component to produce the high frequency part of the PDS.
The additive Lorentzian component plays the same as the role in the additive
PDS model. Therefore, for our multiplicative model, the high-frequency part
does not need to be convolved to the QPO signal, but is simply added together.
This additive component appearing in the PDS model looks to destroy the
additivity of the phase lag brought about by the time domain convolution. We
ignored the contribution of this additive component in our previous phase lag
calculations. We make a simulation to investigate the effect of this
additional additive component on the total FDPLS. Due to the dependence of the
FDPLS of additive signals on the PDS of the individual components, we need to
know the PDS parameters of each component. Specifically, we first obtain the
PDS parameters for each component based on the results of the fit of the
multiplicative PDS model (shown in the right panels of Fig. 4), and then
simulate four signals based on the mean of these parameters in two energy
bands: $x_{1}$ (the convolved signal of energy band 1), $y_{1}$ (the
additional additive signal of energy band 1), $x_{2}$ (the convolved signal of
energy band 2) and $y_{2}$ (the additional additive signal of energy band 2).
The PDS of these four signals are shown in the upper panel of Fig. 8. We then
set the FDPLS model of $x_{1}$ and $x_{2}$ to
$\phi(f)=-0.5e^{-\frac{(f-0.4)^{2}}{0.1^{2}}}-0.3e^{-\frac{(f-0.8)^{2}}{0.15^{2}}}+1$
and set the FDPLS model of $y_{1}$ and $y_{2}$ to 1. Finally we calculate the
FDPLS of $x_{1}$ and $x_{2}$ and the FDPLS of $x_{1}+y1$ and $x_{2}+y_{2}$,
respectively. By comparing these two FDPLS, the effect of the additional
additive component on the total FDPLS can be known. As can be seen from the
lower panel of Fig. 8, the additional additive component has almost no effect
on the total FDPLS, except for a slight dilution of the phase lag of the
harmonic component. It is therefore reasonable to ignore the effect of the
additional additive component on the FDPLS of the convolution components in
our previous analysis.

Figure 8: The effect of an additional additive components on the total FDPLS.
Upper panel: PDS of simulated signals. Lower panel: FDPLS with or without an
additional additive component.

Traditionally, it is mostly assumed that the observed components are additive
in the time domain, and it has also been suggested that it might be more
reasonable to multiply these components in the time domain based on the
fluctuation propagation model (e.g. Ingram & van der Klis (2013)). However,
neither of these two models can explain the phase lag correction in Ma21. In
order to explain the correction of phase lags in Ma21, we propose a
convolution model instead of the additive and multiplicative models in the
time domain, which is supported by comparison between simulations and data on
both PDS and FDPLS. This suggests that the convolution model can explain the
behaviour of the phase lag observed in MAXI J1820+070, in which case the phase
lag correction method applied in Ma21 is correct.

Finally, it is worth pointing out that our current convolution model still has
limitations. For examples, it is not yet possible to explain the energy
dependence of the phase lag using the convolution model, and the relationship
of individual components to specific physical processes needs further
development. However, it is certain that at least part of the time domain
signal is filtered by the system before it reaches the observer (e.g. both the
accretion disk and corona/jet play the role of low-pass filters to some
extent), and these response processes are necessarily accompanied by time
domain convolution operations.

## Acknowledgments

This research has made use of the data from the _Insight_ -HXMT mission, a
project funded by China National Space Administration and the Chinese Academy
of Sciences. This work is supported by the National Natural Science Foundation
of China under grants 12133007, U1838201, U1938201.

## Data Availability

The data used in this paper can be found in the _Insight_ -HXMT website
(http://hxmtweb.ihep.ac.cn/).

## References

  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Belloni & Hasinger (1990) Belloni T., Hasinger G., 1990, A&A, 227, L33
  * Belloni et al. (2002) Belloni T., Psaltis D., van der Klis M., 2002, ApJ, 572, 392
  * Bu et al. (2015) Bu Q.-c., Chen L., Li Z.-s., Qu J.-l., Belloni T. M., Zhang L., 2015, ApJ, 799, 2
  * Cabanac et al. (2010) Cabanac C., Henri G., Petrucci P. O., Malzac J., Ferreira J., Belloni T. M., 2010, MNRAS, 404, 738
  * Cui (1999) Cui W., 1999, ApJ, 524, L59
  * Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&ARv, 15, 1
  * Huppenkothen et al. (2019) Huppenkothen D., et al., 2019, ApJ, 881, 39
  * Ingram (2016) Ingram A. R., 2016, Astronomische Nachrichten, 337, 385
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Ingram & Motta (2019) Ingram A. R., Motta S. E., 2019, New Astron. Rev., 85, 101524
  * Ingram & van der Klis (2013) Ingram A., van der Klis M., 2013, MNRAS, 434, 1476
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Kara et al. (2019) Kara E., et al., 2019, Nature, 565, 198
  * Kawamuro et al. (2018) Kawamuro T., et al., 2018, The Astronomer’s Telegram, 11399, 1
  * Leahy et al. (1983) Leahy D. A., Darbro W., Elsner R. F., Weisskopf M. C., Sutherland P. G., Kahn S., Grindlay J. E., 1983, ApJ, 266, 160
  * Lin et al. (2000) Lin D., Smith I. A., Böttcher M., Liang E. P., 2000, ApJ, 531, 963
  * Ma et al. (2021) Ma X., et al., 2021, Nature Astronomy, 5, 94
  * Miyamoto et al. (1991) Miyamoto S., Kimura K., Kitamoto S., Dotani T., Ebisawa K., 1991, ApJ, 383, 784
  * Morgan et al. (1997) Morgan E. H., Remillard R. A., Greiner J., 1997, ApJ, 482, 993
  * Motta (2016) Motta S. E., 2016, Astronomische Nachrichten, 337, 398
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Newville et al. (2014) Newville M., Stensitzki T., Allen D. B., Ingargiola A., 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Zenodo, doi:10.5281/zenodo.11813
  * Poutanen (2001) Poutanen J., 2001, in White N. E., Malaguti G., Palumbo G. G. C., eds, American Institute of Physics Conference Series Vol. 599, X-ray Astronomy: Stellar Endpoints, AGN, and the Diffuse X-ray Background. pp 310–325 (arXiv:astro-ph/0002505), doi:10.1063/1.1434644
  * Psaltis et al. (1999) Psaltis D., Belloni T., van der Klis M., 1999, ApJ, 520, 262
  * Qu et al. (2004) Qu J. L., Chen Y., Wu M., Chen L., Song L. M., 2004, Ap&SS, 293, 441
  * Qu et al. (2010) Qu J. L., Lu F. J., Lu Y., Song L. M., Zhang S., Ding G. Q., Wang J. M., 2010, ApJ, 710, 836
  * Rapisarda et al. (2014) Rapisarda S., Ingram A., van der Klis M., 2014, MNRAS, 440, 2882
  * Rapisarda et al. (2016) Rapisarda S., Ingram A., Kalamkar M., van der Klis M., 2016, MNRAS, 462, 4078
  * Timmer & Koenig (1995) Timmer J., Koenig M., 1995, A&A, 300, 707
  * Torres et al. (2019) Torres M. A. P., Casares J., Jiménez-Ibarra F., Muñoz-Darias T., Armas Padilla M., Jonker P. G., Heida M., 2019, ApJ, 882, L21
  * Uttley et al. (2011) Uttley P., Wilkinson T., Cassatella P., Wilms J., Pottschmidt K., Hanke M., Böck M., 2011, MNRAS, 414, L60
  * Wijnands et al. (1999) Wijnands R., Homan J., van der Klis M., 1999, ApJ, 526, L33
  * Zhang et al. (2020) Zhang L., et al., 2020, MNRAS, 494, 1375
  * van der Klis et al. (1987) van der Klis M., Hasinger G., Stella L., Langmeier A., van Paradijs J., Lewin W. H. G., 1987, ApJ, 319, L13


# A simple evaluation of a theta value and the Kronecker limit formula

Fernando Chamizo The author is partially supported by the PID2020-113350GB-I00
grant of the MICINN (Spain) and by “Severo Ochoa Programme for Centres of
Excellence in R&D” (SEV-2015-0554).

(July 18, 2021)

###### Abstract

We evaluate the classic sum $\sum_{n\in\mathbb{Z}}e^{-\pi n^{2}}$. The novelty
of our approach is that it does not require any prior knowledge about modular
forms, elliptic functions or analytic continuations. Even the $\Gamma$
function, in terms of which the result is expressed, only appears as a complex
function in the computation of a real integral by the residue theorem. Another
contribution of this note is to provide a very simple proof of the Kronecker
limit formula.

Keywords: Theta function, Gamma function

Mathematics Subject Classification (2020): Primary 11F67, 11Y60; Secondary
11F27

##  1 Introduction

Our goal is to give a proof of the following result with very few
prerequisites. In general, the special values of theta and allied functions
are related to deep topics in number theory (complex multiplication, class
field theory, modular forms, elliptic functions, etc., cf. [3], [2]) which we
avoid here.

######  Theorem 1.1.

Consider $\theta(z)=\sum_{n\in\mathbb{Z}}e^{\pi in^{2}z}$. Then

$\theta(i)=(2\pi)^{-1/4}\sqrt{\frac{\Gamma(1/4)}{\Gamma(3/4)}}=\frac{\Gamma(1/4)}{\pi^{3/4}\sqrt{2}}$

with $\Gamma$ the classical Gamma function
$\Gamma(s)=\int_{0}^{\infty}t^{s-1}e^{-t}\;dt$.

We will prove the first equality, the second follows from the relation
$\Gamma(s)\Gamma(1-s)=\pi\csc(\pi s)$ that do not use elsewhere. In fact the
$\Gamma$ function only appears as a complex function in the computation of an
integral (Lemma 3.2) and, beyond that, we barely use its defining integral
representation for $s>1$.

Except for a special case of the Jacobi triple product identity and the well
known formula for the number of representations as a sum of two squares (both
separated in §2 and admitting elementary proofs, not included here), the proof
is completely self-contained. The techniques only involve basic real and
complex variable methods. No modular properties of $\theta$ and $\eta$ and no
functional equations of any $L$-function or Eisenstein series nor their
analytic continuations, are required.

Our argument includes a proof of a version of the (first) Kronecker limit
formula (Proposition 3.1) simpler than the ones we have found in the
literature (cf. [8]) which may have independent interest. We address the
reader to the interesting paper [4] for the history and relevance of this
formula.

##  2 Two auxiliary results

We first recall the factorization of the $\theta$ function.

######  Lemma 2.1.

For $|q|<1$,

$\sum_{n=-\infty}^{\infty}q^{n^{2}}=\prod_{n=1}^{\infty}\big{(}1-q^{2n}\big{)}\big{(}1+q^{2n-1}\big{)}^{2}.$

The next result is the classic formula for $r(n)$, the number of
representations of $n$ as a sum of two squares, in terms of the nontrivial
character $\chi$ modulo $4$ (i.e., $\chi(n)=(-1)^{(n-1)/2}$ for $n$ odd and
zero for $n$ even).

######  Lemma 2.2.

For $n\in\mathbb{Z}^{+}$ and $s>1$, we have

$r(n)=4\sum_{d\mid n}\chi(n)\quad\text{or
equivalently,}\quad\sum_{n=1}^{\infty}r(n)n^{-s}=4\zeta(s)L(s)$

with $\zeta(s)=\sum_{n=1}^{\infty}n^{-s}$ the Riemann zeta function and
$L(s)=\sum_{n=1}^{\infty}\chi(n)n^{-s}$.

We will say some words about their proofs.

Lemma 2.1 comes from the Jacobi triple product identity which admits
elementary combinatorial proofs (see [6, §8.3] and [1]) but arguably, even
today, the conceptually most enlightening proof is the classic one based on
complex analysis [10, §10.1]. It uses the invariance by two translations of
certain entire function to conclude that it is a constant, which is computed
with a beautiful argument due to Gauss [9, §78].

Lemma 2.2 can be derived from the triviality of some spaces of modular forms
or from some properties of elliptic functions [10, §10.3.1], [9, §84]. A less
demanding proof, requiring quadratic residues and almost nothing else, is to
use the representations of an integer by the quadratic forms in a class [6,
§12.4]. A longer alternative is to show that $\mathbb{Z}[i]$ is a UFD and
deduce the result from $r(p)=4(1+\chi(p))$ for $p$ prime, which is essentially
Fermat two squares theorem [5, Art.182] (see [11] for a “one-sentence” proof
of the latter).

##  3 The Kronecker limit formula and the theta evaluation

We first state a compact version of the Kronecker limit formula and provide a
proof only requiring the residue theorem and the very easy [7, p.23] and well
known result $(s-1)\zeta(s)\to 1$ as $s\to 1^{+}$.

The Epstein zeta function $\zeta(s,Q)$ associated to a positive definite
binary quadratic form $Q$ and the Dedekind $\eta$ function are defined by:

$\zeta(s,Q)=\sum_{\vec{n}\in\mathbb{Z}^{2}\setminus\\{\vec{0}\\}}\big{(}Q(\vec{n})\big{)}^{-s}\quad\text{
and }\quad\eta(z)=e^{\pi iz/12}\prod_{n=1}^{\infty}\big{(}1-e^{2\pi
inz}\big{)}.$

We assume $s>1$ and $\Im z>0$ to assure the convergence.

######  Proposition 3.1.

Let $Q(x,y)=ax^{2}+bxy+cy^{2}$ be a real form with $D=4ac-b^{2}>0$ and $a>0$.
Then

$\lim_{s\to
1^{+}}\Big{(}\frac{\sqrt{D}}{4\pi}\zeta(s,Q)-\zeta(2s-1)\Big{)}=\log\frac{\sqrt{a/D}}{|\eta(z_{Q})|^{2}}\qquad\text{with}\quad
z_{Q}=\frac{-b+i\sqrt{D}}{2a}.$

###### Proof.

Let $f(s)=-\int_{-\infty}^{\infty}Q(x,1)^{-s}\;dx$. The limit in the statement
equals $L_{1}-L_{2}$ with

$L_{1}=\lim_{s\to
1^{+}}\frac{\sqrt{D}}{4\pi}\Big{(}\zeta(s,Q)+2\zeta(2s-1)f(s)\Big{)},\
L_{2}=\lim_{s\to 1^{+}}\zeta(2s-1)\Big{(}\frac{\sqrt{D}}{2\pi}f(s)+1\Big{)}.$

L’Hôpital’s rule shows $L_{2}=\frac{\sqrt{D}}{4\pi}f^{\prime}(1)$ because
$(2s-2)\zeta(2s-1)\to 1$ (and the residue theorem assures
$f(1)=-2\pi/\sqrt{D}$). Then the result follows if we prove

(3.1) $L_{1}=-\log|\eta(z_{Q})|^{2}\qquad\text{and}\qquad
f^{\prime}(1)=-\frac{4\pi}{\sqrt{D}}\log\sqrt{\frac{a}{D}}.$

We have $f^{\prime}(1)=\int_{-\infty}^{\infty}(\log p)/p$ with
$p(x)=ax^{2}+bx+c$. With the change of variables $2ax+b=\sqrt{D}\tan(t/2)$ we
obtain

$f^{\prime}(1)=-\frac{2}{\sqrt{D}}\int_{-\pi}^{\pi}\log\frac{2|\cos(t/2)|}{\sqrt{D/a}}=-\frac{2}{\sqrt{D}}\Re\int_{C}\log\Big{(}\frac{1+z}{\sqrt{D/a}}\Big{)}\frac{dz}{iz}$

with $C$ the unit circle, where we have used
$\log\big{(}2|\cos(t/2)|\big{)}=\Re\log(1+z)$ with $z=e^{it}$. Cauchy’s
integral formula gives the second identity in (3.1).

When we sum $Q(m,n)^{-s}$ the contribution of $n=0$ is $2a^{-s}\zeta(2s)$. Let
$g_{s}(z)=Q(z,1)^{-s}+Q(z,-1)^{-s}$. For $n\neq 0$,
$Q(m,n)=n^{2}Q(m/n,1)=n^{2}Q(-m/n,-1)$. Then by the residue theorem in the
band $B_{\epsilon}=\\{|\Im z|<\epsilon\\}$ with $0<\epsilon<\Im z_{Q}$,

$\zeta(s,Q)-2\frac{\zeta(2s)}{a^{s}}=\sum_{n=1}^{\infty}\frac{1}{n^{2s}}\sum_{m\in\mathbb{Z}}g_{s}\big{(}\frac{m}{n}\big{)}=\sum_{n=1}^{\infty}\frac{-1}{2n^{2s-1}}\int_{\partial
B_{\epsilon}}g_{s}(z)i\cot(\pi nz)\;dz.$

As $g_{s}$ is even, $\int_{\partial B_{\epsilon}}=-2\int_{L_{\epsilon}}$ with
$L_{\epsilon}=\\{\Im z=\epsilon\\}$ oriented to the right and the sum is
$\sum_{n}n^{1-2s}\int_{L_{\epsilon}}$. Note that
$\int_{L_{\epsilon}}g_{s}=\int_{L_{0}}g_{s}=-2f(s)$. Then adding
$2\zeta(2s-1)f(s)$ is equivalent to replace $i\cot(\pi nz)$ by $i\cot(\pi
nz)-1$ in $\int_{L_{\epsilon}}$. The expansion $i\cot
w-1=2e^{2iw}/(1-e^{2iw})=2(e^{2iw}+e^{4iw}+\dots)$ assures an exponential
decay and we have

$L_{1}=\frac{\sqrt{D}}{4\pi}\Big{(}2\frac{\zeta(2)}{a}+\sum_{n,k=1}^{\infty}\frac{2}{n}\int_{L_{\epsilon}}g_{1}(z)e^{2\pi
inkz}\;dz\Big{)}.$

Substitute $\zeta(2)=\pi^{2}/6$ and note that
$g_{1}(z)=\big{(}a(z-z_{Q})(z-\bar{z}_{Q})\big{)}^{-1}+\big{(}a(z+z_{Q})(z+\bar{z}_{Q})\big{)}^{-1}$.
The residue theorem in $\\{\Im z>\epsilon\\}$ gives promptly

$L_{1}=\frac{\pi\sqrt{D}}{12a}+\sum_{n,k=1}^{\infty}\frac{1}{n}\big{(}e^{2\pi
nkiz_{Q}}+e^{-2\pi
nki\bar{z}_{Q}}\big{)}=\frac{\pi\sqrt{D}}{12a}-\sum_{k=1}^{\infty}\log\big{|}1-e^{2\pi
kiz_{Q}}\big{|}^{2}$

where the second equality comes from
$\log(1-w)+\log(1-\bar{w})=\log|1-w|^{2}$. The sum is
$\log\big{(}|\eta(z_{Q})|^{2}|e^{-\pi iz_{Q}/6}|\big{)}$ and the proof of
(3.1) is complete. ∎

The evaluation of an integral will be play a role in the final step of our
proof of Theorem 1.1. We proceed again employing the residue theorem.

######  Lemma 3.2.

Let

$I=\frac{1}{\pi}\int_{0}^{\infty}\frac{\log t}{\cosh
t}\;dt\qquad\text{then}\quad\exp(I)=\frac{\Gamma(3/4)}{\Gamma(1/4)}\sqrt{2\pi}.$

###### Proof.

Consider $f(z)=i\sec(2\pi z)\log\Gamma(1/2+z)$ on the vertical band
$B=\big{\\{}|\Re z|<1/2\big{\\}}$. It defines a meromorphic function (for
certain branch of the logarithm because $\Gamma$ does not vanish) with simple
poles at $z_{\pm}=\pm 1/4$. Clearly the residues satisfy $2\pi
i\text{Res}(f,z_{\pm})=\pm\log\Gamma(1/2+z_{\pm})$. This function is
integrable along $\partial B$ and the residue theorem shows

$\log\frac{\Gamma(3/4)}{\Gamma(1/4)}=\int_{\partial
B}f=\int_{-\infty}^{\infty}\frac{\log\Gamma(1+it)}{\cosh(2\pi
t)}\;dt-\int_{-\infty}^{\infty}\frac{\log\Gamma(it)}{\cosh(2\pi t)}\;dt.$

Using $\Gamma(1+it)=it\Gamma(it)$ and taking real parts to avoid
considerations about the branch of the logarithm,

$\log\frac{\Gamma(3/4)}{\Gamma(1/4)}=\int_{-\infty}^{\infty}\frac{\log|t|}{\cosh(2\pi
t)}\;dt=\frac{1}{\pi}\int_{0}^{\infty}\frac{\log(t/2\pi)}{\cosh
t}\;dt=I-\int_{0}^{\infty}\frac{\log(2\pi)}{\pi\cosh t}\;dt.$

The last integral is $\log\sqrt{2\pi}$ just changing $t=\log u$. ∎

###### Proof of Theorem 1.1.

The identity $\theta(z)=\prod_{n=1}^{\infty}\big{(}1-e^{2\pi
inz}\big{)}\big{(}1+e^{\pi i(2n-1)z}\big{)}^{2}$ follows from Lemma 2.1 with
$q=e^{\pi iz}$. Some elementary manipulations with the definition of $\eta$
show $\theta(z)=\eta^{2}\big{(}\frac{1}{2}z+\frac{1}{2}\big{)}/\eta(z+1)$. Let
$Q=x^{2}+y^{2}$ and $Q^{\prime}=2x^{2}-2xy+y^{2}$ with $z_{Q}=i$ and
$z_{Q^{\prime}}=\frac{1+i}{2}$. We have $\zeta(s,Q)=\zeta(s,Q^{\prime})$
because $Q^{\prime}=x^{2}+(x-y)^{2}$. Then Proposition 3.1 implies
$|\eta(z_{Q^{\prime}})|^{2}/|\eta(z_{Q})|^{2}=\sqrt{2}$ and, noting
$\theta(i)=|\theta(i)|$ and $|\eta(z+1)|=|\eta(z)|$,

$\theta(i)=\theta(z_{Q})=\Big{|}\frac{\eta(z_{Q^{\prime}})}{\eta(z_{Q}+1)}\Big{|}^{2}|\eta(z_{Q}+1)|=\sqrt{2}|\eta(z_{Q})|.$

Recalling Lemma 3.2, Theorem 1.1 is equivalent to
$I=-\log\big{(}2|\eta(z_{Q})|^{2}\big{)}$. By Lemma 2.2, we have
$\zeta(s,Q)=4\zeta(s)L(s)$ and, since Proposition 3.1, we must prove

$\lim_{s\to 1^{+}}\Big{(}\frac{2}{\pi}\zeta(s)L(s)-\zeta(2s-1)\Big{)}=I.$

It is known $\zeta(s)\sim(s-1)^{-1}+\gamma$ as $s\to 1$ with $\gamma$ the
Euler-Mascheroni constant, and it admits a short elementary proof [7, p.23].
Using111For a quick proof write
$\Gamma(s)=\lim_{n\to\infty}\int_{0}^{n}\big{(}1-\frac{x}{n}\big{)}^{n}x^{s-1}\;dx$
to obtain, by repeated partial integration,
$\lim\frac{n!n^{s}}{s(s+1)\cdots(s+n)}$ (Gauss’ definition of $\Gamma$). The
derivative of its logarithm at $s=1$ gives finally
$\Gamma^{\prime}(1)=-\gamma=\lim\big{(}\log
n-\frac{1}{1}-\frac{1}{2}-\cdots-\frac{1}{n}\big{)}$.
$\Gamma^{\prime}(1)=-\gamma$ we have $\zeta(s)-2\Gamma(s)\zeta(2s-1)\to 0$ and
the previous limit is

$\lim_{s\to
1^{+}}\Big{(}\frac{4}{\pi}\Gamma(s)L(s)-1\Big{)}\zeta(2s-1)=\lim_{s\to
1^{+}}\frac{4\Gamma(s)L(s)-\pi}{2\pi(s-1)}=\frac{2}{\pi}\frac{d}{ds}\Big{|}_{s=1}\big{(}\Gamma(s)L(s)\big{)}$

by L’Hôpital’s rule. It only remains to show that this derivative is $\pi
I/2$. Plainly $\Gamma(s)n^{-s}=\int_{0}^{\infty}t^{s-1}e^{-nt}\;dt$. Then

$\Gamma(s)L(s)=\int_{0}^{\infty}t^{s-1}\big{(}e^{-t}-e^{-3t}+e^{-5t}-e^{-7t}+\dots\big{)}\;dt=\int_{0}^{\infty}\frac{t^{s-1}}{2\cosh
t}\;dt$

and Lemma 3.2 implies the result differentiating under the integral sign. ∎

### Acknowledgments

I am deeply indebted to E. Valenti.

## References

  * [1] G. E. Andrews.  A simple proof of Jacobi’s triple product identity.  Proc. Amer. Math. Soc., 16:333–334, 1965.
  * [2] F. Chamizo and D. Raboso.  Modular forms and almost integers _(Spanish)_.  Gac. R. Soc. Mat. Esp., 13(3):539–555, 2010.
  * [3] D. A. Cox.  Primes of the form $x^{2}+ny^{2}$.  A Wiley-Interscience Publication. John Wiley & Sons, Inc., New York, 1989\.  Fermat, class field theory and complex multiplication.
  * [4] W. Duke, Ö. Imamoḡlu, and Á. Tóth.  Kronecker’s first limit formula, revisited.  Res. Math. Sci., 5(2):Paper No. 20, 21, 2018.
  * [5] C. F. Gauss.  Disquisitiones arithmeticae.  Springer-Verlag, New York, 1986.  Translated and with a preface by Arthur A. Clarke, Revised by W. C. Waterhouse, C. Greither and A. W. Grootendorst and with a preface by Waterhouse.
  * [6] L. K. Hua.  Introduction to number theory.  Springer-Verlag, Berlin-New York, 1982.  Translated from the Chinese by P. Shiu.
  * [7] H. Iwaniec.  Lectures on the Riemann zeta function, volume 62 of University Lecture Series.  American Mathematical Society, Providence, RI, 2014.
  * [8] Y. Motohashi.  A new proof of the limit formula of Kronecker.  Proc. Japan Acad., 44:614–616, 1968.
  * [9] H. Rademacher.  Topics in analytic number theory.  Die Grundlehren der mathematischen Wissenschaften, Band 169. Springer-Verlag, New York-Heidelberg, 1973.  Edited by E. Grosswald, J. Lehner and M. Newman.
  * [10] E. M. Stein and R. Shakarchi.  Complex analysis, volume 2 of Princeton Lectures in Analysis.  Princeton University Press, Princeton, NJ, 2003.
  * [11] D. Zagier.  A one-sentence proof that every prime $p\equiv 1\pmod{4}$ is a sum of two squares.  Amer. Math. Monthly, 97(2):144, 1990.

Departamento de Matemáticas and ICMAT, Universidad Autónoma de Madrid, 28049
Madrid, Spain

Email address<EMAIL_ADDRESS>


#  Shape the Future of ITS –
Optimization and scheduling for a large scale urban transportation system –
in a fast-changing world

Yi Zhang Yi Zhang is affiliated with with the Institute for Infocomm Research
(I2R), Agency for Science, Technology and Research (A*STAR), Singapore 138632.
Emails: yzhang120@e.ntu.edu.sg(zhang$\textunderscore$yi@i2r.a-star.edu.sg).

##  I Introduction

The rapid growth of the population facilitates the urban sprawl, automobile
production and leads to heavy traffic congestions, pollution, noises and
traffic fatalities. Therefore, it is indispensable to develop an eco-friendly
sustainable transport with high performance of the existing traffic network.
Although many transformational products have been developed and used on the
road, such as autonomous vehicles (AV) and drones, and their quality is
continually upgraded, 100% replacement of traditional cars to AVs requires not
only the safety and security evaluation but also corresponding physical
infrastructures, e.g., non-signalized intersection, charging stations. Also,
we still have a long way to go in order to completely commercialize AVs with
level 5 automation. Thus, it is predictable that traffic signals will still
play an important role for a certain long time in this hybrid transportation
world, which involves both human-driven cars and AVs under different automated
levels. At the period of transition, a feasible and implementable approach is
more helpful and practical to reshape the transportation system in the near
future.

The concept of the pedestrian/transit-oriented transport, which is designed
for the human body instead of the car body, has been proposed by researchers
[2]: “cities with very high walking, cycling and transit mode share (i.e., 75%
or more) typically have high density, mixed use urban centres at or above
100-200 people per hectare and are supported by a transportation strategy that
prioritizes pedestrians first, then cyclists followed by transit users.”
However, vehicle traffic still gets excessively high attention in many cities’
guidelines, and walkability is seldom taken into account. Also, walking is the
access to public transport (PT), thereby PT services could accordingly be
promoted if the pedestrian walking safety and pleasantness can be satisfied.
PT system has been studied for several decades due to its large ridership and
sustainability on economic efficiency, environmental protection and social
equity. Proper bus dispatching and operation can attract more passengers,
which encourages commuters to change their travel mode from private
automobiles to public buses to further alleviate the traffic congestion and
air pollution.

However, PT system, a typical ride-sharing transport serving a variety of
access needs and providing an equal social value, now faces the challenges
from the mobility-on-demand (MOD) system, which is famed for its easy
transactions and convenient access via mobile phones. Therefore, we need to
explore and create a win-win cooperation model between both PT and AV systems
at current transitional period. The flexible demand-driven pattern in AVs
could make up the shortages due to fixed routes in the PT system. On the other
hand, PT system is still coping with the high-volume transfer tasks, which may
not be feasible for AVs due to the safety concern.

The objective of this essay is to propose a set of coordinated technological
solutions to transform existing transport system to a more intelligent
interactive system by adopting optimization and control methods implementable
in the near future, thereby improving public services and quality of life for
residents. In this essay, three different application scenes that closely
related to people’s daily life are discussed. We firstly propose a traffic
light scheduling strategy via model predictive control (MPC) method, with the
aim to fairly minimize both pedestrians’ and vehicles’ delay. After that, a
combined dispatching-operation system is proposed to further increase the
control flexibility, and corresponding implementation solution for boarding
control is also illustrated. Finally, a possible scheme to combine both PT and
AV systems is proposed to improve existing PT system.

##  II Working Packages and Tasks

###  II-A Traffic light scheduling for vehicle-pedestrian mixed-flow networks

Most traffic signal controllers differentiate vehicles and pedestrians, and
focus on their great contributions to vehicle flows, this is reasonable when
pedestrian volume is low. However, in the downtown areas where large number of
pedestrians interfere with vehicular traffic, optimizing traffic signals only
for vehicles may create more conflicts between both traffic participants and
potentially reduce the economic interest, since pedestrians in CBD areas are
usually potential customers of nearby shopping malls. According to SGS
Economics and Planning [1], optimized pedestrian flow can increase $1.3
billion a year for Melbourne CBD area. On the other hand, rich data can be
obtained with the help of the advanced sensor equipments, such as V2X, 5G
communication, Lidar and so on. Powerful machine learning algorithms can be
effectively utilized to help predict the traffic flows based on the large
amount of database, which better serve the optimized signal controller.

In view of this, we propose this signal controller with the aim to fairly
minimize both vehicles’ and pedestrians’ delay [4][7]. Fig. 1 illustrates the
framework of the proposed real-time traffic light scheduling strategy
implemented in simulator software VISSIM. The macroscopic flow model for both
pedestrians and vehicles are developed, and the impacts of the signal on
pedestrian crossing capacity is captured in the model. Benefitted from the
advanced sensors, signal priority level could also be incorporated into the
model to enable higher priority for public buses. The mixed-flow model is then
solved by adopting the commercial optimization solver or evolutionary
algorithms. After that, the optimized signal phases and duration are sent to
the simulator VISSIM, which mimics the real urban traffic environment.
Meanwhile, the traffic information, such as traffic volumes, turning ratios,
is stored in the back-end database, which is used to re-train the AI models by
adopting the machine learning algorithms to predict the required information.
The predicted traffic parameters and current real-time information are all
sent to the controller side, which is solved in a rolling horizon manner, and
the whole process continues with the evolve of the traffic system to form a
closed-loop control strategy.

Figure 1: Framework of traffic light scheduling for vehicle-pedestrian mixed-
flow networks

###  II-B A combined dispatching-operation strategy for public bus management
incorporated with the boarding control

The bus operation systems in most today’s studies are proposed based on known
bus dispatching time or schedule headways/frequencies. On the other hand, the
bus dispatching systems seldom consider operation control for on-road buses.
It is understandable that separately deal with these two problems could reduce
the computational complexity significantly compared with the complicated
combined model. However, the increasing application and perfection of
telematics (e.g., Automatic Vehicle Location (AVL), Automatic Passenger
Counting (APC), etc.) in the bus management system, the collection of the
real-time information becomes possible, also, the constantly upgraded computer
with more computing power is keeping breaking the records, thus, it is
predictable that to solve a combined optimization problem is implementable
[5]. Decision variables, such as bus dispatching time, bus speed between any
two adjacent stops, bus dwell time at each stop, OD-based boarding volume,
could all be described into a holistic model to enable this combined
dispatching-operation bus management system [6]. Also, boarding control has
been studied in some literature and proved its high efficiency in improving
the bus service quality. However, it still has not been widely implemented in
the real world.

Fig. 2 gives an graph illustration of future bus stop, which makes the
boarding control applicable. Imagine a message via mobile apps or screen board
located at the bus stop could be sent or displayed to passengers when a bus is
approaching the stop, and only passengers at the front queue will be selected
to enter the designated boarding area at the bus stop. The identification of
the front-queue passenger can be realized by the camera or Lidar installed at
the bus stop. When the bus reaches the bus stop, only passengers who are
waiting at the designated area are allowed to board the bus. This requires
infrastructural enhancement at the bus stop, not only the area re-design but
also the advanced sensors.

Figure 2: Graphic illustration of future bus stop

###  II-C Autonomous bus fleet management for mobility-on-demand system

To increase the vehicle utilization and reduce the carbon dioxide emissions,
the on-line carsharing services have been carried out in many cities in the
past few years. Also, the adoption of the electric vehicles could potentially
reduce emissions to promote a sustainable environment, meanwhile, the forth-
coming commercialization of AVs has constantly reshaping the transportation
system. All of the above emerging concepts and technologies have combined and
led to an autonomous mobility-on-demand system (AMOD). AMOD system is not a
new concept and has been proposed since 2014 [3], however, the implementation
and research studies about AMOD system mainly focus on private cars and taxi
services, and this may bring challenges on the traditional public transport,
which has been regarded as the symbol of equity and accessibility. Therefore,
it is strongly necessary to identify and explore the synergistic possibilities
between AMOD and public transports, especially under this transitional period.

Fig. 3 illustrates one possible case capturing how AVs could support the
existing public bus lines. The left subfigure describes two traditional bus
lines served by corresponding buses dispatched from the terminal. The orange
bar represents the demand size at each bus stop. All control variables
mentioned in Section II-B could be adopted and tuned to facilitate the
operation of the traditional bus service. Meanwhile, we could also dispatch
AVs to serve bus stops on multiple bus lines but with high demands due to its
flexible on-demand characteristics. AVs are normally commercialized as
electric cars with limited battery capacity, thus, it is inevitable to develop
a smart AV routing strategy catering for the demands of the users as well as
the battery capacity level, as illustrated in Fig. 4. As AVs do not follow a
fixed bus route and on-board passengers need to be sent to their destined
stops, thus, on-demand autonomous bus fleet dispatching strategy inevitably
requires the boarding control at each bus stop, and this enables the system to
consider not only the routing decision but also the volume dynamics, which is
different with current one passenger pick-up and drop-off problem. The
optimization involving in routing, subsequent bus stop selection and
corresponding boarding volume makes the problem more challenging but also
interesting.

Figure 3: Application of autonomous buses on supporting traditional public bus
lines Figure 4: Autonomous bus routing for multi-loading and on-road charging

##  III Conclusion

In this essay, technological solutions from a pedestrian/transit-oriented
perspective have been provided. Firstly, an adaptive traffic signal control
framework is presented, with consideration of a macroscopic mixed-flow
optimization model, VISSIM simulation platform, historical database and
machine learning-enabled AI prediction model. The experiment results in our
paper [4] have demonstrated that the proposed strategy in a Manhattan-shaped
network can strike a good balance between pedestrians’ and vehicle drivers’
needs. Therefore, we are optimistic to implement this strategy in the real
world to enhance the existing traffic signals in the near future. Next, the
bus dispatching and its on-road operation, especially the boarding control,
are combined together to minimize the passenger delay time and the operating
bus vacancy, which makes our strategy more flexible and adaptive to meet the
passenger demand. In our previous study [6], a multi-bus dispatching and
boarding control strategy can reduce roughly 50% of remaining passenger
volumes when compared with the timetable-based fixed schedule, which is quite
promising and exciting. Moreover, to achieve a synergistic objective between
PT system and AV system, a scheme is proposed, where AVs are dispatched to
pick up passengers at high-demand stops to support the PT lines. This scheme
shall be further studied in my future work to incorporate multiple optimized
variables, including AV routing, volume dynamics and boarding limits. Overall,
the essay presents three application scenes in the intelligent transportation
field. The suggested strategies and guidelines can assist in developing an
intelligent future urban transport system in order to provide citizens a
smarter, safer and more interactive transportation experience.

## References

  * [1] SGS Economics and Planning.  CBD pedestrian analysis, technical report city of melbourn.  Technical report, 2014.
  * [2] Jennie Moore, Kirstin Mille, Richard Registe, Sarah Campbell, Julian Zelazny, and Sven Eberlein.  International ecocity standards(brochure).  Technical report, 2017.
  * [3] Kevin Spieser, Kyle Treleaven, Rick Zhang, Emilio Frazzoli, Daniel Morton, and Marco Pavone.  Toward a systematic approach to the design and evaluation of automated mobility-on-demand systems: A case study in singapore.  In Road vehicle automation, pages 229–245. Springer, 2014.
  * [4] Yi Zhang, Kaizhou Gao, Yicheng Zhang, and Rong Su.  Traffic light scheduling for pedestrian-vehicle mixed-flow networks.  IEEE Transactions on Intelligent Transportation Systems, 20(4):1468–1483, 2018.
  * [5] Yi Zhang, Rong Su, Yicheng Zhang, and Nadeesha Sandamali Gammana Guruge.  A multi-bus dispatching strategy based on boarding control.  IEEE Transactions on Intelligent Transportation Systems, 2021.
  * [6] Yi Zhang, Rong Su, Yicheng Zhang, and Bohui Wang.  Dynamic multi-bus dispatching strategy with boarding and holding control for passenger delay alleviation and schedule reliability: A combined dispatching-operation system.  IEEE Transactions on Intelligent Transportation Systems, 2021.
  * [7] Yi Zhang, Yicheng Zhang, and Rong Su.  Pedestrian-safety-aware traffic light control strategy for urban traffic congestion alleviation.  IEEE Transactions on Intelligent Transportation Systems, 22(1):178–193, 2019.

|  Yi Zhang (S’17-M’21) Yi Zhang received her Bachelor degree of Engineering
from Shandong University, China in 2014, and the PhD degree in Electrical and
Electronic Engineering from Nanyang Technological University, Singapore in
2020. She is currently a research scientist at the Institute for Infocomm
Research (I2R) in the Agency for Science, Technology and Research, Singapore
(A*STAR). Her research interests focus on intelligent transportation system,
including urban traffic flow management, model-based traffic signal
scheduling, lane change prediction and bus dispatching and operation
management.
---|---