2019 Vol. X No. XX, 000–000

11institutetext: School of Astronomy and Space Science, Nanjing University,
Nanjing 210093, China<EMAIL_ADDRESS>
22institutetext: The School of Science, Langfang Normal University, Langfang
065000, China
33institutetext: Key Laboratory of Modern Astronomy and Astrophysics, Nanjing
University, Nanjing 210023, China
Received ********; accepted ********

#  The mid-infrared variability of the SDSS optical quasars

Hongtao Wang 112233 Yong Shi† 1133

###### Abstract

Based on the Seventh Data Release (DR7) quasar catalog from the Sloan Digital
Sky Survey, we investigate the variability of optical quasars in W1, W2, W3
and W4 bands of the Wide-field Infrared Survey Explorer (WISE) and the Near-
Earth Object Wide-field Infrared Survey Explorer (NEOWISE). Adopting the
structure function method, we calculate the structure function ($\rm\delta
t$=1 yr) which shows no obvious correlations with the bolometric luminosity,
the black hole mass and the Eddington ratio. The ensemble structure functions
in W1 and W2 bands show that the SF slopes are steeper than those in previous
studies which may be caused by different cadence and observational epoch
number. We further investigate the relation of variability amplitude
$\sigma_{m}$ between mid-infrared band and optical band, but no obvious
correlation is found. No correlation is found between W1-W2 and g-r color. We
think the mid-infrared emission of quasars may be smoothed out by the extended
dust distribution, thus leading to no obvious correlation. For the radio-loud
quasar sub-sample, we further analyze the relation between the variability
amplitude in the mid-infrared band and the radio luminosity at 6 cm, but no
obvious correlations are found, which indicate the mid-infrared emission
contributed from the synchrotron radiation of the relativistic jet is very
weak.

###### keywords:

quasar: variability–quasar: mid-IR–quasar: characteristics

##  1 Introduction

Active galactic nuclei (AGNs) are powered by the supermassive black holes
(SMBHs) that accrete matter and release enormous energy (Rees 1984).
Surrounding the central SMBH, there is usually an accretion disc and a dusty
torus (Urry & Padovani 1995), which account for the observed big blue bump
(BBB) and mid-infrared bump in their spectral energy distribution (SED),
respectively. The BBB is considered as the thermal emission from the accretion
disc (Shakura & Sunyaev 1973) and the mid-infrared bump is believed to
originate from the dusty torus (Koshida et al. 2014; Mandal et al. 2018).

Quasars have long been recognized as intrinsically variable sources since
their discovery ((Matthews & Sandage 1963; Wagner & Witzel 1995; Ulrich,
Maraschi & Urry 1997). The aperiodic flux variability of AGNs is ubiquitous
from radio to high-energy $\gamma$-ray, and the timescales vary from minutes
to years (Ulrich et al. 1997). The variability provides us many invaluable
clues to study AGNs, such as estimating the black hole mass (Peterson et al.
2004; Vestergaard & Peterson 2006 and reference therein), constraining the
structure of broad line regions (Pancoast et al. 2011) and identifying AGNs
(Padovani et al. 2017) etc.

In mid-infrared (mid-IR) band, the emission is dominated by the dusty torus
(see, e.g., Nenkova et al. 2008), which reprocesses the UV/optical emission
from accretion disks into the mid-IR emission. Non-thermal emission from jets
sometimes also contributes to some fractions of the mid-IR luminosity in
radio-loud AGNs. Based on 4$\sim$5 epoch observations, Kozłowski et al. (2016)
analyzed the ensemble mid-IR variability of $\sim$ 1500 quasars by the single
power-law structure function (SF), and found the SF slope in the mid-IR band
was steeper than that in the optical band. Moreover, the SF slope show no
correlations with the luminosity in mid-IR band and rest-frame wavelength,
while the amplitude does show a correlation with the rest-frame wavelength and
an anti-correlation with the mid-IR luminosity.

With the coming of WISE, there is a new opportunity to carry out studies of
mid-IR variability. The WISE mission is divided into two main cycles. The
first stage is the full solid hydrogen cryogenic phase with W1, W2, W3 and W4
channels and the asteroid-hunting phase with only W1 and W2 channels from
January 2010 to February 2011. All observed data sets during this stage are
included in the ‘AllWISE MultiEpoch Photometry (MEP) Database’. The second
stage is near earth object WISE Reactivation mission during 2014 to 2018 with
W1 and W2 bands. The data sets in this stage are included in the ‘NEOWISE-R
Single Exposure (L1b) Source Table’. The total combined baseline is $\sim$ 9
years in W1 and W2 bands, which offers a chance to investigate the variability
characteristics of quasars in mid-infrared band. Comparing with Spitzer
studies (Kozłowski, et al. 2016 and reference therein), WISE/NEOWISE offers
much higher cadence and larger sample size which provide a great opportunity
for us to systematically investigate the mid-infrared variability of quasars.

The paper is organized as follows. The data sets and sample selections are
described in Section 2. The analytical method and main results are presented
in Section 3. The discussion and conclusion of this work are given in Section
4 and Section 5, respectively. Throughout this paper, we assume a flat
$\Lambda$CDM cosmology with $\rm\Omega_{M}=0.30$, $\rm\Omega_{\lambda}=0.70$
and $\rm H_{0}=70.0~{}km~{}s^{-1}Mpc^{-1}$.

##  2 Sample and Data

In order to investigate the variability of quasars in mid-infrared band, we
utilize the DR7 quasar sample from the Sloan Digital Sky Survey (Shen et al.
2011), including 105,783 quasars brighter than $\rm M_{i}=-22.0$ mag. This
sample has at least one broad emission line with the FWHM larger than $\rm
1000~{}km~{}s^{-1}$, and also contains many available parameters, e.g.,
redshift, black hole mass, Eddington ratio and luminosity in 5100 Å etc. We
cross-match the SDSS DR7 quasars catalogue with the mid-IR catalogue, ‘AllWISE
MultiEpoch Photometry Database’ and ‘NEOWISE-R Single Exposure (L1b) Source
Database’(update to April 2019) using a search radius of 5 arcsec. All the
105,783 SDSS quasars have WISE counterparts. The angular resolution is 6.1,
6.4, 6.5and 12.0in W1,W2,W3 and W4 bands, respectively. Our search radius of
5is thus in general larger than the position accuracy of WISE sources (angular
resolution divided by the signal-to-noise). We exclude the epoch photometry
without errors and merge observations within one day.

Some bad photometric measurements are further excluded in W1,W2,W3 and W4
bands for these objects with the criterias as follows.

  1. 1.

The number of point spread function (PSF) components in the profile fit for a
object (nb) must be less than 3, and the best quality single-exposure image (
qi fact=1) are not actively deblended (na$=$0).

  2. 2.

The signal noise ratio (SNR) is larger than 10. The reduced $\chi^{2}$ of
single-exposure profile-fit is less than $5$ in both W1 and W2 bands (i.e.,
w1rchi2$<5$ and w2rchi2$<5$).

  3. 3.

Both the contamination and confusion flag (cc flags) are labeled as 0000. The
angular distance from the nominal boundaries of the South Atlantic Anomaly
(saa sep) is larger than 5.0. Moon masking flag (moon masked) is setted as 00
to avoid contamination by scattered moonlight. Overall frameset quality score
(qual frame) is larger than 0 to ensure the framesets with good quality.

We also apply a 5$\sigma$ criteria ($\sigma$ is the standard deviation of the
light curve) to clip each light curve for removing the photometric outliers
which may be caused by cosmic-rays. Finally, 103,762 objects at least two
epochs in W1 band and W2 band are remain to investigate the long-term
variability. The distribution of bolometric luminosity and redshift of the
final sample is shown in Figure 1. In our final sample, the minimum, median
and maximum time span of the observation is $\sim$ 0.5 years, 4.5 years and
8.5 years, respectively. The minimum, median and maximum epoch number are 2, 6
and 143, respectively. We also utilize the same criteria to construct the
light curve in W3 and W4 bands but with much shorter baseline which is 9
months and 7 months, respectively.

Figure 1:  The distribution of bolometric luminosity and redshift of the final
quasar sample.

.

##  3  The methods to characterize the variability

###  3.1  Structure Function Method

The structure function (SF) method is a typical technique to analyze the
variability of AGNs. It calculates the rms amplitude as a function of time lag
$\delta t$ between different observations. SF is not affected by windows
problem and aliasing due to sparse sampling (Hughes et al. 1992), and suitable
for irregular data sets. Among different definitions of SFs (Simonetti, et al.
1985; Vanden Berk et al. 2004; Sumi et al. 2005; Bauer et al. 2009; MacLeod et
al. 2012), we utilize the one by Kozłowski et al. (2016)

$SF_{\rm obs}(\delta t)=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y(t)-y(t+\delta
t))^{2}},\\\ $ (1)

where $y(t)$ is the observed flux at time t, N is the total observation pairs
in the time lag bin $\delta t$. SF is calculated by three processes as
follows. First, we calculate the time lag and magnitude variation of each
pairs in the data sets, and sort the results by time lag. Second, we correct
the time lag in observed frame $|t_{i}-t_{j}|$ into rest frame by $\delta
t=|t_{i}-t_{j}|(1+z)^{-1}$. Third, we divide the $\delta t$ into several bins
and calculate the SF in each bin by expression (1). The bin size is depend on
the number in the data sets. The target is to obtain the well-distributed SF
as far as possible. Emmanoulopoulos et al.(2010) pointed out that the limited
time sampling of individual object could give rise to spurious features in the
SF result, but it could be relieved by the ensemble SF (Sun, et al. 2015). The
SF in each bin is calculated by the bootstrap method 1000 resampling and
choose the $1\sigma$ value as the error of SF.

###  3.2  The variability amplitude

Besides the SF, the variability amplitude ($\sigma_{m}$) is also usually used
to investigate the variability of quasars. The amplitude measures the variance
of the light curve after correcting the measurement uncertainty (Sesar et al.
2007; Ai et al.2010; Ai et al.2013; Rakshit & Stalin 2017; Rakshit et al.
2019). The $\sigma_{m}$ is calculated as the following processes. The standard
deviation ($\Sigma$) of the observed signal in the data sets is first given by

$\Sigma=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(m_{i}-<m>)^{2}},$ (2)

in which $m_{i}$ is the magnitude of $i$-th epoch, and $<m>$ is the mean
amplitude, n is the number of detections. The error $\epsilon$ is then
calculated from the individual error with the following expression

$\epsilon^{2}=\frac{1}{N}\sum_{i=i}^{N}{\epsilon_{i}^{2}+\epsilon_{s}^{2}}.$
(3)

Here $\epsilon_{i}$ is the uncertainty of $i$-th measurement value and
$\epsilon_{s}$ is the systematic uncertainty. Jarrett et al. (2011) showed the
systematic uncertainties of 0.024 mag and 0.028 mag in $\rm W$1 and $\rm W$2,
respectively. At last, the variability amplitude can be written as

$\sigma_{m}={\sqrt{\Sigma^{2}-\epsilon^{2}},\quad\Sigma>\epsilon},$ (4)

otherwise, the $\sigma_{m}=0$.

##  4  The results

###  4.1  The ensemble structure function in W1,W2,W3 and W4 bands

We constructed the ensemble structure functions in four WISE bands as shown in
Figure 2. We fit the results by four-parameter expression (3) in Wang&
Shi(2019) by the least square method, then obtain the slope $\beta=$ $2.75\pm
0.19$ and $2.90\pm 0.31$, respectively. The $\beta$ values suggest the
ensemble SFs in W1 and W2 bands deviate from the damped random walk model
($\beta=1$). Kozłowski et al.(2016) analyzed the mid-IR variability of 1500
AGNs based on the The Decadal IRAC Boötes Survey (DIBS) and its predecessor
the Spitzer Deep, Wide-Field Survey (SDWFS) by single power-law SF, which
suggests the slope to be $\gamma\approx 0.45$ ($\beta\approx 0.90$). It is
obvious that our SF slopes are steeper than that in Kozłowski et al.(2016).
While each of their objects only have 5 epoch observations, nearly 2/5 of our
objects have more than 10 epochs with on average denser cadence. As verified
by Wang & Shi (2019), the SF slope is sensitive to the sampling of a light
curve. As shown in the lower panel of Figure 2, the SFs in W3 and W4 bands do
not show rising trends due to short time baselines.

Figure 2:  The result of ensemble structure function in W1 (a), W2 (b), W3 (c)
and W4 (d) bands

###  4.2  The relation between the ensemble structure functions in W1 and W2
bands and redshift, bolometric luminosity, Eddington ratio as well as the
black hole mass

In Figure 3, we divide the sample into 10 bins in redshift (the bin size of
redshift = 0.55) and 10 bins in bolometric luminosity (the bin size of
luminosity = 0.45) to have in total 100 cells. We calculate the ensemble
structure function in each cell and use the SF ($\delta t$ = 1 yr) as the
variability amplitude. The timescale at one year is suitable because it is far
away from the noise region in the short time lag and the observations are rich
at this timescale. The number of quasars in each cell is presented in Figure
4. As shown in Figure 3, the SF ($\delta t$ = 1 yr) in the W1 band is about
0.1 mag, which is close the SF at $\delta t$ = 2 yr in Kozłowski et al. (2010,
2016). No correlations are seen for this SF ($\delta t$ = 1 yr) with either
redshift or the bolometric luminosity. The SF ($\delta t$ = 1 yr) in the W2
band show overall similar behavior to that in W1 band.

We also analyze the relation between the SF ($\delta t$ = 1 yr) and the black
hole mass as well as the Eddington ratio. The distribution of the Eddington
ratio and the black hole mass is shown in Figure 5. We also divide it into 10
bins in the black hole mass (the bin size = 0.3) from 7.3 to 10.3 and 10 bins
in the Eddington ratio (the bin size = 0.3) from -2.3 to 0.7. In each cell,
the SF ($\delta t$ = 1 yr) in ensemble structure function is presented as the
variability amplitude. As shown in Figure 6, no correlations are found for the
SF ($\delta t$ = 1 yr) with either Eddington ratio or the black hole mass.

Figure 3:  The analytical results in W1 band (upper panel) and W2 band (lower
panel) by structure function.

.

Figure 4:  The object number in W1 and W2 bands.

.

Figure 5:  The distribution of Eddington ratio and the black hole mass of the
final quasar sample.

.

Figure 6:  The relation between the variability amplitude (the SF ($\delta t$
= 1 yr)) in W1 band (upper panel) and W2 band (lower panel) and Edding ratio
as well as black hole mass.

.

###  4.3  The relation between mid-infrared and optical variability

Figure 7:  The relation between mid-infrared and optical variability. The
upper is the $\sigma_{m}$ in W1 band and the lower is the $\sigma_{m}$ in W2
band.

.

Based on the variability amplitude ($\sigma_{m}$) in Section 3.2, we further
investigate the relation between mid-infrared and optical variability
amplitude. We obtain the optical light curve from the SDSS Stripe 82 data sets
from MacLeod et al. (2012), including $\sim$9000 quasars with $\sim$60 epoch
observations and the time baseline of $\sim$ 9 years which is similar to that
in mid-IR band. We choose r-band data sets to calculate the $\sigma_{m}$ due
to its higher SNR, and more available data sets than that in u, g, i and z
bands. At last, we obtain 3943 quasars with measurements of the variability
amplitude in both W1 band and r band, 2856 quasars with both W2 band and r
band. The correlations of the variability amplitudes between mid-IR and
optical bands are presented in Figure 7. The variability amplitude in W1 band
is in general less than 0.4 mag and shows no relationship with that in the
r-band. Again, the W2 band shows similar behaviors. The most optical variable
objects shows very small amplitudes of mid-IR variability in which the torus
may largely evaporate due to the past high phase, thus the rest dust is far
away from the black hole (Jiang et al.2017,2019). We assume the lag is simply
determined by the torus inner radius ($\rm
R_{sub}=0.5L_{46}^{0.5}(1800K/T_{sub})^{2.6}$pc), if the bolometric luminosity
is $\rm 10^{46}ergs^{-1}$ and the temperature 1800 K, thus the time lag is
about 2 yrs. The optical data sets from Stripe 82 coverage is from September
1998 to November 2007, and the mid-IR data sets is from January 2010 to April
2019. Hence, the short time lag and different time coverage may weaken the
correlation between optical band and mid-infrared band. The lack of apparent
relationships of the variability amplitude between W1/W2 and optical bands may
also be caused by several factors such as the time lag between the optical and
mid-IR emission, the smoothing effect in the mid-IR due to its extended size.

###  4.4  The color relation between mid-infrared band and optical band

In order to further check the relation between mid-infrared band and optical
band, we investigate the average color relation between W1-W2 and g-r. The
color index g-r is from the SDSS Stripe 82 data sets. The result is presented
in Figure 8. Most of the g-r is smaller than 0.5 mag. Most of the W1-W2 is
between 0.8 mag and 1.6 mag. The W1-W2 become flat trend with increasing g-r,
which suggest no obvious correlation between W1-W2 and g-r in the DR7 quasar
sample. Yang et al. (2018) analyzed the color variation of 21 changing-look
AGNs, confirmed bluer-when-brighter trend in optical band and redder-when-
brighter trend in mid-infrared band. No correlation between two bands may be
also caused by the time lag between two bands and the smoothing effect in the
mid-IR band.

Figure 8:  The relationship between mid-infrared color (W1-W2) and optical
color (g-r).

.

###  4.5  The relation between the variability amplitude in mid-infrared and
radio luminosity at 6 cm

In order to further examine whether there are some contributions from
synchrotron emission in the mid-infrared band, we select a sub-sample of
radio-loud quasars (the radio loudness above 10) from the Faint Images of the
Radio Sky at Twenty-Centimeters (FIRST) radio survey (White et al. 1997). The
radio loudness is defined as the ratio of the radio flux density at rest frame
6 cm (1.4 GHz) and the optical flux density at rest frame 2500 Å (Shen et al.
2011). In order to further enhance the reliability, we exclude objects of
$\sigma_{m}=0$ or the errors of $\sigma_{m}$ larger than 0.1 mag in the sub-
sample. At last, we obtain 2949 objects that have variability amplitude
$\sigma_{m}$ in W1 band and 2008 objects in W2 band. The relationships between
luminosity in 1.4 GHz and variability amplitude $\sigma_{m}$ in W1/W2 band are
presented in Figure 9. No correlations are seen for both plots. We further
select the objects between the redshift $0.0<z<0.5$, 660 and 548 objects in W1
band and W2 band, respectively. The results are shown in Figure 10 and also no
correlations are found. Jiang et al. (2012) analyzed the mid-infrared
variability of three radio-loud narrow-line Seyfert 1 galaxies using WISE data
sets, and confirmed that they are blazars with jets close to our line of
sight. Their infrared emission probably contains the synchrotron emission from
the relativistic jet. Mao et al. (2018) showed $\gamma$-ray blazars had larger
mid-IR variability than non-$\gamma$-ray blazars. Bisogni et al. (2019) showed
the infrared bump cannot be induced by synchrotron emission form the blazar-
like objects, and found no difference between radio quiet and radio-loud
objects in the stacked mid-infrared spectral energy distribution. Rakshit et
al. (2019) found a positive correlation between the variability amplitude and
1.4 GHz radio power by analyzing 63 radio-detected Narrow Line Seyfert 1
galaxies. From the SED view, the contribution from jet component in radio loud
quasar is mainly in radio band (larger than $\sim$100 $\rm\mu m$), very few
component in mid-IR band (Shang et al.2011). The synchrotron peak of SED in
radio-loud NLS1 mainly distribute from $10^{12}$ GHz (100 $\mu m$) to
$10^{14}$ GHz (3 $\mu m$)(Abdo et al.(2009); Yao et al.(2015); Paliya et
al.(2019)), then it’s likely to be observed in mid-IR band. In contrast, our
sub-sample is composed of radio-loud quasars, and find no obvious correlations
between mid-infrared variability and radio luminosity, which are different
results of Jiang et al. (2012) and Rakshit et al. (2010). This may be because
of different galaxy types: narrow-line Seyfert in theirs vs. radio-loud
quasars in ours.

Figure 9:  The relation between the luminosity at 1.4 GHz and variation
amplitude in mid-infrared W1 band (upper panel) and W2 band (lower panel).

.

Figure 10:  The relation between the luminosity at 1.4 GHz and variation
amplitude in mid-infrared W1 band (upper panel) and W2 band (lower panel) in
the redshift $0.0<z<0.5$.

.

##  5 Discussion

Kozłowski et al. (2010) investigated the variability of objects with four
epochs from the Spitzer Deep Wide-Field Survey (SDWFS) in 3.6 and 4.5 $\mu$m
bands by structure function method. They found the variability amplitude was
$\sim$0.1 mag on the rest-frame timescale of two years and was anti-correlated
with the luminosity. Kozłowski et al. (2016) further added the fifth epoch
from the follow-up the Decadal IRAC Bo$\rm\ddot{o}$tes Survey (DIBS) of the
SDWFS, extending the time baseline to $\sim$ 10 years, and confirmed the above
anti-correlation. In this work, we find that the SF ($\delta t$ = 1 yr) is
about 0.1 mag, which is consistent with the results in Kozłowski et al. (2010,
2016). But the SF slopes in our result are steeper than that in Kozłowski et
al. (2010, 2016) which may be due to different cadence and number of
observational epoch. We don’t find obvious correlation between the bolometric
luminosity and SF ($\delta t$ = 1 yr) in W1 and W2 bands.

The Eddington ratio is usually considered to be the main driver of variability
in quasars. Kelly (2009) analyzed about 100 radio quiet quasars and obtained
the anti-correlation ($-0.25\pm 0.22$) with the relation of log $\sigma^{2}$
and log L/LEdd. Rakshit& Stalin (2017) used the data sets from Catalina Real
Time Transient Survey to analyze the optical variabilty of 5080 Seyfert
galaxies, and found the variability amplitude ($\sigma_{m}$) was anti-
correlated with logL/LEdd, with the correlation coefficient $-0.16$ and p
value $10^{-30}$. Rakshit et al.(2019) analyzed the mid-variability of 492
NLS1 galaxies using WISE data sets, and found the correlations of $\sigma_{m}$
in W1 and W2 bands were anti-correlated with log L/LEdd with the correlation
coefficient (p value) are $-0.14(10^{-3})$ and $-0.21(10^{-6})$, respectively.
The similar trends with Eddington ratio indicate the variability in mid-
infrared band may be correlated with that in optical band.

However, our results in Section 4.3 and 4.4 indicate that the correlations are
not obvious. The variation in the mid-infrared is smoothed out by the extended
dust distribution, hence the intrinsic relation of variation amplitude
$\sigma_{m}$ between optical and infrared band may be smoothed out.

Rakshit et al.(2019) found a positive correlation between the variability
amplitude $\sigma_{m}$ in mid-IR W1/W2 band and radio luminosity from 63
radio-detected narrow-line Seyfert 1 galaxies. Their work inspires us to
investigate the relation between mid-infrared variability and radio luminosity
for radio-loud quasar sample. The results indicate no obvious correlations. It
means the mid-infrared emission from jet contribution is weak in radio-loud
quasars which is consistent with the view of Bisogni et al. (2019), and
implies that the emission of radio-loud quasars in mid-infrared band is
dominated by the thermal emission from dusty torus.

##  6 Conclusions

In this work, based on the multi-epoch observations from WISE and NEOWISE, we
investigate the variability characteristics of the SDSS DR7 quasars in
W1,W2,W3 and W4 bands. The results are listed as follows.

(1). By means of SF method, we find that SF ($\delta t$ = 1 yr) shows no
obvious correlation with the bolometric luminosity, the black hole mass and
the Eddington ratio. The SF slopes in W1 and W2 bands are steeper than that in
Kozłowski et al.(2016) which may be caused by different cadence and
observational number.

(2). We further investigate the relation of variability amplitude between mid-
IR band and optical band, and no obvious correlation is found. The relation
between W1-W2 and g-r is also investigated and shows no correlation . We think
the mid-infrared emission of quasars may be smoothed out by the extended dust
distribution, thus leading to no obvious correlation.

(3). Based on the radio-loud quasar sub-sample, we analyze the relation
between the variability amplitude in mid-infrared band and the radio
luminosity at 6 cm. No obvious correlations are found, which indicates that
the mid-infrared emission from synchrotron radiation of relativistic jet is
weak.

###### Acknowledgements.

We thank the anonymous referees for helpful suggestions that improved the
manuscript. We would like to thank Ge Xue from Nanjing university give us many
useful help in drawing the figure by Python. We also thank Chen Yongyun and Yu
Xiaoling from Nanjing university for their useful comments and discussion when
preparing this paper. This publication makes use of data products from the
Wide-field Infrared Survey Explorer, which is a joint project of the
University of California, Los Angeles, and the Jet Propulsion
Laboratory/California Institute of Technology, funded by the National
Aeronautics and Space Administration, and also makes use of data products from
the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is
a project of the Jet Propulsion Laboratory/California Institute of Technology.
NEOWISE is funded by the National Aeronautics and Space Administration. This
work is supported by the National Key R$\&$D Program of China (No.
2017YFA0402704, No. 2018YFA0404502), the National Natural Science Foundation
of China (NSFC grants 11733002 and 11773013), the Excellent Youth Foundation
of the Jiangsu Scientific Committee (BK20150014). W.H.T acknowledges support
from the Science and Technology supporting Program in Langfang city
(No.2018011005).

## References

  * [Abdo et al.(2009)] Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 707, L142
  * [Ai et al.(2010)] Ai, Y. L., Yuan, W., Zhou, H. Y., et al. 2010, ApJ, 716, L31
  * [Ai et al.(2013)] Ai, Y. L., Yuan, W., Zhou, H., et al. 2013, AJ, 145, 90
  * [Bauer et al.(2009)] Bauer, A., Baltay, C., Coppi, P., et al. 2009, ApJ, 696, 1241
  * [Bisogni et al.(2019)] Bisogni, S., Lusso, E., Marconi, A., et al. 2019, MNRAS, 485, 1405
  * [Emmanoulopoulos et al.(2010)] Emmanoulopoulos, D., McHardy, I. M., & Uttley, P. 2010, MNRAS, 404, 931
  * [Hughes et al.(1992)] Hughes, P. A., Aller, H. D., & Aller, M. F. 1992, ApJ, 396, 469
  * [Jarrett et al.(2011)] Jarrett, T. H., Cohen, M., Masci, F., et al. 2011, ApJ, 735, 112
  * [Jiang et al.(2012)] Jiang, N., Zhou, H.-Y., Ho, L. C., et al. 2012, ApJ, 759, L31
  * [Jiang et al.(2017)] Jiang, N., Wang, T., Yan, L., et al. 2017, ApJ, 850, 63
  * [Jiang et al.(2019)] Jiang, N., Wang, T., Mou, G., et al. 2019, ApJ, 871, 15
  * [Kelly et al.(2009)] Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
  * [Koshida et al.(2014)] Koshida, S., Minezaki, T., Yoshii, Y., et al. 2014, ApJ, 788, 159
  * [Kozłowski et al.(2010)] Kozłowski, S., Kochanek, C. S., Stern, D., et al. 2010, ApJ, 716, 530
  * [Kozłowski et al.(2016)] Kozłowski, S., Kochanek, C. S., Ashby, M. L. N., et al. 2016, ApJ, 817, 119
  * [MacLeod et al.(2012)] MacLeod, C. L., Ivezić, Ž., Sesar, B., et al. 2012, ApJ, 753, 106
  * [Mandal et al.(2018)] Mandal, A. K., Rakshit, S., Kurian, K. S., et al. 2018, MNRAS, 475, 5330
  * [Matthews & Sandage(1963)] Matthews, T. A., & Sandage, A. R. 1963, ApJ, 138, 30
  * [Mao et al.(2018)] Mao, L., Zhang, X., & Yi, T. 2018, Ap&SS, 363, 167
  * [McHardy et al.(2006)] McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., & Fender, R. P. 2006, Nature, 444, 730
  * [Nenkova et al.(2008)] Nenkova, M., Sirocky, M. M., Nikutta, R., et al. 2008, ApJ, 685, 160
  * [Padovani et al.(2017)] Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&A Rev., 25, 2
  * [Paliya et al.(2019)] Paliya, V. S., Parker, M. L., Jiang, J., et al. 2019, ApJ, 872, 169
  * [Pancoast et al.(2011)] Pancoast, A., Brewer, B. J., & Treu, T. 2011, ApJ, 730, 139
  * [Peterson(2004)] Peterson, B. M. 2004, The Interplay Among Black Holes, Stars and ISM in Galactic Nuclei, 15
  * [Rakshit, & Stalin(2017)] Rakshit, S., & Stalin, C. S. 2017, ApJ, 842, 96
  * [Rakshit et al.(2019)] Rakshit, S., Johnson, A., Stalin, C. S., et al. 2019, MNRAS, 483, 2362
  * [Rees(1984)] Rees, M. J. 1984, ARA&A, 22, 471
  * [Sesar et al.(2007)] Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, AJ, 134, 2236
  * [Shakura & Sunyaev(1973)] Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * [Shang et al.(2011)] Shang, Z., Brotherton, M. S., Wills, B. J., et al. 2011, ApJS, 196, 2
  * [Shen et al.(2011)] Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45
  * [Simonetti et al.(1985)] Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46
  * [Sumi et al.(2005)] Sumi, T., Woźniak, P. R., Eyer, L., et al. 2005, MNRAS, 356, 331
  * [Sun et al.(2015)] Sun, M., Trump, J. R., Shen, Y., et al. 2015, ApJ, 811, 42
  * [Ulrich et al.(1997)] Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445
  * [Urry & Padovani(1995)] Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
  * [Vanden Berk et al.(2004)] Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692
  * [Vestergaard, & Peterson(2006)] Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689
  * [Wagner & Witzel(1995)] Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163
  * [Wang, & Shi(2019)] Wang, H., & Shi, Y. 2019, Ap&SS, 364, 27
  * [White et al.(1997)] White, R. L., Becker, R. H., Gregg, M. D., et al. 1997, American Astronomical Society Meeting Abstracts 191, 103.05
  * [Wright et al.(2010)] Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * [Yang et al.(2018)] Yang, Q., Wu, X.-B., Fan, X., et al. 2018, ApJ, 862, 109
  * [Yao et al.(2015)] Yao, S., Yuan, W., Zhou, H., et al. 2015, MNRAS, 454, L16


KUNS-2885, YITP-21-75

Statistical analysis method for the worldvolume hybrid Monte Carlo algorithm

Masafumi Fukuma1*** E-mail address<EMAIL_ADDRESS>,
Nobuyuki Matsumoto2††† E-mail address<EMAIL_ADDRESS>and Yusuke
Namekawa3‡‡‡ E-mail address<EMAIL_ADDRESS>

1Department of Physics, Kyoto University, Kyoto 606-8502, Japan 2RIKEN/BNL
Research center, Brookhaven National Laboratory, Upton, NY 11973, USA 3Yukawa
Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

We discuss the statistical analysis method for the worldvolume hybrid Monte
Carlo (WV-HMC) algorithm [arXiv:2012.08468], which was recently introduced to
substantially reduce the computational cost of the tempered Lefschetz thimble
method. In the WV-HMC algorithm, the configuration space is a continuous
accumulation (worldvolume) of deformed integration surfaces, and sample
averages are considered for various subregions in the worldvolume. We prove
that, if a sample in the worldvolume is generated as a Markov chain, then the
subsample in the subregion can also be regarded as a Markov chain. This
ensures the application of the standard statistical techniques to the WV-HMC
algorithm. We particularly investigate the autocorrelation times for the
Markov chains in various subregions, and find that there is a linear relation
between the probability to be in a subregion and the autocorrelation time for
the corresponding subsample. We numerically confirm this scaling law for a
chiral random matrix model.

###### Contents

  1. 1 Introduction
  2. 2 Stochastic process in a subregion and the integrated autocorrelation time
    1. 2.1 Stochastic process in a subregion
    2. 2.2 Integrated autocorrelation time for a subchain
    3. 2.3 Scaling law for the integrated autocorrelation times
  3. 3 Numerical confirmation of the scaling law
    1. 3.1 Setup
    2. 3.2 Results
  4. 4 Summary and outlook
  5. A Jackknife method for the integrated autocorrelation times

##  1 Introduction

The numerical sign problem has prevented us from the first-principles analysis
of various important systems, such as quantum chromodynamics (QCD) at finite
density [1], quantum Monte Carlo calculations of quantum statistical systems
[2], and the real-time dynamics of quantum fields.

Among various approaches to the sign problem, some utilize the
complexification of dynamical variables. For example, in the complex Langevin
method [3, 4, 5, 6], one considers the Langevin equation in the complexified
configuration space. In the path optimization method [7, 8, 9, 10], with the
aid of machine learning one looks for an optimized integration surface for
which the average phase factor is maximal. In the Lefschetz thimble method
[11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23], one deforms the
integration surface according to the antiholomorphic gradient flow. The
deformed surface asymptotes to a union of Lefschetz thimbles, each of which
gives a constant value to the imaginary part of the action and thus is free
from the sign problem. Although there can appear the ergodicity problem due to
the existence of infinitely high potential barriers between thimbles, this
ergodicity problem can be diminished by tempering the system with the flow
time [19]. This tempered Lefschetz thimble method (TLTM) thus solves the sign
and ergodicity problems simultaneously. Moreover, the computational cost of
TLTM has recently been reduced significantly by developing the worldvolume
tempered Lefschetz thimble method (WV-TLTM) [23], which we are going to review
now.

Let ${\mathbb{R}}^{N}=\\{x\\}$ be the configuration space and $S(x)$ the
action (allowed to be complex-valued). Our main interest is to numerically
estimate the expectation values of observables $\mathcal{O}(x)$,

$\displaystyle\langle\mathcal{O}\rangle\equiv\frac{\int_{{\mathbb{R}}^{N}}dx\,e^{-S(x)}\,\mathcal{O}(x)}{\int_{{\mathbb{R}}^{N}}dx\,e^{-S(x)}}.$
(1.1)

Under the assumption that $e^{-S(z)}$ and $e^{-S(z)}\,\mathcal{O}(z)$ are
entire functions of $z\in{\mathbb{C}}^{N}$, Cauchy’s theorem allows us to
continuously deform the integration surface without changing the value of
integral. By expressing the deformation as a flow $z_{t}(x)$ $(t\geq 0)$ with
$z_{0}(x)=x$, the deformed integration surface at flow time $t$ can be written
as $\Sigma_{t}=\\{z_{t}(x)\,|\,x\in{\mathbb{R}}^{N}\\}$, and thus we have the
equality

$\displaystyle\langle\mathcal{O}\rangle=\frac{\int_{\Sigma_{t}}dz_{t}\,e^{-S(z_{t})}\,\mathcal{O}(z_{t})}{\int_{\Sigma_{t}}dz_{t}\,e^{-S(z_{t})}}.$
(1.2)

Since the numerator and the denominator are both independent of $t$, we can
integrate each of them over an arbitrary interval $[T_{0},T_{1}]$ with an
arbitrary function $W(t)$. This rewrites (1.2) to the integration over the
worldvolume $\mathcal{R}\equiv\bigcup_{T_{0}\leq t\leq T_{1}}\Sigma_{t}$:

$\displaystyle\langle\mathcal{O}\rangle=\frac{\int_{T_{0}}^{T_{1}}dt\,e^{-W(t)}\int_{\Sigma_{t}}dz_{t}\,e^{-S(z_{t})}\,\mathcal{O}(z_{t})}{\int_{T_{0}}^{T_{1}}dt\,e^{-W(t)}\int_{\Sigma_{t}}dz_{t}\,e^{-S(z_{t})}}=\frac{\int_{\mathcal{R}}Dz\,K(z)\,e^{-W(t(z))}\,e^{-S(z)}\,\mathcal{O}(z)}{\int_{\mathcal{R}}Dz\,K(z)\,e^{-W(t(z))}\,e^{-S(z)}},$
(1.3)

where $Dz$ is the induced volume element on $\mathcal{R}$, and $K(z)$ is the
Jacobian: $dt\,dz_{t}=K(z)\,Dz$. The weight factor $e^{-W(t)}$ is chosen so
that the probability for a configuration to appear at time $t$ is (almost)
independent of $t$. This setting is especially necessary when the whole range
of $t$ is relevant to simulations, as in the WV-TLTM. We further rewrite (1.3)
to the ratio of reweighted averages:

$\displaystyle\langle\mathcal{O}\rangle=\frac{\int_{\mathcal{R}}Dz\,e^{-V(z)}\,A(z)\,\mathcal{O}(z)}{\int_{\mathcal{R}}Dz\,e^{-V(z)}\,A(z)}=\frac{\langle
A(z)\,\mathcal{O}(z)\rangle_{\mathcal{R}}}{\langle
A(z)\rangle_{\mathcal{R}}}.$ (1.4)

Here, the reweighted average

$\displaystyle\langle\mathcal{O}(z)\rangle_{\mathcal{R}}\equiv\frac{\int_{\mathcal{R}}Dz\,e^{-V(z)}\,\mathcal{O}(z)}{\int_{\mathcal{R}}Dz\,e^{-V(z)}}$
(1.5)

is defined for the weight $e^{-V(z)}\equiv e^{-{\rm Re}\,S(z)-W(t(z))}$, and
$A(z)$ is the reweighting factor

$\displaystyle
A(z)\equiv\frac{e^{-S(z)-W(t(z))}dt\,dz_{t}}{e^{-V(z)}Dz}=e^{-i\,{\rm
Im}\,S(z)}K(z),$ (1.6)

whose explicit form can be found in Ref. [23].

The reweighted average (1.5) is estimated by the average over a sample
generated by the hybrid Monte Carlo (HMC) algorithm with the potential $V(z)$.
We will generally call a HMC algorithm on an accumulation of integration
surfaces the worldvolume hybrid Monte Carlo algorithm (the WV-HMC algorithm),
which includes the WV-TLTM.

In the case of the WV-TLTM, the interval $[T_{0},T_{1}]$ should include the
region with large $t$ so as to solve the sign problem, and at the same time
include the region with small $t$ so as to reduce the ergodicity problem.
However, the small-$t$ region may be contaminated by the sign problem, and the
large-$t$ region by the ergodicity problem. Thus, in order to reduce the
contributions from these potentially contaminated regions, it was proposed in
Ref. [23] to estimate the observables from a subsample in a subregion with
$[\tilde{T}_{0},\tilde{T}_{1}]$ ($\subset[T_{0},T_{1}]$) (see Fig. 1), where
$\tilde{T}_{0}$ and $\tilde{T}_{1}$ are chosen such that the estimated values
vary only slightly against the changes of $\tilde{T}_{0}$ and
$\tilde{T}_{1}$.111$\tilde{T}_{0,1}$ were written as $\hat{T}_{0,1}$ in Ref.
[23].

Figure 1:  The worldvolume $\mathcal{R}$ and a subregion
$\tilde{\mathcal{R}}$.

In Ref. [23] the standard error analysis was employed as if the subsample
itself is generated as a Markov chain with ergodicity and detailed balance.
However, this is not obvious and needs a justification. In this paper, we
prove that, if consecutive configurations in the worldvolume are generated as
a Markov chain with ergodicity and detailed balance, then the subset
consisting of the configurations belonging to a subregion can also be regarded
as a Markov chain with ergodicity and detailed balance intact. We particularly
investigate the integrated autocorrelation times for the Markov chains in
various subregions, and find that there is a linear relation between the
probability to be in a subregion and the integrated autocorrelation time for
the corresponding subsample. We numerically confirm this scaling law for a
chiral random matrix model (the Stephanov model [24, 25]).

This paper is organized as follows. In section 2, we first prove that the
subset consisting of the configurations in a subregion is a Markov chain. We
then argue that there should be a linear relation between the probability to
be in a subregion and the integrated autocorrelation time for the
corresponding subsample. Section 3 demonstrates this scaling by explicit
numerical calculations for the Stephanov model. Section 4 is devoted to
conclusion and outlook. In appendix A we explain our Jackknife method for
estimating the statistical errors of the integrated autocorrelation times.

##  2 Stochastic process in a subregion and the integrated autocorrelation
time

Let $\mathcal{R}=\\{z\\}$ be the full configuration space (the worldvolume in
the case of the WV-TLTM). Suppose that we are given a Markov chain
$\\{z_{k}\\}$ $(k=0,1,2,\ldots)$ in $\mathcal{R}$ with the transition
probability $P(z_{k+1}|z_{k})$ that satisfies ergodicity as well as the
detailed balance condition with respect to the unique equilibrium distribution
$p_{\rm eq}(z)$:

$\displaystyle P(z^{\prime}|z)\,p_{\rm eq}(z)=P(z|z^{\prime})\,p_{\rm
eq}(z^{\prime}).$ (2.1)

Only in this section, we write $\int_{\mathcal{R}}dz\,p_{\rm
eq}(z)\,\mathcal{O}(z)$ simply by $\langle\mathcal{O}\rangle$ (instead of
$\langle\mathcal{O}\rangle_{\mathcal{R}}$).

###  2.1 Stochastic process in a subregion

We now look at a subregion $\tilde{\mathcal{R}}$ in $\mathcal{R}$, whose
complement we denote by
$\tilde{\mathcal{R}}^{c}\equiv\mathcal{R}\backslash\tilde{\mathcal{R}}$.222In
the WV-TLTM, $\mathcal{R}=\bigcup_{t=T_{0}}^{T_{1}}\Sigma_{t}$ and
$\tilde{\mathcal{R}}=\bigcup_{t=\tilde{T}_{0}}^{\tilde{T}_{1}}\Sigma_{t}$ with
$T_{0}\leq\tilde{T}_{0}<\tilde{T}_{1}\leq T_{1}$.  Then, from the Markov chain
$\\{z_{k}\\}$ $(k=0,1,2,\ldots)$ in $\mathcal{R}$, we can extract a
subsequence $\\{\tilde{z}_{\ell}\\}$ $(\ell=0,1,2,\ldots)$ that consists of
the configurations belonging to $\tilde{\mathcal{R}}$. We first notice that
this sequence is a Markov chain with the following transition probability from
$\tilde{z}\in\tilde{\mathcal{R}}$ to
$\tilde{z}^{\prime}\in\tilde{\mathcal{R}}$:

$\displaystyle\tilde{P}(\tilde{z}^{\prime}|\tilde{z})=P(\tilde{z}^{\prime}|\tilde{z})$
$\displaystyle+\int_{\tilde{\mathcal{R}}^{c}}dw\,P(\tilde{z}^{\prime}|w)P(w|\tilde{z})$
$\displaystyle+\int_{\tilde{\mathcal{R}}^{c}}dw_{2}\,dw_{1}\,P(\tilde{z}^{\prime}|w_{2})P(w_{2}|w_{1})P(w_{1}|\tilde{z})$
$\displaystyle+\cdots.$ (2.2)

Since $P(z^{\prime}|z)$ is ergodic by assumption, and thus since a
configuration which has left $\tilde{\mathcal{R}}$ will eventually reenter
$\tilde{\mathcal{R}}$ at a finite number of steps,
$\tilde{P}(\tilde{z}^{\prime}|\tilde{z})$ is ergodic and satisfies the
probability conservation:

$\displaystyle\int_{\tilde{\mathcal{R}}}d\tilde{z}^{\prime}\,\tilde{P}(\tilde{z}^{\prime}|\tilde{z})=1\quad(\forall\tilde{z}\in\tilde{\mathcal{R}}).$
(2.3)

Furthermore, using the expression (2.2), one can easily show that
$\tilde{P}(\tilde{z}^{\prime}|\tilde{z})$ satisfies the following equality:

$\displaystyle\tilde{P}(\tilde{z}^{\prime}|\tilde{z})\,p_{\rm
eq}(\tilde{z})=\tilde{P}(\tilde{z}|\tilde{z}^{\prime})\,p_{\rm
eq}(\tilde{z}^{\prime})\quad(\tilde{z}^{\prime},\tilde{z}\in\tilde{\mathcal{R}}),$
(2.4)

from which we find that the unique equilibrium distribution $\tilde{p}_{\rm
eq}(\tilde{z})$ for $\tilde{P}(\tilde{z}^{\prime}|\tilde{z})$ is given by

$\displaystyle\tilde{p}_{\rm eq}(\tilde{z})=\frac{p_{\rm
eq}(\tilde{z})}{\int_{\tilde{\mathcal{R}}}d\tilde{z}^{\prime}\,p_{\rm
eq}(\tilde{z}^{\prime})}.$ (2.5)

We thus have proved that the estimation of the expectation value
$\int_{\tilde{\mathcal{R}}}d\tilde{z}\,\tilde{p}_{\rm
eq}(\tilde{z})\,\mathcal{O}(\tilde{z})$ with the subsample from a subregion
$\tilde{\mathcal{R}}$ can be statistically analyzed as if it is a Markov
chain.

###  2.2 Integrated autocorrelation time for a subchain

Let $\\{z_{k}\\}$ $(k=0,1,2,\ldots)$ again be a Markov chain in $\mathcal{R}$
with the transition probability $P(z^{\prime}|z)$. Denoting
$\mathcal{O}(z_{k})$ by $\mathcal{O}_{k}$, we define the integrated
autocorrelation time of $\mathcal{O}$ by333 This normalization gives the
effective sample size as $N^{\rm eff}_{\rm conf}=N_{\rm conf}/\tau_{\rm int}$
(see, e.g., Ref. [26]).

$\displaystyle\tau_{\rm int}(\mathcal{O})\equiv
1+2\sum_{k=1}^{\infty}\frac{C_{k}(\mathcal{O})}{C_{0}(\mathcal{O})},$ (2.6)

where
$C_{k}(\mathcal{O})\equiv\langle\mathcal{O}_{0}\mathcal{O}_{k}\rangle_{c}\equiv\langle(\mathcal{O}_{0}-\langle\mathcal{O}\rangle)(\mathcal{O}_{k}-\langle\mathcal{O}\rangle)\rangle$
is the autocorrelation. Similarly, we define the integrated autocorrelation
time for the subchain $\\{\tilde{z}_{\ell}\\}$ $(\ell=0,1,2,\ldots)$ in
$\tilde{\mathcal{R}}$, and denote it by $\tilde{\tau}_{\rm int}(\mathcal{O})$.
Note that $\tilde{\tau}_{\rm int}(\mathcal{O})$ is generically smaller than
$\tau_{\rm int}(\mathcal{O})$, because one-step transition with $\tilde{P}$
can correspond to transitions of multiple steps with $P$.

The ratio $\tilde{\tau}_{\rm int}(\mathcal{O})/\tau_{\rm int}(\mathcal{O})$
can be evaluated as follows, when both the numerator and the denominator are
not too small. We first write by $\epsilon$ and $\tilde{\epsilon}$,
respectively, the average Monte Carlo times evolved in one-step transition
with $P$ and $\tilde{P}$.444 This $\epsilon$ corresponds to the Langevin time
for Langevin algorithms and to the molecular dynamics time multiplied by the
average acceptance rate for HMC algorithms.  We then note that
$\tilde{\epsilon}$ can be written with $\epsilon$ by using the probability $p$
for a configuration in $\tilde{\mathcal{R}}$ to stay in $\tilde{\mathcal{R}}$
at the next step and the probability $q$ for a configuration in
$\tilde{\mathcal{R}}^{c}$ to stay in $\tilde{\mathcal{R}}^{c}$ as well:

$\displaystyle\tilde{\epsilon}$
$\displaystyle=p\,\epsilon+(1-p)(1-q)\,2\epsilon+(1-p)q(1-q)\,3\epsilon+\cdots$
$\displaystyle=\frac{2-p-q}{1-q}\,\epsilon$ $\displaystyle\Bigl{(}$
$\displaystyle=\frac{(1-p)+(1-q)}{1-q}\,\epsilon\geq\epsilon\Bigr{)}.$ (2.7)

Furthermore, since autocorrelations should be the same at large Monte Carlo
time scales, we have the equality

$\displaystyle\tau_{\rm int}(\mathcal{O})\,\epsilon=\tilde{\tau}_{\rm
int}(\mathcal{O})\,\tilde{\epsilon}.$ (2.8)

Note that this renormalization-group-like argument holds only when $\tau_{\rm
int}(\mathcal{O})$ and $\tilde{\tau}_{\rm int}(\mathcal{O})$ are not too
small. Combining Eq. (2.7) and Eq. (2.8), we obtain the desired result:

$\displaystyle\frac{\tilde{\tau}_{\rm int}(\mathcal{O})}{\tau_{\rm
int}(\mathcal{O})}=\frac{\epsilon}{\tilde{\epsilon}}=\frac{1-q}{2-p-q}\leq 1.$
(2.9)

###  2.3 Scaling law for the integrated autocorrelation times

We now apply the preceding arguments to the case where the configuration space
is the worldvolume of the WV-TLTM. We argue that there must be a linear
relation between the probability to be in a subregion $\tilde{\mathcal{R}}$
and the integrated autocorrelation time $\tilde{\tau}_{\rm int}$ for the
corresponding subchain. This claim will be confirmed numerically in the next
section.

To simplify discussions, we assume that the integrated autocorrelation time
for the flow time $t$ [i.e., $\tau_{\rm int}(\mathcal{O})$ with
$\mathcal{O}(z)=t(z)$] is sufficiently small, $\tau_{\rm int}(t)\simeq 1$.
This can be easily realized, if necessary, by removing consecutive
configurations from the chain at appropriate intervals. The smallness of
$\tau_{\rm int}(t)$ means that the probability $p$ simply expresses the
probability for a configuration to be in $\tilde{\mathcal{R}}$. We now recall
that the distribution of $t$ is uniform in equilibrium for the WV-TLTM [see
the discussion below Eq. (1.3)]. Thus, $p$ is given by

$\displaystyle p=\frac{\tilde{T}_{1}-\tilde{T}_{0}}{T_{1}-T_{0}}.$ (2.10)

A similar statement holds for the probability $q$ which now expresses the
probability for a configuration to be in $\tilde{\mathcal{R}}^{c}$, and thus
we obtain the relation $p+q=1$. Then, Eq. (2.7) leads to the relation
$\tilde{\epsilon}=\epsilon/p$, and thus, combined with Eq. (2.10), we obtain
the following scaling law for the integrated autocorrelation times:

$\displaystyle
p=\frac{\tilde{T}_{1}-\tilde{T}_{0}}{T_{1}-T_{0}}=\frac{\tilde{\tau}_{\rm
int}(\mathcal{O})}{\tau_{\rm int}(\mathcal{O})}.$ (2.11)

We now consider the numerical estimation of $\langle\mathcal{O}\rangle$ using
a subsample belonging to the interval $[\tilde{T}_{0},\tilde{T}_{1}]$. Since
Cauchy’s theorem ensures $\langle\mathcal{O}\rangle$ to be independent of the
choice of $[\tilde{T}_{0},\tilde{T}_{1}]$, the estimate should not vary
largely against the changes around an appropriately chosen pair
$[\tilde{T}_{0},\tilde{T}_{1}]$. Furthermore, the statistical error
$\delta\langle\mathcal{O}\rangle$ also hardly depends on the choice of
$[\tilde{T}_{0},\tilde{T}_{1}]$. To see this, let us write the number of
configurations in the interval $[\tilde{T}_{0},\tilde{T}_{1}]$ as $N_{\rm
conf}(\tilde{T}_{0},\tilde{T}_{1})$. Since the distribution of $t$ is almost
uniform, the ratio $N_{\rm conf}(\tilde{T}_{0},\tilde{T}_{1})/N_{\rm
conf}(T_{0},T_{1})$ almost equals $p$, and thus we obtain the relation

$\displaystyle\frac{N_{\rm conf}(\tilde{T}_{0},\tilde{T}_{1})}{N_{\rm
conf}(T_{0},T_{1})}\simeq\frac{\tilde{\tau}_{\rm int}(\mathcal{O})}{\tau_{\rm
int}(\mathcal{O})}.$ (2.12)

This in turn means that the effective number of configurations in a subsample
does not depend on the choice of the interval $[\tilde{T}_{0},\tilde{T}_{1}]$,
because

$\displaystyle N_{\rm conf}^{\rm
eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})\equiv\frac{N_{\rm
conf}(\tilde{T}_{0},\tilde{T}_{1})}{\tilde{\tau}_{\rm
int}(\mathcal{O})}\simeq\frac{N_{\rm conf}(T_{0},T_{1})}{\tau_{\rm
int}(\mathcal{O})}\quad(\forall\,\tilde{T}_{0},\,\tilde{T}_{1}).$ (2.13)

When $N_{\rm conf}$ is sufficiently large (as we always assume), the
statistical error is given by the formula
$\delta\langle\mathcal{O}\rangle=\sigma_{\mathcal{O}}/\sqrt{N_{\rm conf}^{\rm
eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})}$ with a constant
$\sigma_{\mathcal{O}}$ ($=\sqrt{\langle\mathcal{O}^{2}\rangle_{c}}$). Since
$N_{\rm conf}^{\rm eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})$ is almost
independent of the choice of $[\tilde{T}_{0},\tilde{T}_{1}]$, so is
$\delta\langle\mathcal{O}\rangle$.

##  3 Numerical confirmation of the scaling law

We numerically confirm the scaling law (2.11) for the WV-TLTM [23] applied to
a chiral random matrix model (the Stephanov model [24, 25]).

###  3.1 Setup

The partition function of the Stephanov model for $N_{f}$ quarks with equal
mass $m$ is given by

$\displaystyle Z_{n}^{N_{f}}=e^{n\mu^{2}}\int
d^{2}X\,e^{-S(X,X^{\dagger})}\equiv e^{n\mu^{2}}\int d^{2}X\,e^{-n\,{\rm
tr}\,X^{\dagger}X}\,{\det}^{N_{f}}(D+m),$ (3.1)

where $X=(X_{ij}=x_{ij}+i\,y_{ij})$ is an $n\times n$ complex matrix. The
$2n\times 2n$ matrix $D$ expresses the Dirac operator in the chiral
representation,

$\displaystyle D\equiv\left(\begin{array}[]{cc}0&i\,(X+C)\\\
i\,(X^{\dagger}+C)&0\end{array}\right),$ (3.4) $\displaystyle
i\,C\equiv\left(\begin{array}[]{cc}(\mu+i\tau)\,1_{n/2}&0\\\
0&(\mu-i\tau)\,1_{n/2}\end{array}\right),$ (3.7)

where $\mu$ and $\tau$ represent the chemical potential and the temperature,
respectively. The chiral condensate and the number density are defined,
respectively, by

$\displaystyle\langle\bar{\psi}\psi\rangle$
$\displaystyle\equiv\frac{1}{2n}\,\frac{\partial}{\partial m}\,\ln
Z_{n}^{N_{f}},\quad\langle\psi^{\dagger}\psi\rangle\equiv\frac{1}{2n}\,\frac{\partial}{\partial\mu}\,\ln
Z_{n}^{N_{f}}.$ (3.8)

We will set the parameters to $N_{f}=1$, $n=2$, $\mu=0.6$, $\tau=0$,
$m=0.004$.

We generate a sample $\\{z_{j}\\}$ $(j=1,\ldots,N_{\rm conf})$ with the HMC
algorithm using the potential $V(z)$, and estimate the reweighted average
$\langle\mathcal{O}(z)\rangle_{\mathcal{R}}$ [see Eq. (1.5)] by the sample
average

$\displaystyle\overline{\mathcal{O}}\equiv\frac{1}{N_{\rm
conf}}\sum_{j=1}^{N_{\rm conf}}\mathcal{O}(z_{j}).$ (3.9)

The estimator of the integrated autocorrelation time $\tau_{\rm
int}(\mathcal{O})$ is given by

$\displaystyle\overline{\tau}_{\rm int}(\mathcal{O};k_{\rm max})\equiv
1+2\sum_{k=1}^{k_{\rm
max}}\frac{\overline{C}_{k}(\mathcal{O})}{\overline{C}_{0}(\mathcal{O})}.$
(3.10)

Here, $\overline{C}_{k}(\mathcal{O})$ is an estimator of the autocorrelation
$C_{k}(\mathcal{O})=\langle\mathcal{O}_{0}\mathcal{O}_{k}\rangle_{\mathcal{R},\,c}$,
whose explicit form is given in appendix A. We have truncated the summation at
$k_{\rm max}$ to avoid summing up statistical fluctuations around zero (see,
e.g., Ref. [28]). The statistical error $\delta\overline{\tau}_{\rm
int}(\mathcal{O};k_{\rm max})$ is estimated by a Jackknife method that is
described in appendix A. Values of $k_{\rm max}$ and bin sizes used in the
estimations of $\overline{\tau}_{\rm int}(\mathcal{O};k_{\rm max})$ are
summarized in Table 1.

$\mathcal{O}$ | $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---|---|---
$t$ | $k_{\rm max}$ | 75 | 60 | 60 | 60 | 60 | 60 | 50 | 50 | 50 | 50
| bin size | 140 | 140 | 140 | 140 | 140 | 120 | 120 | 100 | 100 | 100
${\rm Re}\,(\bar{\psi}\psi A)$ | $k_{\rm max}$ | 900 | 250 | 250 | 250 | 250 | 300 | 300 | 300 | 400 | 500
| bin size | 1100 | 800 | 900 | 500 | 900 | 500 | 600 | 500 | 500 | 400
${\rm Re}\,(\psi^{\dagger}\psi A)$ | $k_{\rm max}$ | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500
| bin size | 1000 | 1300 | 1200 | 1300 | 1100 | 1100 | 1100 | 1100 | 1100 | 1300
${\rm Re}\,A$ | $k_{\rm max}$ | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500
| bin size | 1000 | 1000 | 900 | 700 | 700 | 600 | 600 | 500 | 500 | 500
$\mathcal{O}$ | $i$ | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19
$t$ | $k_{\rm max}$ | 50 | 50 | 50 | 50 | 50 | 40 | 30 | 25 | 25 | 15
| bin size | 100 | 100 | 100 | 100 | 120 | 80 | 100 | 120 | 100 | 100
${\rm Re}\,(\bar{\psi}\psi A)$ | $k_{\rm max}$ | 500 | 400 | 500 | 500 | 300 | 250 | 400 | 250 | 400 | 300
| bin size | 400 | 400 | 300 | 300 | 300 | 300 | 400 | 400 | 400 | 400
${\rm Re}\,(\psi^{\dagger}\psi A)$ | $k_{\rm max}$ | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 250 | 400
| bin size | 1100 | 1100 | 1300 | 1300 | 1100 | 1100 | 1100 | 1100 | 1100 | 1100
${\rm Re}\,A$ | $k_{\rm max}$ | 500 | 500 | 500 | 500 | 500 | 500 | 500 | 250 | 250 | 250
| bin size | 500 | 400 | 500 | 400 | 400 | 400 | 400 | 500 | 400 | 500
Table 1: Values of $k_{\rm max}$ and bin sizes.

In the numerical simulation with the HMC, we set $T_{0}=0$ and $T_{1}=0.02$.
The HMC updates are performed with the molecular dynamics time increment
$\Delta s=0.001$ and the step number $n_{\rm HMC}=100$ with the average
acceptance rate more than $0.99$. We employ 20 independent sets of
configurations, each set consisting of $3\times 10^{6}$ configurations in
$[T_{0},T_{1}]$. The observables are measured at every 6 iterations of the HMC
algorithm, so that we have 20 independent samples of size $N_{\rm
conf}=5\times 10^{5}$. We fix $\tilde{T}_{1}=T_{1}$ and vary $\tilde{T}_{0}$
as $\tilde{T}_{0}^{(i)}\equiv(i/20)T_{1}$ ($i=0,...,19$), for each of which an
independent set of configurations is used.

###  3.2 Results

We now demonstrate the scaling law (2.11) from explicit numerical
calculations. Recall that the argument for the scaling is based on the
smallness of the integrated autocorrelation time of $t$ ($\tau_{\rm
int}(t)\simeq 1$) and the uniformity of the distribution of $t$. Figure 2
shows that the condition $\tilde{\tau}_{\rm int}(t)\simeq 1$ is satisfied for
all $p=(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$.

Figure 2:  Results of $\tilde{\tau}_{\rm int}(t)$.

Figure 3 is the histogram of $p=(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$,
which is almost flat as required.

Figure 3:  Histogram of $t$. The number of configurations in a bin is the
average over 20 independent samples, each consisting of $N_{\rm conf}=5\times
10^{5}$. The function $W(t)$ is tuned such that the histogram becomes almost
flat for these eight bins.

This is realized by tuning the functional form of $W(t)$ as in Ref. [23].

Figure 4 exhibits the scaling law (2.11) for three operators,
$\mathcal{O}={\rm Re}\,(\bar{\psi}\psi A)$, $\mathcal{O}={\rm
Re}\,(\psi^{\dagger}\psi A)$, and $\mathcal{O}={\rm Re}\,A$. We see that the
scaling is satisfied for $\mathcal{O}={\rm Re}\,(\bar{\psi}\psi A)$ in the
region $(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})\geq 0.35$, and for
$\mathcal{O}={\rm Re}\,(\psi^{\dagger}\psi A)$ and $\mathcal{O}={\rm Re}\,A$
in the region $(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})\geq 0.20$.
Deviations from the scaling at small
$(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$ should be due to
$\tilde{\tau}_{\rm int}(\mathcal{O})=O(1)$ [see the comment below Eq. (2.8)].

Figure 4:  $\tilde{\tau}_{\rm int}(\mathcal{O})$ against
$(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$. Top left: $\mathcal{O}={\rm
Re}\,(\bar{\psi}\psi A)$. Top right: $\mathcal{O}={\rm
Re}\,(\psi^{\dagger}\psi A)$. Bottom: $\mathcal{O}={\rm Re}\,A$. Data points
in the shaded region give $\tilde{\tau}_{\rm int}(\mathcal{O})=O(1)$ and thus
are excluded from the linear fitting [see the comment below Eq. (2.8)].

We perform the $\chi^{2}$-fit to these data points with

$\displaystyle\chi^{2}(\mathcal{O})\equiv\sum_{i=0}^{i_{\rm
max}}\frac{[\tilde{\tau}^{(i)}_{\rm
int}(\mathcal{O})-c\,(\tilde{T}_{1}-\tilde{T}_{0}^{(i)})/(T_{1}-T_{0})]^{2}}{[\delta\tilde{\tau}_{\rm
int}^{(i)}(\mathcal{O})]^{2}},$ (3.11)

where $\tilde{\tau}_{\rm int}^{(i)}(\mathcal{O})$ is the integrated
autocorrelation time for the subsample with the interval
$[\tilde{T}_{0},\tilde{T}_{1}]=[\tilde{T}_{0}^{(i)},T_{1}]$. The fit results
are the following: For $\mathcal{O}={\rm Re}\,(\bar{\psi}\psi A)$ with $i_{\rm
max}=13$, $c=3.1$ and $\chi^{2}/{\rm DOF}=1.0$. For $\mathcal{O}={\rm
Re}\,(\psi^{\dagger}\psi A)$ with $i_{\rm max}=16$, $c=16.4$ and
$\chi^{2}/{\rm DOF}=1.2$. For $\mathcal{O}={\rm Re}\,A$ with $i_{\rm max}=16$,
$c=15.5$ and $\chi^{2}/{\rm DOF}=1.1$. These data support the scaling law
(2.11).

We plot in Fig. 5 the values of $N_{\rm conf}^{\rm
eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})$ for various operators
$\mathcal{O}$ [see Eq. (2.13)]. The statistical errors are estimated with the
Jackknife method. We observe, as expected, that $N_{\rm conf}^{\rm
eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})$ takes almost the same values
for each $\mathcal{O}$ in the range where we observe the scaling law.

Figure 5:  The values of $N_{\rm conf}^{\rm
eff}(\mathcal{O};\tilde{T}_{0},\tilde{T}_{1})$ [see Eq. (2.13)]. Top left:
$\mathcal{O}={\rm Re}\,(\bar{\psi}\psi A)$. Top right: $\mathcal{O}={\rm
Re}\,(\psi^{\dagger}\psi A)$. Bottom: $\mathcal{O}={\rm Re}\,A$. The regions
where the scaling law (2.11) is broken are shaded.

Finally, in Fig. 6 we plot the expectation values of the chiral condensate and
the number density, together with their statistical errors.

Figure 6:  The estimate (top) and its statistical error (bottom) of the
expectation value $\langle\mathcal{O}\rangle$ against
$(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$. Left:
$\mathcal{O}=\bar{\psi}\psi$ (the chiral condensate). Right:
$\mathcal{O}=\psi^{\dagger}\psi$ (the number density). The regions where the
scaling law (2.11) is broken in the numerator and/or the denominator are
shaded for the statistical errors.

The statistical errors are estimated with the Jackknife method with the bin
sizes fixed to 500. We see that both the means and the statistical errors take
almost constant values in a region where
$(\tilde{T}_{1}-\tilde{T}_{0})/(T_{1}-T_{0})$ is not small. The deviation of
the means should be attributed to the complex geometry at large flow times,
which requires larger statistics and better control of systematic errors (such
as those from numerical integrations of the flow equation and of the Hamilton
equations accompanied by the projection from ${\mathbb{C}}^{N}$ to
$\mathcal{R}$). The deviation of the statistical errors is due to the
violation of the scaling law (2.11) for either of the numerator or the
denominator (or both) in the ratio of the reweighted averages.

##  4 Summary and outlook

In this paper, we have established the statistical analysis method for the WV-
HMC algorithm, whose major use is intended for the WV-TLTM [23]. We proved
that, if consecutive configurations in the worldvolume are generated as a
Markov chain with ergodicity and detailed balance, then the subset consisting
of the configurations belonging to a subregion can also be regarded as a
Markov chain with ergodicity and detailed balance intact. We particularly
investigated the integrated autocorrelation times for the Markov chains in
various subregions, and found that there is a linear relation between the
probability to be in a subregion and the integrated autocorrelation time for
the corresponding subsample. We numerically confirmed this scaling law for the
Stephanov model.

Now with this statistical analysis method at hand, we can safely apply the WV-
TLTM to large-scale simulations of the systems that have serious sign
problems, such as finite density QCD, strongly correlated electron systems,
frustrated spin systems, and the real-time dynamics of quantum fields. A study
along this line is now in progress and will be reported elsewhere.

## Acknowledgments

The authors thank Issaku Kanamori, Yoshio Kikukawa and Jun Nishimura for
useful discussions. This work was partially supported by JSPS KAKENHI (Grant
Numbers 18J22698, 20H01900, 21K03553) and by SPIRITS (Supporting Program for
Interaction-based Initiative Team Studies) of Kyoto University (PI: M.F.).
N.M. is supported by the Special Postdoctoral Researchers Program of RIKEN.

##  Appendix A Jackknife method for the integrated autocorrelation times

In this appendix, we give a Jackknife method to estimate the integrated
autocorrelation times $\tau_{\rm int}(\mathcal{O})$.

Let $\\{x_{j}\\}$ $(j=1,\ldots,N_{\rm conf})$ be a set of consecutive
configurations generated as a Markov chain. We estimate the expectation value
$\langle\mathcal{O}\rangle$ by the sample average

$\displaystyle\overline{\mathcal{O}}\equiv\frac{1}{N_{\rm
conf}}\sum_{j=1}^{N_{\rm conf}}\mathcal{O}(x_{j}).$ (A.1)

The estimator of the integrated autocorrelation time $\tau_{\rm
int}(\mathcal{O})\equiv
1+2\sum_{k=1}^{\infty}C_{k}(\mathcal{O})/C_{0}(\mathcal{O})$ is given by

$\displaystyle\bar{\tau}_{\rm int}(\mathcal{O};k_{\rm max})\equiv
1+2\sum_{k=1}^{k_{\rm
max}}\frac{\overline{C}_{k}(\mathcal{O})}{\overline{C}_{0}(\mathcal{O})},$
(A.2)

where $\overline{C}_{k}(\mathcal{O})$ is the estimator of the autocorrelation
$C_{k}(\mathcal{O})\equiv\langle\mathcal{O}_{0}\,\mathcal{O}_{k}\rangle_{c}$.
The summation is truncated at $k_{\rm max}$ to avoid summing up statistical
fluctuations around zero (see, e.g., Refs. [27, 28]). The value of $k_{\rm
max}$ should not be set very large compared to $\tau_{\rm int}(\mathcal{O})$,
otherwise contributions from statistical fluctuations around zero may dominate
the error. There has been known an explicit formula for the statistical error
$\delta\tau_{\rm int}(\mathcal{O})$ when $1\ll k_{\rm max}\ll N_{\rm conf}$
(more precisely, as $N_{\rm conf}\to\infty$, $k_{\rm max}\to\infty$ and
$k_{\rm max}/N_{\rm conf}\to 0$) [26] (see also Refs. [27, 28]),555 This is
given by $\delta\tau_{\rm int}(\mathcal{O})=\tau_{\rm
int}(\mathcal{O})\,\sqrt{2(2k_{\rm max}+1)/N_{\rm conf}}$ for the estimator of
the autocorrelation, $\overline{C}_{k}(\mathcal{O})=(N_{\rm
conf}-k)^{-1}\,\sum_{n=1}^{N_{\rm
conf}-k}(\mathcal{O}(x_{n})-\overline{\mathcal{O}})(\mathcal{O}(x_{n+k})-\overline{\mathcal{O}})$
[26].  but this may not be applicable to the case when $\tau_{\rm
int}(\mathcal{O})=O(1)$, for which the condition $k_{\rm max}\gg 1$ cannot be
met. Therefore, in this paper we adopt the Jackknife method for the estimation
of $\delta\tau_{\rm int}(\mathcal{O})$.

In order to apply a resampling method of Jackknife, we introduce a sample of
multidimensional observables
$\\{X_{r}=(X_{r,k}^{\mathcal{O}\mathcal{O}},X_{r}^{\mathcal{O}})\\}$
($r=1,\ldots,N_{\rm conf}-k_{\rm max}$) with

$\displaystyle X_{r,k}^{\mathcal{O}\mathcal{O}}\equiv\frac{1}{k_{\rm
max}-k+1}\sum_{i=0}^{k_{\rm
max}-k}\mathcal{O}(x_{r+i})\mathcal{O}(x_{r+k+i})\quad(k=0,\ldots,k_{\rm
max}),$ (A.3) $\displaystyle X_{r}^{\mathcal{O}}\equiv\frac{1}{k_{\rm
max}+1}\sum_{i=0}^{k_{\rm max}}\mathcal{O}(x_{r+i}).$ (A.4)

Since $\langle
X_{r,k}^{\mathcal{O}\mathcal{O}}\rangle=\langle\mathcal{O}(x_{0})\,\mathcal{O}(x_{k})\rangle$
and $\langle X_{r}^{\mathcal{O}}\rangle=\langle\mathcal{O}(x)\rangle$, the
autocorrelations
$C_{k}(\mathcal{O})=\langle\mathcal{O}(x_{0})\mathcal{O}(x_{k})\rangle_{c}$
can be estimated by

$\displaystyle\overline{C}^{(X)}_{k}(\mathcal{O})\equiv\frac{1}{N_{\rm
conf}-k_{\rm max}}\sum_{r=1}^{N_{\rm conf}-k_{\rm
max}}X_{r,k}^{\mathcal{O}\mathcal{O}}-\bigg{(}\frac{1}{N_{\rm conf}-k_{\rm
max}}\sum_{r=1}^{N_{\rm conf}-k_{\rm max}}X_{r}^{\mathcal{O}}\bigg{)}^{2}.$
(A.5)

Since $\tau_{\rm int}(\mathcal{O})$ is a function of $C_{k}(\mathcal{O})$, it
can be estimated with $\overline{C}^{(X)}_{k}(\mathcal{O})$ as

$\displaystyle\tau_{\rm int}(\mathcal{O})\approx 1+2\sum_{k=1}^{k_{\rm
max}}\frac{\overline{C}^{(X)}_{k}(\mathcal{O})}{\overline{C}^{(X)}_{0}(\mathcal{O})}.$
(A.6)

The statistical error $\delta\tau_{\rm int}(\mathcal{O})$ can then be
estimated with the standard Jackknife method.

## References

  * [1] G. Aarts, Introductory lectures on lattice QCD at nonzero baryon number. J. Phys. Conf. Ser. 706, no. 2, 022004 (2016) [arXiv:1512.05145 [hep-lat]].
  * [2] L. Pollet, Recent developments in Quantum Monte-Carlo simulations with applications for cold gases, Rep. Prog. Phys. 75, 094501 (2012) [arXiv:1206.0781 [cond-mat.quant-gas]].
  * [3] G. Parisi, “On complex probabilities,” Phys. Lett. B 131, 393 (1983).
  * [4] J.R. Klauder, “Coherent State Langevin Equations for Canonical Quantum Systems With Applications to the Quantized Hall Effect,” Phys. Rev. A 29, 2036 (1984).
  * [5] G. Aarts, F. A. James, E. Seiler and I. O. Stamatescu, “Adaptive stepsize and instabilities in complex Langevin dynamics,” Phys. Lett. B 687, 154-159 (2010) [arXiv:0912.0617 [hep-lat]].
  * [6] J. Nishimura and S. Shimasaki, “New Insights into the Problem with a Singular Drift Term in the Complex Langevin Method,” Phys. Rev. D 92, no.1, 011501 (2015) [arXiv:1504.08359 [hep-lat]].
  * [7] Y. Mori, K. Kashiwa and A. Ohnishi, “Toward solving the sign problem with path optimization method,” Phys. Rev. D 96, no.11, 111501 (2017) [arXiv:1705.05605 [hep-lat]].
  * [8] Y. Mori, K. Kashiwa and A. Ohnishi, “Application of a neural network to the sign problem via the path optimization method,” PTEP 2018, no.2, 023B04 (2018) [arXiv:1709.03208 [hep-lat]].
  * [9] A. Alexandru, P. F. Bedaque, H. Lamm and S. Lawrence, “Finite-Density Monte Carlo Calculations on Sign-Optimized Manifolds,” Phys. Rev. D 97, no.9, 094510 (2018) [arXiv:1804.00697 [hep-lat]].
  * [10] F. Bursa and M. Kroyter, “A simple approach towards the sign problem using path optimisation,” JHEP 12, 054 (2018) [arXiv:1805.04941 [hep-lat]].
  * [11] E. Witten, “Analytic continuation of Chern-Simons theory,” AMS/IP Stud. Adv. Math. 50, 347-446 (2011) [arXiv:1001.2933 [hep-th]].
  * [12] M. Cristoforetti, F. Di Renzo and L. Scorzato, “New approach to the sign problem in quantum field theories: High density QCD on a Lefschetz thimble,” Phys. Rev. D 86, 074506 (2012) [arXiv:1205.3996 [hep-lat]].
  * [13] M. Cristoforetti, F. Di Renzo, A. Mukherjee and L. Scorzato, “Monte Carlo simulations on the Lefschetz thimble: Taming the sign problem,” Phys. Rev. D 88, no. 5, 051501(R) (2013) [arXiv:1303.7204 [hep-lat]].
  * [14] H. Fujii, D. Honda, M. Kato, Y. Kikukawa, S. Komatsu and T. Sano, “Hybrid Monte Carlo on Lefschetz thimbles \- A study of the residual sign problem,” JHEP 1310, 147 (2013) [arXiv:1309.4371 [hep-lat]].
  * [15] H. Fujii, S. Kamata and Y. Kikukawa, “Lefschetz thimble structure in one-dimensional lattice Thirring model at finite density,” JHEP 11, 078 (2015) [erratum: JHEP 02, 036 (2016)] [arXiv:1509.08176 [hep-lat]].
  * [16] H. Fujii, S. Kamata and Y. Kikukawa, “Monte Carlo study of Lefschetz thimble structure in one-dimensional Thirring model at finite density,” JHEP 12, 125 (2015) [erratum: JHEP 09, 172 (2016)] [arXiv:1509.09141 [hep-lat]].
  * [17] A. Alexandru, G. Başar and P. Bedaque, “Monte Carlo algorithm for simulating fermions on Lefschetz thimbles,” Phys. Rev. D 93, no. 1, 014504 (2016) [arXiv:1510.03258 [hep-lat]].
  * [18] A. Alexandru, G. Başar, P. F. Bedaque, G. W. Ridgway and N. C. Warrington, “Sign problem and Monte Carlo calculations beyond Lefschetz thimbles,” JHEP 1605, 053 (2016) [arXiv:1512.08764 [hep-lat]].
  * [19] M. Fukuma and N. Umeda, “Parallel tempering algorithm for integration over Lefschetz thimbles,” PTEP 2017, no. 7, 073B01 (2017) [arXiv:1703.00861 [hep-lat]].
  * [20] A. Alexandru, G. Başar, P. F. Bedaque and N. C. Warrington, “Tempered transitions between thimbles,” Phys. Rev. D 96, no. 3, 034513 (2017) [arXiv:1703.02414 [hep-lat]].
  * [21] M. Fukuma, N. Matsumoto and N. Umeda, “Applying the tempered Lefschetz thimble method to the Hubbard model away from half filling,” Phys. Rev. D 100, no. 11, 114510 (2019) [arXiv:1906.04243 [cond-mat.str-el]].
  * [22] M. Fukuma, N. Matsumoto and N. Umeda, “Implementation of the HMC algorithm on the tempered Lefschetz thimble method,” [arXiv:1912.13303 [hep-lat]].
  * [23] M. Fukuma and N. Matsumoto, “Worldvolume approach to the tempered Lefschetz thimble method,” PTEP 2021, no.2, 023B08 (2021) [arXiv:2012.08468 [hep-lat]].
  * [24] M. A. Stephanov, “Random matrix model of QCD at finite density and the nature of the quenched limit,” Phys. Rev. Lett. 76, 4472 (1996) [hep-lat/9604003].
  * [25] M. A. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov and J. J. M. Verbaarschot, “On the phase diagram of QCD,” Phys. Rev. D 58, 096007 (1998) [hep-ph/9804290].
  * [26] N. Madras and A. D. Sokal, “The pivot algorithm: A highly efficient Monte Carlo method for self-avoiding walk,” J. Stat. Phys. 50, 109 (1988).
  * [27] U. Grenander and M. Rosenblatt, “Statistical Spectral Analysis of Time Series Arising from Stationary Stochastic Processes,” Ann. Math. Statist. 24 (4) 537 - 558 (1953).
  * [28] M. B. Priestley, “Spectral analysis and time series,” Academic Press, London (1981).