# Characterizing the emission region property of blazars

Junhui Fan Center for Astrophysics, Guangzhou University, Guangzhou 510006,
China  Key Laboratory for Astronomical Observation and Technology of
Guangzhou, Guangzhou, 510006, China  Astronomy Science and Technology Research
Laboratory of Department of Education of Guangdong Province,
Guangzhou, 510006, China  Hubing Xiao Shanghai Key Lab for Astrophysics,
Shanghai Normal University, Shanghai, 200234, China  Wenxin Yang  Center for
Astrophysics, Guangzhou University, Guangzhou 510006, China  Key Laboratory
for Astronomical Observation and Technology of Guangzhou, Guangzhou, 510006,
China  Astronomy Science and Technology Research Laboratory of Department of
Education of Guangdong Province,
Guangzhou, 510006, China  Lixia Zhang  Center for Astrophysics, Guangzhou
University, Guangzhou 510006, China  Key Laboratory for Astronomical
Observation and Technology of Guangzhou, Guangzhou, 510006, China  Astronomy
Science and Technology Research Laboratory of Department of Education of
Guangdong Province,
Guangzhou, 510006, China  Anton A. Strigachev  Institute of Astronomy and
National Astronomical Observatory, Bulgarian Academy of Sciences,
72 Tsarigradsko shosse Boulevard, 1784 Sofia, Bulgaria  Rumen S. Bachev
Institute of Astronomy and National Astronomical Observatory, Bulgarian
Academy of Sciences,
72 Tsarigradsko shosse Boulevard, 1784 Sofia, Bulgaria  Jianghe Yang
Department of Physics and Electronics Science, Hunan University of Arts and
Science, Changde 415000, China  Center for Astrophysics, Guangzhou University,
Guangzhou 510006, China

###### Abstract

The studies and constraints on the emission region are crucial to the blazar
radiation mechanism. Yet the previous works mainly focus on individual
sources. In this work, we make use of the largest and the latest spectral
energy distribution (SED) fitting results in the literature to statistically
study the blazar emission region property in the framework of leptonic one-
zone. Our results reveal (1) FSRQs show lower electron energy ($\gamma_{\rm
p}\lesssim 1.6\times 10^{3}$) than BL Lacs and tend to have a stronger
magnetic field ($B$) and smaller electron-to-magnetic energy ratio ($U_{\rm
e}/U_{\rm B}$) than BL Lacs; (2) we find the electro-magnetic equipartition
would rather happen in the jets of BL Lacs than happen in the jets of FSRQs;
(3) there are 682 blazars with a magnetic field weaker critical value of
generating the Kelvin-Helmholtz instability, thus one-third of the blazars in
our sample are able to produce this instability; (4) the distance ($d_{\rm
em}$) between the emission region and the central black hole (BH) is in the
scale of $\sim$0.1 pc, the location of the emission region may be evenly
distributed inside and outside the broad line region (BLR).

##  1 Introduction

Active galactic nuclei (AGNs), one of the most popular extragalactic objects
in astronomy, emit radiation in the overall electromagnetic spectrum. The
ultimate energy source of AGNs is believed to be the gravitational potential
of supermassive black holes (BHs), which are embedded in its center, rather
than the nuclear fusion of stars (Lynden-Bell, 1969). An accretion disk can be
formed surrounding the BHs by matter losing angular momentum before falling
onto the BHs (Rees, 1984; Cao & Spruit, 2013). AGNs are divided into radio-
loud and radio-quiet ones according to their relative radio emission intensity
compared to their optical emission intensity (Strittmatter et al., 1980;
Kellermann et al., 1989; Xiao et al., 2022c). This dichotomy is mainly caused
by the presence of a strong relativistic and collimated jet in the radio-loud
AGNs (Urry & Padovani, 1995). The jet is so powerful that dominates the entire
emission of radio-loud AGNs, but the mechanism of jet launching is still
controversial. Blandford et al. (1977) suggested that the jet is powered by
extracting the rotational energy of BH (B-Z process), while Blandford & Payne
(1982) suggests that the jet is powered by the rotational energy of accretion
disk (B-P process) (see also in Xiong & Zhang, 2014; Xiao et al., 2022b; Zhang
et al., 2022).

Blazars, as an extreme subclass of radio-loud AGNs, exhibit distinctive
observational properties. The rapid and large amplitude variability, high and
variable polarization, strong and variable $\gamma$-ray emissions, and
apparent superluminal motion of blazars have been observed and studied (Wills
et al., 1992; Urry & Padovani, 1995; Villata et al., 2006; Fan, 2002; Fan et
al., 2014, 2021; Gupta et al., 2016; Xiao et al., 2019, 2020, 2022d; Abdollahi
et al., 2020). These properties are believed as consequences of a Doppler
beaming effect due to a small viewing angle ($\theta$) between the jet axis
and line of sight (eg., Ghisellini et al., 1993; Fan et al., 2013; Pei et al.,
2016; Xiao et al., 2020). The Doppler beaming effect is usually pronounced as
the time dilation and intensity amplification through a Doppler factor
($\delta=[\Gamma(1-\beta{\rm cos}\theta)]^{-1}$, where $\Gamma$ is the bulk
Lorentz factor and $\beta$ is the jet speed in units of the speed of light,
$c$), which can be estimated indirectly (Ghisellini et al., 1993; Readhead,
1994; Fan, 2005; Fan et al., 2013; Liodakis et al., 2018). There are two
subclasses of blazars, namely flat spectrum radio quasars (FSRQ) and BL
Lacertae objects (BL Lacs). The former one characterizes an optical spectrum
with strong emission lines (rest-frame equivalent width, ${\rm EW>5\AA}$),
while the latter one demonstrates no or weak emission features (${\rm
EW<5\AA}$) (Urry & Padovani, 1995; Scarpa & Falomo, 1997). Broadband studies
illustrate a typical two-hump spectral energy distribution (SED), the lower
energy bump ranges from radio to X-ray band and peaks at the infrared (IR) to
X-ray band, which is believed to be the synchrotron emission of the
relativistic electrons in the jet, the higher energy bump starts from the
X-ray band to the $\gamma$-ray band and peaks at the X-ray to GeV $\gamma$-ray
band (Abdo et al., 2010; Fan et al., 2016; Paliya et al., 2021; Yang et al.,
2022, 2023). The radiation mechanism of the higher energy bump is still
controversial, the leptonic model suggests that the higher energy bump is
attributed to the inverse Compton (IC) process (Blandford & Koenigl, 1979;
Sikora et al., 1994; Sokolov & Marscher, 2005; H. E. S. S. Collaboration et
al., 2015; Xue et al., 2019; Tan et al., 2020; Wang et al., 2022a), while the
hadronic model interprets it through the proton synchrotron radiation and
secondary particle cascade (Mücke & Protheroe, 2001; Dimitrakoudis et al.,
2012; Diltz et al., 2015; Cerruti et al., 2015; IceCube Collaboration et al.,
2018; Xue et al., 2021; Wang et al., 2022b).

Broadband SED modelling is an efficient approach to investigate the jet
properties, eg., constraining the magnetic field ($B$) of the emission region,
the Doppler factor, the emission region size ($R$), and the electron energy
distribution (EED) (Massaro et al., 2004; Tramacere et al., 2011; Ghisellini
et al., 2014; Ghisellini & Tavecchio, 2015; Chen, 2017). However, there are
degeneracies between these critical parameters (Kubo et al., 1998; Ghisellini
et al., 2014). Additional methods are employed to further constrain these
parameters, for instance, the variability time scale is used to give an upper
limit of $R\leq c\Delta t\delta/(1+z)$ due to the causality; the Owens Valley
Radio Observatory (OVRO) made efforts of estimating $\delta$ through radio
lightcurve of blazar flare/outburst (Liodakis et al., 2018); the Monitoring of
Jets in AGN with VLBA Experiments (MOJAVE) program monitors the radio
brightness, polarization variation and apparent motion (Lister et al., 2015,
2019), and the apparent speed is an indicator of $\delta$ (Zhang & Fan, 2008;
Xiao et al., 2019, 2020); The polarization observation is critical to the
study of blazar magnetic field. The observed polarization shows a frequency-
dependence of rotation measure (RM) in 3C 273 and suggests the magnetic field
may be structured helically on a large scale (Wardle, 2018; Hovatta et al.,
2019). Recently, constraining magnetic field strength with simultaneously
observed optical and radio linear polarization and circular polarization was
suggested by Liodakis et al. (2022), who proposed a formula that the magnetic
field strength is proportional to the square of circular polarization and
inversely proportional to the linear polarization. Based on this, they
rejected that the high-energy emission models requiring high magnetic field
strength and a low positron fraction.

It is also accessible to determine magnetic field, Doppler factor, emission
region size, etc, via SED features. Kubo et al. (1998) considered a population
of relativistic electrons, which forms an EED of broken power-law with a
breakpoint at $\gamma_{\rm p}$, and assumed that the lower-energy component
SED peaks at a frequency corresponding to that radiated by the electrons with
$\gamma_{\rm p}$. Coupling the assumption with the synchrotron and IC
radiation mechanism, they calculated the strength of the magnetic field and
the electron Lorentz factor with the energy and intensity of the two SED
peaks. Similarly, Chen (2018) applied the same method with a log-parabola EED
and obtained formulas to calculate these parameters describing jet properties.
Benefiting from the observations of Fermi Large Area Telescope (LAT), a large
number of blazars have been discovered (Abdollahi et al., 2020), and the study
of blazar comes to its era of prosperity.

In this work, taking the advantage of the latest and largest blazar SED
fitting results of Yang et al. (2022) and Yang et al. (2023), we aim to
investigate the property of blazar emission region and put constraints on
those critical parameters. This paper is organized as follows: In section 2 we
present our sample; Our method, analysis and results will be presented in
Section 3; The discussions in Section 4; Our conclusions will be given in
Section 5.

##  2 Sample

To study the blazar emission region and constrain its relevant parameters, we
need to connect these parameters with observational quantities. We collected a
sample of Fermi blazars with available synchrotron peak frequency
($\log\nu_{\rm sy}$) and the corresponding luminosity ($\log L_{\rm sy}$) from
Yang et al. (2022) and the IC peak frequency ($\log\nu_{\rm IC}$) and
luminosity ($\log L_{\rm IC}$) from Yang et al. (2023). In total, we have 2708
sources, including 759 FSRQs, 1141 BL Lacs, and 808 blazar candidates of
uncertain type (BCUs). There are 1791 sources in our sample with available
redshift, including 750 FSRQs with an average value of $\langle z_{\rm
F}\rangle=1.201\pm 0.647$, 843 BL Lacs with average value $\langle z_{\rm
B}\rangle=0.528\pm 0.501$, and 198 BCUs with average value $\langle z_{\rm
U}\rangle=0.919\pm 0.738$. In this sample, there are 512 blazars are
associated with sources in Liodakis et al. (2018) and thus with available
Doppler factor ($\delta$), with an average value of $\langle\delta_{\rm
F}\rangle=17.39\pm 13.42$ for FSRQs, an average value of $\langle\delta_{\rm
B}\rangle=11.26\pm 10.29$ for BL Lacs, and an average value of
$\langle\delta_{\rm U}\rangle=12.50\pm 13.37$ for BCUs.

Table 1: The magnetic field strength and electron energy for 2708 Fermi blazars. 4FGL Name | Class_1 | Class_2 | z | $\delta$ | $\log\nu_{\rm sy}$ | $\log\nu_{\rm IC}$ | P1 | $\log B$ | $\log\gamma_{\rm p}$ | $\log U_{\rm e}$ | $B_{\rm c}$
---|---|---|---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
J0001.2-0747 | BLL | ISP |  |  | 14.1 | 23.14 | -0.12 | -2.25 | 4.46 | 1.82 | 28.82
J0001.2+4741 | BCU | ISP |  |  | 14.1 | 22.50 | -0.09 | $-1.56\sim 1.64$ | $4.14\sim 2.54$ | $0.68\sim-4.26$ |
J0001.5+2113 | FSRQ | LSP | 1.11 |  | 13.2 | 20.62 | -0.18 | 2.77 | 1.47 | -5.86 | 0.004
J0002.4-5156 | BCU | HSP |  |  | 15.7 | 23.99 | -0.09 | $0.15\sim 0.25$ | $4.08\sim 4.03$ | $-2.29\sim-2.44$ |
J0003.1-5248 | BCU | HSP |  |  | 15.9 | 24.31 | -0.07 | $0.23\sim 0.13$ | $4.14\sim 4.19$ | $-2.24\sim-2.09$ |
… | … | … | … | … | … | … | … | … | … | … | …

##  3 Method and Results

###  3.1 The basics of one-zone leptonic model

In the leptonic frame, the higher energy emission results from the inverse
Compton scattering of internal or external soft photons, namely the
synchrotron self-Compton (SSC) and external Compton (EC) process. The SSC
process is believed to dominate the radiation in BL Lacs jets, while the EC
process dominates over the SSC process for FSRQs jets. The synchrotron power
is mainly produced by those electrons with Lorentz factor of $\gamma_{\rm p}$
and contribute most to the synchrotron peak, and the synchrotron peak
frequency in the observer frame is given by

$\nu_{\rm sy}=3.7\times 10^{6}\,\gamma_{\rm
p}^{2}\,B\,\frac{\delta}{1+z}\,{\rm Hz},$ (1)

where $B$ in units of Gs (Tavecchio et al., 1998). Correspondingly, the
synchrotron peak luminosity in the observer frame is expressed as (Chen, 2018)

$L_{\rm sy}=4\pi\,\frac{16}{9}\pi R^{3}\,\frac{\sigma_{\rm T}c}{8\pi}\,U_{\rm
B}\,N_{\rm 0}\,\gamma_{\rm p}^{3}\,\delta^{4},$ (2)

where $\sigma_{\rm T}$ is the Thomson cross section, $U_{\rm B}$ is the
magnetic field energy density ($U_{\rm B}=B^{2}/8\pi$), $N_{\rm 0}$ is the
normalization parameter of the EED. Following Chen (2018), a three-parameter
log-parabolic function is employed to describe the EED in this work

$N(\gamma)=N_{\rm 0}\left(\frac{\gamma}{\gamma_{\rm
p}}\right)^{-3}10^{-2b\log(\gamma/\gamma_{\rm p})},$ (3)

$b$ is the curvature and has a relation with the synchrotron bump, $b\simeq
5P_{\rm 1}$ (see, eg., Massaro et al., 2006; Chen, 2014). This EED is only
phenomenologically assumed to follow the log-parabola shape of SED, without
taking into account the evolution due to injection and cooling effects.

Moreover, if the electrons in the Thomson regime, the peak frequency of the
SSC component is

$\nu_{\rm SSC}=\frac{4}{3}\,\gamma_{\rm p}^{2}\,\nu_{\rm sy}.$ (4)

In the case of the EC process, soft photons are fed externally, the peak
frequency is given by

$\nu_{\rm EC}=\frac{4}{3}\,\gamma_{\rm p}^{2}\,\nu_{\rm
ext}\,\frac{\Gamma\delta}{1+z},$ (5)

where, $\nu_{\rm ext}$ is the frequency of external photons, $\nu_{\rm
ext}=2.46\times 10^{15}\,{\rm Hz}$ for the case of external photons coming
from the broad line region (BLR) and $\nu_{\rm ext}=7.7\times 10^{13}\,{\rm
Hz}$ for the case of external photons coming from the dusty torus (DT)
(Tavecchio & Ghisellini, 2008; Ghisellini & Tavecchio, 2015).

###  3.2 The magnetic field and the Lorentz factor of electrons

The magnetic field and the electron energy are crucial to the radiation model
of blazar jets. In a leptonic model, one can estimate these two parameters,
$B$ and $\gamma_{\rm p}$, through equations 1, 4 and 5 in both SSC and EC
processes. We calculated $B$ and $\gamma_{\rm p}$ for the 1900 blazars,
including 1141 BL Lacs and 759 FSRQs, the SSC process is assumed to be the
emission case of BL Lacs and the EC process is considered for the FSRQs. An
average replacement is used for those sources without available redshift or
Doppler factors, $\langle z_{\rm F}\rangle=1.21$ and $\langle\delta_{\rm
F}\rangle=17.47$ for FSRQs, $\langle z_{\rm B}\rangle=0.528$ and
$\langle\delta_{\rm B}\rangle=11.26$ for BL Lacs. In addition, for the EC
process in FSRQs, we consider the external photons from the DT for those HSPs,
while the external photons with an origination of the BLR for those LSPs and
ISP. In this case, we have 754 FSRQs using the model of external soft photons
coming from the BLR, this is consistent with the assumption of soft photon
origin for FSRQs in literature (e.g., Tan et al., 2020), and also consistent
with the fact that the FSRQs show significant broad emission lines and the
emission lines contribute to the EC component significantly (Xiao et al.,
2022a). The rest of the 5 FSRQs are considered as HSPs, and the HSPs are
naturally considered as TeV candidates (Zhu et al., 2023). The TeV emission
could be severely absorbed by interacting with BLR soft photons, thus we
assume these 5 FSRQs with soft photons from the DT. The results of $\log B$
and $\log\gamma_{\rm p}$ are listed in columns (9) and (10) of Table 1, and
are displayed in Figure 1. The Gaussian fit is applied to the distributions,
and the fitting results give a mean value $\log B^{\rm B}=-0.51$ with a
standard deviation of 1.66 for BL Lacs and a mean value $\log B^{\rm F}=1.76$
with a standard deviation of 0.85 for FSRQs; a mean value $\log\gamma_{\rm
p}^{\rm B}=4.17$ with a standard deviation of 0.32 for BL Lacs and a mean
value $\log\gamma_{\rm p}^{\rm F}=2.01$ with a standard deviation of 0.43 for
FSRQs. A two-sample Kolmogorov–Smirnov (K–S) test, is employed to test whether
each parameter for these two subclasses is from the same parent distribution.
In the K-S test, the probability smaller than the critical value ($p=0.05$)
would be used to reject the null hypothesis, which is that the two
distributions are coming from the same parent distribution. Our results of K-S
tests show $p\sim 0$ for both $\log B$ and $\log\gamma_{\rm p}$ distributions,
suggest the BL Lac $\log B^{\rm B}$ distribution and the FSRQ $\log B^{\rm F}$
distribution come from different parent distributions, as well as for the BL
Lac $\log\gamma_{\rm p}^{\rm B}$ distribution and the FSRQ $\log\gamma_{\rm
p}^{\rm F}$ distribution.

Figure 1: The distributions of $B$ and $\gamma_{\rm p}$. The left panel gives
the distributions of the two parameters, the dashed curve stands for the
Gaussian fit of the histogram. The right panel gives the corresponding
cumulative probability distributions (CPD). The red colour stands for FSRQs
and the blue colour stands for the BL Lacs throughout this paper

.

Besides, we calculate the range of both $\log B$ and $\log\gamma_{\rm p}$ for
the 808 BCUs in our sample, the limits are obtained by assuming the seed
photons come from either the SSC or from the EC. During the calculation, the
average replacement is used for those BCUs without available redshift and
Doppler factors. The results are also listed in columns (9) and (10) of Table
1.

###  3.3 The emission region of blazars

The jet emission region size ($R$) is usually constrained by a causality
reason, which is the variability time scale, and expressed as

$R\leq c\,\Delta t\,\frac{\delta}{1+z}.$ (6)

$R$ is easily obtained if we assume a variability timescale $\Delta t$ for
those blazars with available $z$ and $\delta$. Throughout this paper, we
assume $\Delta t=1\,{\rm day}$ (Nalewajko, 2013; Fan et al., 2013) to estimate
$R$. Consider a constant and symmetric jet geometry, if the full jet cross
section is responsible for the emission region diameter, then the distance
from the central SMBH to the emission region $d_{\rm em}$ and the emission
region size $R$ are correlated as

$R=d_{\rm em}\tan\\!\phi,$ (7)

where $\phi$ is the semi-aperture opening angle of the jet (Acharyya et al.,
2021). $\phi$ is an undetectable quantity, the value of $\tan\phi=0.1$ is
fixed to study the blazar jet power in Ghisellini & Tavecchio (2009) and
$\tan\phi\lesssim 0.25$ is suggested in Dermer et al. (2009). In this work, we
can assume this opening angle $\phi$ to be close to the viewing angle $\theta$
for blazars ($\theta\simeq\phi$) due to the fact that the blazar jet is
pointing at observers. Actually, $\theta$ can be estimated by using the
apparent velocity of the resolved jet components and expressed as

$\tan\\!\phi\simeq\tan\\!\theta=\frac{2\beta_{\rm app}}{\beta^{2}_{\rm
app}+\delta^{2}-1},$ (8)

where $\beta_{\rm app}$ is the apparent velocity of the jet component
(Kellermann et al., 2004; Liodakis et al., 2018; Xiao et al., 2019), which is
usually observed via the Very Long Baseline Interferometry (VLBI) technique
(Lister et al., 2009, 2019, 2021). We manage to collect $\theta$ from Liodakis
et al. (2018) for 182 blazars (35 BL Lacs and 147 FSRQs) of our sample and
list them in column (5) of Table 2, and show them in Figure 2. We notice that
there are 142 of 147 FSRQs with $\tan\\!\theta<0.25$, taking 96.9%, there are
30 of 35 BL Lacs with $\tan\\!\theta<0.25$, taking 85.7%, . Then we calculate
the $d_{\rm em}$ for these 182 sources and list the results in column (6) of
Table 2, and show them in Figure 3. The distribution of the emission region
distance gives a mean value $\log d_{\rm em}^{\rm B}=17.37$ with a standard
deviation of 0.91 for BL Lacs and $\log d_{\rm em}^{\rm F}=17.49$ with a
standard deviation of 0.88 for FSRQs. A K-S test result of $0.31$ suggests
that the $\log d_{\rm em}^{\rm B}$ and $\log d_{\rm em}^{\rm F}$ could come
from the same distribution.

Figure 2: The distribution of $\tan\phi$ Figure 3: The distribution of $\log d_{\rm em}$ Table 2: The location of the emission region for 182 blazars. 4FGL Name | Class_1 | z | $\delta$ | $\theta$ | $\log d_{\rm em}$ | $\log d_{\rm em}\,{\rm(G15)}$ | $\log L_{\rm BLR}$
---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
J0006.3-0620 | BLL | 0.346676 | 6.96 | 8.25 | 16.97 |  | 43.60
J0017.5-0514 | FSRQ | 0.227 | 12.02 | 0.72 | 18.31 | 16.17 | 43.77
J0019.6+7327 | FSRQ | 1.781 | 7.84 | 7.32 | 16.75 |  | 45.62
J0051.1-0648 | FSRQ | 1.975 | 5.61 | 7.27 | 16.58 |  | 46.11
J0102.8+5824 | FSRQ | 0.644 | 18.51 | 2.39 | 17.84 | 16.75 | 45.04
… | … | … | … | … | … |  | …

##  4 Discussions

###  4.1 The energy density of the relativistic electron population

In the system of a blackhole-based jet, the particle-field relations depend on
uncertain jet formation, particle acceleration, and radiation mechanisms
(Dermer et al., 2014). It is natural that systems with interacting components
often tend to equipartition. As a synchrotron source, blazar jet contains
relativistic electrons with some energy density $U_{\rm e}$ and a magnetic
field whose energy density is $U_{\rm B}=B^{2}/(8\pi)$. In order to avoid
involving those poorly understood microphysics processes in the jets, Dermer
et al. (2014) assumed a condition that the equipartition between magnetic-
field and non-thermal electron energy densities holds for blazar jets, as it
was used in the analysis of a wide variety of astrophysical systems; see
Burbidge (1959) for the radio lobes, Pacholczyk (1970) for a study of radio
sources, and Beck & Krause (2005) for a theoretical equipartition study of
synchrotron observations. On large scales, the equipartition also corresponds
to the minimum jet power condition as suggested by Ghisellini & Celotti
(2001), in which the minimum jet power is required to produce the radiation we
observe at all wavelengths. The simplest equipartition relation is $U_{\rm
e}=\xi U_{\rm B}$, $\xi$ is close to the unity (Ghisellini & Celotti, 2001;
Dermer et al., 2014). It is believed that the jet is in global equipartition
between $U_{\rm e}$ and $U_{\rm B}$, however, local out-of equipartition is
also allowed to explain peculiar observation phenomenon. The very high energy
(VHE) detection of 3C 279 (MAGIC Collaboration et al., 2008) required either
emission from leptons far out of equilibrium accompanied by poor fits to the
X-ray or synchrotron data in the leptonic frame. In 2018, Fermi observed a
characteristic peak-in-peak variability pattern on time scales in minutes of
3C 279 resulting from magnetic reconnection (Shukla & Mannheim, 2020).

In Figure 1, we notice that the FSRQs show a stronger magnetic field than the
BL Lacs, while the latter shows a higher energy of electrons than the former.
This is consistent with the typical blazar radiation mechanism paradigm, in
which FSRQ contains a photon-rich environment so that the energy of
relativistic electrons in the emission blob could efficiently dissipate via
the IC process. And, BL Lacs is believed to be located in a photon-starving
environment, thus the accelerated relativistic electrons can be well
preserved. Moreover, the distribution of $\log\gamma_{\rm p}$ shows a clear
separation between the FSRQs and the BL Lacs. This result reveals the energy
of the relativistic electrons is distributed in a wide range, which contains
large variance for different types of blazars, the FSRQs taking the side of
lower energy while the BL Lacs taking the side of higher energy. We suggest
using the intersection point value, $\log\gamma_{\rm p}=3.20\pm
0.01\,(\gamma_{\rm p}\simeq 1.6\times 10^{3})$, of the two Gaussian profiles
to divide FSRQs and BL Lacs.

With known $B$, $\gamma_{\rm p}$, and $N(\gamma)$, one can obtain the magnetic
energy density $U_{\rm B}$ and the electron energy density

$U_{\rm e}=\int_{\gamma_{\rm min}}^{\gamma_{\rm max}}\gamma\,m_{\rm
e}c^{2}N(\gamma){\rm d}\gamma,$ (9)

where the minimum Lorentz $\gamma_{\rm min}=10$ and the maximum Lorentz factor
$\gamma_{\rm max}=1\times 10^{6}$ of electrons is assumed, $m_{\rm e}$ is the
rest mass of electron. Figure 4 shows the distribution of $\log U_{\rm e}$ and
the ratio $\log(U_{\rm e}/U_{\rm B})$. The Gaussian fit gives results of a
mean value $\log U_{\rm e}^{\rm B}=-1.47$ with a standard deviation of 3.12
for BL Lacs and $\log U_{\rm e}^{\rm F}=-4.40$ with a standard deviation of
-2.14 for FSRQs; a mean value $\log(U_{\rm e}^{\rm B}/U_{\rm B}^{\rm B})=0.98$
with a standard deviation of 6.23 for BL Lacs and a mean value $\log(U_{\rm
e}^{\rm F}/U_{\rm B}^{\rm F})=-6.45$ with standard deviation of 3.68 for
FSRQs. K–S tests of $\log U_{\rm e}$ and of $\log(U_{\rm e}/U_{\rm B})$ for BL
Lacs and FSRQs, and give both results of $p\sim 0$. The results suggest $\log
U_{\rm e}^{\rm B}$ and $\log U_{\rm e}^{\rm F}$ are from different
distributions, as well as for the distributions of $\log(U_{\rm e}^{\rm
B}/U_{\rm B}^{\rm B})$ and $\log(U_{\rm e}^{\rm B}/U_{\rm B}^{\rm F})$. Based
on our results shown in Figure 1 and 4, we suggest BL Lacs have an averagely
larger energy and energy density of electron distribution, meanwhile, it also
shows a larger electron-to-magnetic energy ratio than FSRQs.

Based on the distribution of $\log(U_{\rm e}/U_{\rm B})$ in Figure 4. It is
clear that most of the FSRQs are away from the value $\log(U_{\rm e}/U_{\rm
B})=0$, which is the condition of equipartition between magnetic-field and
non-thermal electron energy densities, while the BL Lacs have a mean value
near $\log(U_{\rm e}/U_{\rm B})=0$. Thus, our results suggest that the BL Lacs
stay in a ‘quasi-equipartition’ state, while the FSRQs do not.

Figure 4: The distributions of $\log U_{\rm e}$ and $\log(U_{\rm e}/U_{\rm
B})$

###  4.2 The Kelvin–Helmholtz instability

Variation, one of the characterizing properties of blazars, has been observed
across all frequencies and time scales (eg., Urry, 1996; Dermer, 1999; Fan,
1999; Singh & Meintjes, 2020; Webb et al., 2021; Otero-Santos et al., 2022;
Amaya-Almazán et al., 2022). Blazar variability time scales are observed from
years to months, to days, and even to minutes (Wagner & Witzel, 1995; Fan et
al., 1998, 2018; Aharonian et al., 2007; Albert et al., 2007). The basic idea
of the mechanism for these different time scales of variabilities is that a
disturbance created near the black hole travels outward with a Lorentz factor
$\Gamma$ and radiates energy at a distance $\gtrsim\Gamma^{2}r_{\rm g}$, where
$r_{\rm g}$ is the gravitational radius. This scenario works for most of the
cases of variability in optical, X-ray, and even $\gamma$-ray bands, together
with those models referring to the jet spiral structure, precession, or
geometric effects in the jet (Gopal-Krishna & Wiita, 1992; Camenzind &
Krockenberger, 1992), all the blazar variability seems mostly well explained.
Minutes time scale variability observed in the TeV band (Aharonian et al.,
2007; Albert et al., 2007) demands every efficient particle
acceleration/dissipation mechanism (Shukla & Mannheim, 2020; Wang et al.,
2022a) or extremely small emission zone. Although the variability mechanisms
have been explored and discussed by many people, the initial seed instability
of arising blazar variabilities, in different time scales, was few discussed.

A two-fluid model, which was proposed by (Sol et al., 1989), coupled with the
Kelvin–Helmholtz instability could have a chance. In this scenario, two fluids
(one fluid is a nonrelativistic jet with electron-proton plasma, and another
one is a relativistic jet with an electron-positron plasma) with different
speeds and compositions, the Kelvin–Helmholtz instability occurs at the
junction of the two jet components and generate significant disturbances
(Romero, 1995; Cai et al., 2022). But, this instability only happens when the
magnetic field is weaker than the critical magnetic field strength

$B_{\rm c}=[4\pi Nm_{\rm e}c^{2}(\Gamma^{2}-1)]^{1/2}\Gamma^{-1},$ (10)

where $N$ is the particle number density that we obtained via equation 3,
$\Gamma$ is assumed to be equal to $\delta$. We calculate the $B_{\rm c}$ for
FSRQs and BL Lacs in our sample, and an average value replacement is applied
to those sources without available $\delta$ during the calculation, the
results are listed in column (12) of Table 1. Based on our calculation, there
are 682 sources (taking 36% sources of the FSRQs and BL Lacs) that have a
magnetic field weaker than its critical one. The result suggests that there is
a significant number of Fermi blazars could show the Kelvin–Helmholtz
instability in the emission zone. There are possible consequences of these
instabilities. The instability could be amplified due to accumulation, thus
showing detectable variability and the variability time scale is dependent on
the energy dissipation efficiency. Or the instability could be generated in a
manner of kink, and thus show detectable violent variability. For instance,
these kink instabilities can disrupt the jet (Porth & Komissarov, 2015;
Barniol Duran et al., 2017) and trigger magnetic reconnection (Giannios &
Spruit, 2006; Shukla & Mannheim, 2020). The magnetic reconnection accelerates
particles, and these particles dissipate rapidly in a very short timescale,
e.g., the 3C 279 gamma-ray flare (Shukla & Mannheim, 2020).

With the assumption of two-fluid model Kelvin–Helmholtz instability, our
results suggest that more than one-third of the blazars are able to show this
instability. While the result is still modest if we consider observed
variabilities are initially arising from the seed instabilities because almost
all the observed blazars show strong variabilities although these
variabilities happen in different bands and at different epochs. Possible
reasons are the following (1) the SED fitting results of particle number
density, that we applied to calculate $B_{\rm c}$, are biased by high states
of blazars; (2) other mechanisms of seed instability, that were not considered
in this work, may involve. However, one should keep in mind that even if all
the blazars can generate instabilities, the strong and violent variability
only happens occasionally because most of the instabilities fade before they
are amplified to generate significant variabilities.

###  4.3 The location of the emission region

The location of the emission region is one of the essential properties in
blazars. In the framework of the one-zone leptonic model, emissions in
different bands are believed to be radiated in the same region. However, the
location of the blazar emission region is controversial. The most common
method to determine the location of the emission region is by assuming a
constant jet geometry and taking the full jet cross section as the emission
region diameter, then obtaining the distance between the emission region and
the SMBH through the basic trigonometric relations. Foschini et al. (2011)
used the $\sim$2 yr of Fermi-LAT observations to study the locations of
several blazars and found the location should be within the BLR. While the
expected spectral cut-off, arising from the photon-photon pair annihilation of
$\gamma$-rays with the helium Lyman recombination continuum within the BLR
(Poutanen & Stern, 2010), at the GeV band is not always observed. Not even the
cut-off could also be explained by the other consequence, e.g., a break on the
EED suggested by Dermer et al. (2015). Furthermore, the detection of TeV
emission from blazars suggests the TeV emission region should locate outside
the BLR because of the severe attenuation, the interactions between TeV
photons and the photons in the BLR, would not allow us to observe the TeV
emission from blazars.

In this work, we made our research on the location of the blazar emission
region. We applied a timescale of 1 day, which was suggested by Nalewajko
(2013) for as a typical variability timescale in the source frame in the Fermi
$\gamma$-ray band, this timescale has been used in many works (e.g.,
Ghisellini et al., 1998; Nalewajko, 2013; Fan et al., 2013; Chen, 2018; Pei et
al., 2022). One can set an upper limit on the emission region size via
inequality 6, thus its distance from the SMBH can be estimated via equation 7.
An approximation, which is assuming the viewing angle to be equal to the jet
semi-opening angle, is proposed and employed in this work. We notice that most
of our sources with semi-opening angles of $\tan\\!\phi<0.25$, see the upper
panel of Fig. 3. Our result is partly consistent with the range for the semi-
opening angle, which is suggested as $0.1<\tan\\!\phi<0.25$ (Ghisellini &
Tavecchio, 2009; Dermer et al., 2009). Meanwhile, there are 116 of 147 FSRQs
(taking 78.9%) and 20 of 35 BL Lacs (taking 57.1%) with $\tan\\!\phi\leq 0.1$.
In this case, our result suggests that maybe the viewing angle is smaller than
the actual semi-opening angle and the line of sight should lie within the jet
cone.

As we can see from the bottom panel of Fig. 3, the distributions show that the
emission region is located at a distance of $1.68\times 10^{14}\,{\rm cm}$ to
$6.61\times 10^{19}\,{\rm cm}$, corresponds to $5.4\times 10^{-5}\,{\rm pc}$
to $21.3\,{\rm pc}$ for all blazars. This large variance of 6 orders of
magnitude directly arises from the variance of the viewing angle, from 0.04
for 4FGL J2035.4+1056 to 84.12 for 4FGL J0113.7+0225. The mean values of $\log
d_{\rm em}^{\rm B}=17.37$ (0.076 pc) for BL Lacs and $\log d_{\rm em}^{\rm
F}=17.49$ (0.1 pc) for FSRQs are very close to each other, and the result of
K-S test confirmed that they are from the same distribution. Figure 5 shows a
comparison of $\log d_{\rm em}$ from the present work and from Ghisellini &
Tavecchio (2015), in which they estimated the $\log d_{\rm em}$ for 221
sources and gave an average value of $\sim 1\times 10^{17}\,{\rm cm}$ (0.03
pc), the comparison also shows that the $\log d_{\rm em}$ from the present
work is larger than those from their work. The comparison result could arise
from an underestimated semi-opening angle and thus an overestimated distance
from the SMBH. In addition, different methods we employed to estimate the
$\log d_{\rm em}$ could also arise the discrepancy.

Figure 5: The comparison of $\log d_{\rm em}$ from this work and from
Ghisellini & Tavecchio (2015)

Nevertheless, we check the emission region’s relative location with the BLR
and with the DT. The distances from the BLR and from the DT to the SMBH are
assumed to scale with the square root of the accretion disk luminosity,

$d_{\rm BLR}=10^{17}L_{\rm disk,\,45}^{1/2}\,{\rm cm},$ (11)

and

$d_{\rm DT}=2.5\times 10^{18}L_{\rm disk,\,45}^{1/2}\,{\rm cm},$ (12)

where the $L_{\rm disk,\,45}$ is the accretion disk luminosity $L_{\rm disk}$
in units of $10^{45}\,{\rm erg\cdot s^{-1}}$ (Ghisellini & Tavecchio, 2008).
The accretion disk luminosity can be estimated by the BLR luminosity $L_{\rm
disk}\simeq 10L_{\rm BLR}$ (Calderone et al., 2013). We manage to collect the
$L_{\rm BLR}$ from Xiao et al. (2022b) for 143 of the 182 blazars and list
them in Table 2, we find that there are 66 (taking 46.2%) sources with
emission regions located within the BLR, 63 (taking 44.0%) sources with
emission regions located between the BLR and the DT, and 14 (9.8%) sources
with emission regions located beyond the DT. Locating the emission region much
further out of the DT is not appropriate because there are no important
sources of external photons. The reason for the 14 sources with overestimated
$d_{\rm em}$ is that small semi-opening angles, less than 0.1, are employed.

To sum up, our results suggest about half of the blazars with the emitting
regions $d_{\rm em}<d_{\rm BLR}$ and another half with the emitting regions
$d_{\rm BLR}<d_{\rm em}<d_{\rm DT}$. This is clearly inconsistent with the
result suggested in Ghisellini et al. (2014), in which they suggested 85%
sources with emission region located within the BLR and 15% sources with
emission region located between the BLR and the DT. This discrepancy could be
caused mainly by two reasons, namely, the sample size and the methods of
estimating the $d_{\rm em}$. While the method influence can not be the main
reason because the deviation caused by one particular method should evenly
affect $d_{\rm em}$ for all sources in the sample and should not significantly
change the portion of $d_{\rm em}$ within the BLR or within the DT. The best
way to eliminate this discrepancy is to compile a larger/complete sample of
blazars to estimate the distance from the SMBH to the emission region and
determine the emission region’s relative location to the BLR or the DT.

##  5 Conclusions

In order to study the property of the blazar emission region and put
constraints on the corresponding parameters, we compiled a sample of 2708
Fermi blazars with available broadband SED features. A fraction of the sources
in our sample also has available redshift, Doppler factor, viewing angle, and
BLR luminosity. With the above-mentioned information, we calculated the
magnetic field ($\log B$) and electron energy ($\log\gamma_{\rm p}$) for
certain type blazars (FSRQs and BL Lacs) and provided ranges of $\log B$ and
$\log\gamma_{\rm p}$ for BCUs, electron energy density ($\log U_{\rm e}$), the
energy ratio ($\log(U_{\rm e}/U_{\rm B})$), and the critical magnetic strength
($B_{\rm c}$) to study the jet launching and seed variability. A distance
($d_{\rm em}$) from the SMBH to the emission region is obtained for a
subsample of 182 sources and used to discuss the emission region location
compare to the BLR.

Our main conclusions are as follows: (1) The FSRQs show lower electron energy
than the BL Lacs, we suggest to divided FSRQs and BL Lacs with
$\log\gamma_{\rm p}=3.20\pm 0.01$. Besides, FSRQs show trends of a stronger
magnetic field and smaller electron-to-magnetic energy ratio than the BL Lacs;
(2) Our results suggest that the BL Lacs with a mean value of $\log(U^{\rm
B}_{\rm e}/U^{\rm B}_{\rm B})=0.98$, may fulfil the equipartition between
magnetic field energy density and the electron energy density, while the FSRQs
with a mean value of $\log(U^{\rm F}_{\rm e}/U^{\rm F}_{\rm B})=-6.45$ is away
from the condition $\log(U_{\rm e}/U_{\rm B})\sim 0$; (3) Comparing the $B$
and $B_{\rm c}$, we find 682 blazars with $B$ samller than the $B_{\rm c}$ and
suggest these sources are all candidates of showing Kelvin–Helmholtz
instability. But we note that our result is an old-line and the rest of the
blazars could also have a chance to show this instability; (4) Our result of
the location of the blazar emission region is particularly novel. The result
suggests both FSRQs and BL Lacs with emission region at the distance of
$\sim$0.1 pc. And we find that about half of the blazars with emission regions
within the BLR, while another half sources with the region located between the
BLR and the DT.

J. H. Fan thanks the support from the NSFC (NSFC U2031201, NSFC 11733001,
U2031112), Scientific and Technological Cooperation Projects (2020-2023)
between the People’s Republic of China and the Republic of Bulgaria, Guangdong
Major Project of Basic and Applied Basic Research (Grant No. 2019B030302001),
the science research grants from the China Manned Space Project with NO. CMS-
CSST-2021-A06, and the support for Astrophysics Key Subjects of Guangdong
Province and Guangzhou City. This research was partially supported by the
Bulgarian National Science Fund of the Ministry of Education and Science under
grants KP-06-H38/4 (2019), KP-06-KITAJ/2 (2020) and KP-06-H68/4 (2022). H.B.
Xiao acknowledges the support from the National Natural Science Foundation of
China (NSFC 12203034) and from the Shanghai Science and Technology Fund
(22YF1431500).

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30, doi: 10.1088/0004-637X/716/1/30
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Acharyya et al. (2021) Acharyya, A., Chadwick, P. M., & Brown, A. M. 2021, MNRAS, 500, 5297, doi: 10.1093/mnras/staa3483
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJ, 664, L71, doi: 10.1086/520635
  * Albert et al. (2007) Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862, doi: 10.1086/521382
  * Amaya-Almazán et al. (2022) Amaya-Almazán, R. A., Chavushyan, V., & Patiño-Álvarez, V. M. 2022, ApJ, 929, 14, doi: 10.3847/1538-4357/ac5741
  * Barniol Duran et al. (2017) Barniol Duran, R., Tchekhovskoy, A., & Giannios, D. 2017, MNRAS, 469, 4957, doi: 10.1093/mnras/stx1165
  * Beck & Krause (2005) Beck, R., & Krause, M. 2005, Astronomische Nachrichten, 326, 414, doi: 10.1002/asna.200510366
  * Blandford & Koenigl (1979) Blandford, R. D., & Koenigl, A. 1979, Astrophys. Lett., 20, 15
  * Blandford et al. (1977) Blandford, R. D., McKee, C. F., & Rees, M. J. 1977, Nature, 267, 211, doi: 10.1038/267211a0
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Burbidge (1959) Burbidge, G. R. 1959, ApJ, 129, 849, doi: 10.1086/146680
  * Cai et al. (2022) Cai, J. T., Kurtanidze, S. O., Liu, Y., et al. 2022, ApJS, 260, 47, doi: 10.3847/1538-4365/ac666b
  * Calderone et al. (2013) Calderone, G., Ghisellini, G., Colpi, M., & Dotti, M. 2013, MNRAS, 431, 210, doi: 10.1093/mnras/stt157
  * Camenzind & Krockenberger (1992) Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59
  * Cao & Spruit (2013) Cao, X., & Spruit, H. C. 2013, ApJ, 765, 149, doi: 10.1088/0004-637X/765/2/149
  * Cerruti et al. (2015) Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910, doi: 10.1093/mnras/stu2691
  * Chen (2014) Chen, L. 2014, ApJ, 788, 179, doi: 10.1088/0004-637X/788/2/179
  * Chen (2017) —. 2017, ApJ, 842, 129, doi: 10.3847/1538-4357/aa7744
  * Chen (2018) —. 2018, ApJS, 235, 39, doi: 10.3847/1538-4365/aab8fb
  * Dermer (1999) Dermer, C. D. 1999, Astroparticle Physics, 11, 1, doi: 10.1016/S0927-6505(99)00017-1
  * Dermer et al. (2014) Dermer, C. D., Cerruti, M., Lott, B., Boisson, C., & Zech, A. 2014, ApJ, 782, 82, doi: 10.1088/0004-637X/782/2/82
  * Dermer et al. (2009) Dermer, C. D., Finke, J. D., Krug, H., & Böttcher, M. 2009, ApJ, 692, 32, doi: 10.1088/0004-637X/692/1/32
  * Dermer et al. (2015) Dermer, C. D., Yan, D., Zhang, L., Finke, J. D., & Lott, B. 2015, ApJ, 809, 174, doi: 10.1088/0004-637X/809/2/174
  * Diltz et al. (2015) Diltz, C., Böttcher, M., & Fossati, G. 2015, ApJ, 802, 133, doi: 10.1088/0004-637X/802/2/133
  * Dimitrakoudis et al. (2012) Dimitrakoudis, S., Mastichiadis, A., Protheroe, R. J., & Reimer, A. 2012, A&A, 546, A120, doi: 10.1051/0004-6361/201219770
  * Fan (1999) Fan, J. H. 1999, MNRAS, 308, 1032, doi: 10.1046/j.1365-8711.1999.02767.x
  * Fan (2002) Fan, J.-H. 2002, PASJ, 54, L55, doi: 10.1093/pasj/54.4.L55
  * Fan (2005) —. 2005, Chinese Journal of Astronomy and Astrophysics Supplement, 5, 213
  * Fan et al. (2014) Fan, J.-H., Bastieri, D., Yang, J.-H., et al. 2014, Research in Astronomy and Astrophysics, 14, 1135, doi: 10.1088/1674-4527/14/9/004
  * Fan et al. (1998) Fan, J. H., Xie, G. Z., Pecontal, E., Pecontal, A., & Copin, Y. 1998, ApJ, 507, 173, doi: 10.1086/306301
  * Fan et al. (2013) Fan, J.-H., Yang, J.-H., Liu, Y., & Zhang, J.-Y. 2013, Research in Astronomy and Astrophysics, 13, 259, doi: 10.1088/1674-4527/13/3/002
  * Fan et al. (2016) Fan, J. H., Yang, J. H., Liu, Y., et al. 2016, ApJS, 226, 20, doi: 10.3847/0067-0049/226/2/20
  * Fan et al. (2018) Fan, J. H., Tao, J., Liu, Y., et al. 2018, AJ, 155, 90, doi: 10.3847/1538-3881/aaa547
  * Fan et al. (2021) Fan, J. H., Kurtanidze, S. O., Liu, Y., et al. 2021, ApJS, 253, 10, doi: 10.3847/1538-4365/abd32d
  * Foschini et al. (2011) Foschini, L., Ghisellini, G., Tavecchio, F., Bonnoli, G., & Stamerra, A. 2011, A&A, 530, A77, doi: 10.1051/0004-6361/201117064
  * Ghisellini & Celotti (2001) Ghisellini, G., & Celotti, A. 2001, MNRAS, 327, 739, doi: 10.1046/j.1365-8711.2001.04700.x
  * Ghisellini et al. (1998) Ghisellini, G., Celotti, A., Fossati, G., Maraschi, L., & Comastri, A. 1998, MNRAS, 301, 451, doi: 10.1046/j.1365-8711.1998.02032.x
  * Ghisellini et al. (1993) Ghisellini, G., Padovani, P., Celotti, A., & Maraschi, L. 1993, ApJ, 407, 65, doi: 10.1086/172493
  * Ghisellini & Tavecchio (2008) Ghisellini, G., & Tavecchio, F. 2008, MNRAS, 387, 1669, doi: 10.1111/j.1365-2966.2008.13360.x
  * Ghisellini & Tavecchio (2009) —. 2009, MNRAS, 397, 985, doi: 10.1111/j.1365-2966.2009.15007.x
  * Ghisellini & Tavecchio (2015) —. 2015, MNRAS, 448, 1060, doi: 10.1093/mnras/stv055
  * Ghisellini et al. (2014) Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Nature, 515, 376, doi: 10.1038/nature13856
  * Giannios & Spruit (2006) Giannios, D., & Spruit, H. C. 2006, A&A, 450, 887, doi: 10.1051/0004-6361:20054107
  * Gopal-Krishna & Wiita (1992) Gopal-Krishna, & Wiita, P. J. 1992, A&A, 259, 109
  * Gupta et al. (2016) Gupta, A. C., Agarwal, A., Bhagwan, J., et al. 2016, MNRAS, 458, 1127, doi: 10.1093/mnras/stw377
  * H. E. S. S. Collaboration et al. (2015) H. E. S. S. Collaboration, Abramowski, A., Aharonian, F., et al. 2015, A&A, 573, A31, doi: 10.1051/0004-6361/201321436
  * Hovatta et al. (2019) Hovatta, T., O’Sullivan, S., Martí-Vidal, I., Savolainen, T., & Tchekhovskoy, A. 2019, A&A, 623, A111, doi: 10.1051/0004-6361/201832587
  * IceCube Collaboration et al. (2018) IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * Kellermann et al. (1989) Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989, AJ, 98, 1195, doi: 10.1086/115207
  * Kellermann et al. (2004) Kellermann, K. I., Lister, M. L., Homan, D. C., et al. 2004, ApJ, 609, 539, doi: 10.1086/421289
  * Kubo et al. (1998) Kubo, H., Takahashi, T., Madejski, G., et al. 1998, ApJ, 504, 693, doi: 10.1086/306125
  * Liodakis et al. (2022) Liodakis, I., Blinov, D., Potter, S. B., & Rieger, F. M. 2022, MNRAS, 509, L21, doi: 10.1093/mnrasl/slab118
  * Liodakis et al. (2018) Liodakis, I., Hovatta, T., Huppenkothen, D., et al. 2018, ApJ, 866, 137, doi: 10.3847/1538-4357/aae2b7
  * Lister et al. (2015) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2015, ApJ, 810, L9, doi: 10.1088/2041-8205/810/1/L9
  * Lister et al. (2021) Lister, M. L., Homan, D. C., Kellermann, K. I., et al. 2021, ApJ, 923, 30, doi: 10.3847/1538-4357/ac230f
  * Lister et al. (2009) Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874, doi: 10.1088/0004-6256/138/6/1874
  * Lister et al. (2019) Lister, M. L., Homan, D. C., Hovatta, T., et al. 2019, ApJ, 874, 43, doi: 10.3847/1538-4357/ab08ee
  * Lynden-Bell (1969) Lynden-Bell, D. 1969, Nature, 223, 690, doi: 10.1038/223690a0
  * MAGIC Collaboration et al. (2008) MAGIC Collaboration, Albert, J., Aliu, E., et al. 2008, Science, 320, 1752, doi: 10.1126/science.1157087
  * Massaro et al. (2004) Massaro, E., Perri, M., Giommi, P., & Nesci, R. 2004, A&A, 413, 489, doi: 10.1051/0004-6361:20031558
  * Massaro et al. (2006) Massaro, E., Tramacere, A., Perri, M., Giommi, P., & Tosti, G. 2006, A&A, 448, 861, doi: 10.1051/0004-6361:20053644
  * Mücke & Protheroe (2001) Mücke, A., & Protheroe, R. J. 2001, Astroparticle Physics, 15, 121, doi: 10.1016/S0927-6505(00)00141-9
  * Nalewajko (2013) Nalewajko, K. 2013, MNRAS, 430, 1324, doi: 10.1093/mnras/sts711
  * Otero-Santos et al. (2022) Otero-Santos, J., Acosta-Pulido, J. A., Becerra González, J., et al. 2022, MNRAS, 511, 5611, doi: 10.1093/mnras/stac475
  * Pacholczyk (1970) Pacholczyk, A. G. 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources
  * Paliya et al. (2021) Paliya, V. S., Domínguez, A., Ajello, M., Olmo-García, A., & Hartmann, D. 2021, ApJS, 253, 46, doi: 10.3847/1538-4365/abe135
  * Pei et al. (2022) Pei, Z., Fan, J., Yang, J., Huang, D., & Li, Z. 2022, ApJ, 925, 97, doi: 10.3847/1538-4357/ac3aeb
  * Pei et al. (2016) Pei, Z.-Y., Fan, J.-H., Liu, Y., et al. 2016, Ap&SS, 361, 237, doi: 10.1007/s10509-016-2822-0
  * Porth & Komissarov (2015) Porth, O., & Komissarov, S. S. 2015, MNRAS, 452, 1089, doi: 10.1093/mnras/stv1295
  * Poutanen & Stern (2010) Poutanen, J., & Stern, B. 2010, ApJ, 717, L118, doi: 10.1088/2041-8205/717/2/L118
  * Readhead (1994) Readhead, A. C. S. 1994, ApJ, 426, 51, doi: 10.1086/174038
  * Rees (1984) Rees, M. J. 1984, ARA&A, 22, 471, doi: 10.1146/annurev.aa.22.090184.002351
  * Romero (1995) Romero, G. E. 1995, Ap&SS, 234, 49, doi: 10.1007/BF00627281
  * Scarpa & Falomo (1997) Scarpa, R., & Falomo, R. 1997, A&A, 325, 109
  * Shukla & Mannheim (2020) Shukla, A., & Mannheim, K. 2020, Nature Communications, 11, 4176, doi: 10.1038/s41467-020-17912-z
  * Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153, doi: 10.1086/173633
  * Singh & Meintjes (2020) Singh, K. K., & Meintjes, P. J. 2020, Astronomische Nachrichten, 341, 713, doi: 10.1002/asna.202013731
  * Sokolov & Marscher (2005) Sokolov, A., & Marscher, A. P. 2005, ApJ, 629, 52, doi: 10.1086/431321
  * Sol et al. (1989) Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411, doi: 10.1093/mnras/237.2.411
  * Strittmatter et al. (1980) Strittmatter, P. A., Hill, P., Pauliny-Toth, I. I. K., Steppe, H., & Witzel, A. 1980, A&A, 88, L12
  * Tan et al. (2020) Tan, C., Xue, R., Du, L.-M., et al. 2020, ApJS, 248, 27, doi: 10.3847/1538-4365/ab8cc6
  * Tavecchio & Ghisellini (2008) Tavecchio, F., & Ghisellini, G. 2008, MNRAS, 386, 945, doi: 10.1111/j.1365-2966.2008.13072.x
  * Tavecchio et al. (1998) Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608, doi: 10.1086/306526
  * Tramacere et al. (2011) Tramacere, A., Massaro, E., & Taylor, A. M. 2011, ApJ, 739, 66, doi: 10.1088/0004-637X/739/2/66
  * Urry (1996) Urry, C. M. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 110, Blazar Continuum Variability, ed. H. R. Miller, J. R. Webb, & J. C. Noble, 391, doi: 10.48550/arXiv.astro-ph/9609023
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Villata et al. (2006) Villata, M., Raiteri, C. M., Balonek, T. J., et al. 2006, A&A, 453, 817, doi: 10.1051/0004-6361:20064817
  * Wagner & Witzel (1995) Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163, doi: 10.1146/annurev.aa.33.090195.001115
  * Wang et al. (2022a) Wang, G., Fan, J., Xiao, H., & Cai, J. 2022a, PASP, 134, 104101, doi: 10.1088/1538-3873/ac98e0
  * Wang et al. (2022b) Wang, Z.-R., Liu, R.-Y., Petropoulou, M., et al. 2022b, Phys. Rev. D, 105, 023005, doi: 10.1103/PhysRevD.105.023005
  * Wardle (2018) Wardle, J. 2018, Galaxies, 6, 5, doi: 10.3390/galaxies6010005
  * Webb et al. (2021) Webb, J. R., Arroyave, V., Laurence, D., et al. 2021, Galaxies, 9, 114, doi: 10.3390/galaxies9040114
  * Wills et al. (1992) Wills, B. J., Wills, D., Breger, M., Antonucci, R. R. J., & Barvainis, R. 1992, ApJ, 398, 454, doi: 10.1086/171869
  * Xiao et al. (2022a) Xiao, H., Fan, J., Ouyang, Z., et al. 2022a, ApJ, 936, 146, doi: 10.3847/1538-4357/ac887f
  * Xiao et al. (2020) Xiao, H., Fan, J., Rando, R., Zhu, J., & Hu, L. 2020, Astronomische Nachrichten, 341, 462, doi: 10.1002/asna.202013733
  * Xiao et al. (2022b) Xiao, H., Ouyang, Z., Zhang, L., et al. 2022b, ApJ, 925, 40, doi: 10.3847/1538-4357/ac36da
  * Xiao et al. (2022c) Xiao, H., Zhu, J., Fu, L., Zhang, S., & Fan, J. 2022c, PASJ, 74, 239, doi: 10.1093/pasj/psab121
  * Xiao et al. (2019) Xiao, H., Fan, J., Yang, J., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 129811, doi: 10.1007/s11433-018-9371-x
  * Xiao et al. (2022d) Xiao, H. B., Zhu, J. T., Fan, J. H., et al. 2022d, MNRAS, 517, 4202, doi: 10.1093/mnras/stac2959
  * Xiong & Zhang (2014) Xiong, D. R., & Zhang, X. 2014, MNRAS, 441, 3375, doi: 10.1093/mnras/stu755
  * Xue et al. (2019) Xue, R., Liu, R.-Y., Wang, X.-Y., Yan, H., & Böttcher, M. 2019, ApJ, 871, 81, doi: 10.3847/1538-4357/aaf720
  * Xue et al. (2021) Xue, R., Liu, R.-Y., Wang, Z.-R., Ding, N., & Wang, X.-Y. 2021, ApJ, 906, 51, doi: 10.3847/1538-4357/abc886
  * Yang et al. (2023) Yang, J., Fan, J., Liu, Y., et al. 2023, Science China Physics, Mechanics, and Astronomy, 66, 249511, doi: 10.1007/s11433-022-2062-2
  * Yang et al. (2022) Yang, J. H., Fan, J. H., Liu, Y., et al. 2022, ApJS, 262, 18, doi: 10.3847/1538-4365/ac7deb
  * Zhang et al. (2022) Zhang, L., Liu, Y., & Fan, J. 2022, ApJ, 935, 4, doi: 10.3847/1538-4357/ac7bde
  * Zhang & Fan (2008) Zhang, Y.-W., & Fan, J.-H. 2008, Chinese J. Astron. Astrophys., 8, 385, doi: 10.1088/1009-9271/8/4/02
  * Zhu et al. (2023) Zhu, J. T., Lin, C., Xiao, H. B., et al. 2023, arXiv e-prints, arXiv:2303.10557, doi: 10.48550/arXiv.2303.10557


# Enhanced Next-to-Leading-Order Corrections to Weak Annihilation $B$-Meson
Decays

Cai-Dian Lüa,b<EMAIL_ADDRESS>Yue-Long Shenc corresponding author:
<EMAIL_ADDRESS>Chao Wangd corresponding author:
<EMAIL_ADDRESS>Yu-Ming Wange corresponding author:
<EMAIL_ADDRESS>a Institute of High Energy Physics, CAS, P.O. Box
918(4), Beijing 100049, China
b School of Physics, University of Chinese Academy of Sciences, Beijing
100049, China
c School of Physics and Optoelectronic Engineering, Ocean University of China,
Songling Road 238, Qingdao, Shandong 266100, P.R. China
d Department of Mathematics and Physics, Huaiyin Institute of Technology,
Meicheng East Road 1, Huaian, Jiangsu 223200, P.R. China
e School of Physics, Nankai University, Weijin Road 94, Tianjin 300071, P.R.
China

###### Abstract

We accomplish the analytical computation of the pure weak annihilation non-
leptonic $B$-meson decay amplitudes at leading power in the heavy quark
expansion. The novel observation regarding such fundamental hadronic
quantities is that adding the missing hard-collinear contribution on top of
the hard gluon exchange effect eliminates rapidity divergences entering the
convolution integrals of factorization formulae. Subsequently we identify the
perturbative enhancement mechanism due to the penguin contractions of the
current-current operators from the weak effective Hamiltonian, which yields
the significant impacts on the CP violating observables.

##  I Introduction

It is generally accepted that the exclusive two-body charmless bottom-meson
decays are of fundamental importance for advancing our understanding towards
diverse facets of the strong interaction dynamics governing the flavour-
changing heavy-quark decay processes and for exploring the peculiar
implications of the Cabibbo-Kobayashi-Maskawa (CKM) mechanism for CP violation
in electroweak interactions. The primary challenge of predicting such hadronic
$B$-meson decay observables consists in constructing a field-theoretic
framework to disentangle the calculable perturbative QCD fluctuations from the
non-perturbative soft and collinear physics systematically. In this respect,
the QCD factorization formalism based upon expansions in the small parameter
$\Lambda_{\rm QCD}/m_{b}$ Beneke  _et al._ (1999, 2000); Chay and Kim (2004);
Beneke  _et al._ (2002); Bauer  _et al._ (2004, 2006) has proven to provide
the precise prescriptions for evaluating the appeared non-leptonic matrix
elements of the effective weak Hamiltonian operators. However, the long-
standing obstacle to improve the factorization calculations of the two-body
$B$-meson decay amplitudes arises from the disturbing end-point divergences in
the convolution integrals describing the power-suppressed but
phenomenologically pronounced weak annihilation corrections. In particular,
the inadequacy of achieving the factorization-compatible regularization of the
annihilation amplitudes apparently does not allow to provide quantitative
predictions for the annihilation-dominated decay processes.

Additionally, the abiding lack of higher-order perturbative corrections to the
weak annihilation topologies prevents us from obtaining reliable estimates of
the strong-phase sensitive quantities (for instance, the direct CP asymmetry)
due to a variety of potential enhancements, which do not manifest themselves
at leading-order (LO) in the strong coupling. The process-independent
parametrizations for logarithmically and linearly divergent integrals in the
basic building blocks for the exclusive hadronic $\bar{B}_{q}\to PP$,
$\bar{B}_{q}\to PV$ and $\bar{B}_{q}\to VV$ (with $q=d,\,s$) decay amplitudes
Beneke  _et al._ (2001a); Beneke and Neubert (2003); Beneke  _et al._ (2007);
Bartsch  _et al._ (2008) (see also Cheng and Yang (2008); Cheng and Chua
(2009a, b) and references therein) introduce further assumption for the
current QCD factorization predictions, especially with the finite bottom-quark
mass in practice, which has been examined exploratorily by employing the data-
driven strategy Zhu (2011); Wang and Zhu (2013); Bobeth  _et al._ (2015);
Chang  _et al._ (2015, 2016, 2017). The striking importance of the weak
annihilation mechanism in investigating the non-leptonic $B$-meson decay
amplitudes also motivated the interesting dynamical anatomy Arnesen  _et al._
(2008); Khodjamirian  _et al._ (2005); Feldmann and Hurth (2004) on the basis
of distinct QCD techniques (see Keum  _et al._ (2001); Lu  _et al._ (2001); Li
_et al._ (2004); Ali  _et al._ (2007); Li  _et al._ (2011, 2012, 2013, 2014);
Li and Wang (2015) for an alternative treatment taking into account the non-
vanishing transverse momenta of the active partons). As a consequence,
accomplishing the systematic computation of the weak annihilation topologies
with the non-trivial strong phases will be of the top priority for developing
the conceptual framework of the charmless two-body $\bar{B}_{q}\to M_{1}M_{2}$
decays and for better confronting the anticipated precision measurements of a
large number of CP violating observables at the Belle II experiment
Altmannshofer  _et al._ (2019).

In order to probe factorization properties of the weak annihilation
contributions immaculately, we report on a novel observation of the hard-
collinear gluon exchange (the short-distance configuration with the invariant
mass of ${\cal O}(\sqrt{m_{b}\,\Lambda_{\rm QCD}}))$ in regularizing the
unwanted rapidity divergences analytically, which has been unfortunately
neglected in all previous calculations of the two-body hadronic $B$-meson
decays, by concentrating on the pure annihilation channels such as
$\bar{B}_{s}\to\pi\pi$ and $\bar{B}_{d}\to\phi\phi$. The significance of this
essential insight for consolidating the theory foundation of QCD factorization
in $B$-meson decaying into two light mesons will be further strengthened by
inspecting a sample set of the next-to-leading order (NLO) Feynman diagrams,
which are enhanced by the large Wilson coefficients and/or the multiplication
CKM factors numerically and are promising to yield the sizeable strong phases.
General implications of the improved formalism on the accessible decay
observables will be then discussed with the acceptable models for the leading-
twist light-cone distribution amplitudes of both the bottom meson and the
final-state hadrons.

##  II General analysis

We start by setting up our notation of the effective weak Hamiltonian of the
non-leptonic $\Delta B=1$ transitions in the Standard Model (SM)

$\displaystyle{\cal H}_{\rm eff}$ $\displaystyle=$
$\displaystyle{4\,G_{F}\over\sqrt{2}}\,\sum_{p=u,c}V_{pb}\,V_{pq}^{\ast}\bigg{[}C_{1}(\nu)\,P_{1}^{p}(\nu)+C_{2}(\nu)\,P_{2}^{p}(\nu)$
(1)
$\displaystyle+\sum_{i=3}^{6}C_{i}(\nu)\,P_{i}(\nu)+\sum_{i=3}^{6}C_{iQ}(\nu)\,P_{iQ}(\nu)$
$\displaystyle+C_{7\gamma}(\nu)\,\,P_{7\gamma}(\nu)+C_{8g}(\nu)\,\,P_{8g}(\nu)\bigg{]}\,,$

where we adopt the effective operator basis introduced in Chetyrkin  _et al._
(1997, 1998); Bobeth  _et al._ (2000, 2004); Huber  _et al._ (2006) enabling
the disappearance of Dirac traces with $\gamma_{5}$. The renormalized Wilson
coefficients $C_{i}(\nu)$ with $\nu\sim{\cal O}(m_{b})$ will be evaluated for
our purpose in the next-to-leading-logarithmic (NLL) approximation in the
Chetyrkin-Misiak-Münz basis Chetyrkin  _et al._ (1998).

It proves convenient to cast the exclusive transition amplitude governing the
weak-annihilation bottom-meson decay in the form

$\displaystyle\bar{\cal A}(\bar{B}_{q}\to M_{1}M_{2})=-\langle
M_{1}(p_{1})\,M_{2}(p_{2})|{\cal H}_{\rm eff}|\bar{B}_{q}(p_{B})\rangle$
$\displaystyle=\mp\,i\,{4\,G_{F}\over\sqrt{2}}\,f_{B}\,f_{M_{1}}\,f_{M_{2}}\,\sum_{p=u,c}V_{pb}\,V_{pq}^{\ast}\,\left(\pi\,\alpha_{s}\right)\,{C_{F}\over
N_{c}^{2}}\,$ $\displaystyle\hskip 14.22636pt\times\left[{\cal
T}^{p,\,(0)}+\left({\alpha_{s}\over 4\,\pi}\right)\,{\cal T}^{p,\,(1)}+{\cal
O}(\alpha_{s}^{2})\right]\,,\hskip 14.22636pt$ (2)

where the upper sign in (2) applies for $M_{1}\,M_{2}=PP,\,V_{L}V_{L}$, while
the lower sign when $M_{1}\,M_{2}=PV,\,VP$. We leave out the discussion on
final states consisting of two transversely polarized vector mesons, bearing
in mind that the resulting helicity amplitudes are power-suppressed when
compared with the counterpart longitudinal polarization amplitude Beneke  _et
al._ (2007). The appearing LO quantities ${\cal T}^{p,\,(0)}$ can be
determined by evaluating the first two Feynman diagrams in Figure 1, where we
also display the relevant diagrams generating the dynamically enhanced
contributions to the NLO pieces ${\cal T}^{p,\,(1)}$. It remains important to
remark that the additional LO annihilation diagrams with the gluon radiation
off the final-state partons (not included in Figure 1) lead to the vanishing
effects for the pure annihilation channels in the limit of symmetric
distribution amplitudes and under the assumption of the ${\rm SU(3)}$ flavour
symmetry as noted in Beneke  _et al._ (2000).

Figure 1: Sample Feynman diagrams for the weak annihilation $\bar{B}_{q}\to
M_{1}M_{2}$ decays at LO and NLO. The gluons marked with red colour can carry
either the hard or hard-collinear momentum, while the remaining gluons can
only possess the hard momentum. Symmetric diagrams that follow from (c)-(g) by
exchanging the two off-shell gluons are not displayed.

Taking advantage of the flavour-decomposition strategy Beneke and Neubert
(2003) allows us to write down

$\displaystyle{\cal T}^{p,\,(0)}$ $\displaystyle=$
$\displaystyle\delta_{pu}\,{\cal C}_{M_{1}M_{2}}^{(1)}\,{\cal
B}_{1}(M_{1}M_{2})+\,{\cal C}_{M_{1}M_{2}}^{p,\,(4)}\,{\cal
B}_{4}(M_{1}M_{2})$ (3) $\displaystyle+\,{\cal C}_{M_{1}M_{2}}^{p,\,(4,{\rm
EW})}\,{\cal B}_{4,\rm{EW}}(M_{1}M_{2})\,.$

The prefactors ${\cal C}_{M_{1}M_{2}}^{(p),\,(i,\,(\rm EW))}$ collect the
Clebsh-Gordan coefficients from the flavour structures of the
$\bar{B}_{q}$-meson as well as the final-state mesons and they further absorb
the electric-charge coefficients from $P_{iQ}$. The phenomenologically
dominating effects of ${\cal T}^{p,\,(1)}$ can be categorized in terms of
their topological structures

$\displaystyle{\cal T}^{p,\,(1)}$ $\displaystyle\supset$
$\displaystyle\sum_{i=1}^{6}C_{i}\,{\cal P}_{i}^{p,(1)}(M_{1}M_{2})+C_{8}^{\rm
eff}\,{\cal P}_{8g}^{(1)}(M_{1}M_{2})$ (4)
$\displaystyle+\sum_{i=1}^{6}C_{i}\,{\cal I}_{i}^{p,(1)}(M_{1}M_{2})\,,$

with $C_{8}^{\rm eff}=C_{8g}+C_{3}-C_{4}/6+20\,C_{5}-10C_{6}/3$ Beneke  _et
al._ (2001b). We will dedicate the next section to the analytical computation
of the flavour amplitudes ${\cal B}_{i,({\rm EW})}$ and the NLO building
blocks ${\cal P}_{1...6,8g}^{p,(1)}$ and ${\cal I}_{i}^{p,(1)}$ within QCD
factorization.

##  III QCD factorization for decay amplitudes

Establishing the factorization formulae for the hadronic quantities ${\cal
B}_{i,({\rm EW})}$ can be customarily achieved by investigating the
appropriate partonic amplitudes of the tree diagrams (a) and (b) in Figure 1
in the leading-twist approximation. The yielding contribution due to the gluon
emission from the bottom quark (i.e., the LO diagram (b)) has been
demonstrated to be calculable self-consistently at twist-two order Chernyak
and Zhitnitsky (1982); Beneke  _et al._ (2000). We are then led to anatomize
the infrared behaviour of the effective matrix element represented by the
first diagram in Figure 1 (with an insertion of $P_{2}^{u}$ for the
illustration purpose)

$\displaystyle\langle P_{2}^{u}\rangle^{(0)}_{(a)}$ $\displaystyle=$
$\displaystyle\left(\pi\,\alpha_{s}\right)\,\frac{C_{F}}{N_{c}^{2}}\,\int_{0}^{\infty}d\omega\,\int_{0}^{1}dx\,\int_{0}^{1}dy\,$
(5) $\displaystyle{\rm Tr}\bigg{\\{}{\cal
M}^{B}(v,\omega)\,\left[\gamma_{\perp\nu}\,\left(\not{p}_{1}+\not{q}_{2}-\not{k}\right)\gamma_{\mu}\left(1-\gamma_{5}\right)\right]\,$
$\displaystyle\hskip 14.22636pt{\cal
M}^{M_{1}}(p,y)\,\left[\gamma_{\perp}^{\nu}\right]\,{\cal
M}^{M_{2}}(q,x)\,\left[\gamma^{\mu}\left(1-\gamma_{5}\right)\right]\bigg{\\}}$
$\displaystyle\times\,\frac{1}{[(p_{1}+q_{2}-k)^{2}+i\epsilon][(p_{1}+q_{2})^{2}+i\epsilon]}\,,$

where the explicit expressions of the momentum-space projection operators
${\cal M}^{B},\,{\cal M}^{M_{1}}$ and ${\cal M}^{M_{2}}$ have been derived in
Beneke and Feldmann (2001); Beneke  _et al._ (2001a); Beneke (2002); Beneke
_et al._ (2007) (see Grozin and Neubert (1997) for the original discussion on
$B$-meson distribution amplitudes in heavy quark effective theory (HQET) and
Bell  _et al._ (2013) for an alternative construction of the on-shell light-
cone projectors). Here we have assigned the four-momenta $p_{1,2}$ ($q_{1,2}$)
to the quark and antiquark in the energetic meson $M_{1}$ ($M_{2}$). The
dimensionless variable $y$ ($x$) refers to the momentum fraction carried by
the quark field in the composite $M_{1}$ ($M_{2}$) system in the
(anti)-collinear limit. In addition, the longitudinal momentum component
$\omega$ for the light partonic constituent of the $B$-meson state can be
defined by the kinematic relation $p_{1}\cdot k=\bar{n}\cdot p_{1}\,n\cdot
k/2\equiv\bar{n}\cdot p_{1}\,\omega/2$ with the collinear-$p_{1}$
approximation, where the two light-cone vectors $n_{\mu}$ and $\bar{n}_{\mu}$
fulfill the constraints $n^{2}=\bar{n}^{2}=0$ and $n\cdot\bar{n}=2$. It is
then convenient to simplify the quark/gluon propagators in the last line of
(5) to

$\displaystyle\frac{1}{y\,m_{B}^{2}\,\left[\bar{x}-{\omega/m_{B}}+i\epsilon\right]}\,\frac{1}{y\,\bar{x}\,m_{B}^{2}+i\epsilon}\,,$
(6)

by applying the kinematic properties of $\bar{B}_{q}\to M_{1}M_{2}$, the end-
point behaviours of the twist-two distribution amplitudes for the final-state
mesons Braun and Filyanov (1990); Ball  _et al._ (1998), and the momentum-
fraction scaling rules of the Dirac trace in the second and third lines of
(5).

We can readily verify that the obtained result (6) applies for both the
generic $\bar{x}\sim{\cal O}(1)$ and the end-point region $\bar{x}\sim{\cal
O}(\Lambda_{\rm QCD}/m_{B})$, which correspond to the anti-collinear and anti-
soft-collinear Becher  _et al._ (2004, 2016) momentum $q_{2}$. This
observation is in analogy to the smooth interpolation of the $A$-type
contribution to $\bar{B}_{q}\to\gamma\ell\bar{\ell}$ described by the ${\rm
SCET_{I}}$ factorization formula (3.14) in Beneke  _et al._ (2020). In the
anti-collinear $q_{2}$ case, the resulting convolution integrals in (5) become
$\int_{0}^{\infty}d\omega\,\phi_{B}^{+}(\omega)\,\int_{0}^{1}dx\,\int_{0}^{1}dy\,{\phi_{M_{2}}(x)\,\phi_{M_{1}}(y)/(\bar{x}^{2}\,y)}$
by dropping out the subleading term $\omega/m_{B}$ in (6), which reproduce
exactly the obtained expression in Beneke  _et al._ (2001a) and develop the
rapidity divergence thus spoiling the soft-collinear factorization at leading
twist. The appearance of such singularity as $\bar{x}\to 0$ implies that both
the hard and hard-collinear gluon exchanges (labelled with the red colour in
Figure 1) will bring about the leading-power contribution to the flavour
amplitude ${\cal B}_{1}$, in spite of the phase-space suppression for the size
of the hard-collinear region, which is compensated by the QCD enhancement from
the hard-collinear quark and gluon propagators.

Having in mind the essential importance of including the hard and hard-
collinear gluon exchange effects simultaneously, we proceed to write down the
convolution integrals entering the factorized matrix element $\langle
P_{2}^{u}\rangle^{(0)}_{(a)}$

$\displaystyle\hskip 17.07182pt{\cal G}_{{\cal B}_{1}}$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}d\omega\,\phi_{B}^{+}(\omega)\int_{0}^{1}dx\,\phi_{M_{2}}(x)\int_{0}^{1}dy\,\phi_{M_{1}}(y)$
(7)
$\displaystyle\frac{1}{\bar{x}\,y\,\left(\bar{x}-{\omega/m_{B}}+i\epsilon\right)}\,,$

which can be reduced to a rather compact form ${\cal G}_{{\cal B}_{1}}\approx
18\,\left[\left(\ln{m_{B}\over\lambda_{B}}+\gamma_{\rm
E}-2\right)-i\,\pi\right]$ in the leading-power approximation by employing the
Grozin-Neubert model Grozin and Neubert (1997) for $\phi_{B}^{+}(\omega)$ and
the asymptotic twist-two distribution amplitudes $\phi_{M_{2}}(x)$ and
$\phi_{M_{1}}(y)$. The dimensionful quantity $\lambda_{B}^{-1}$ represents the
inverse moment of the above HQET distribution amplitude Grozin and Neubert
(1997); Korchemsky  _et al._ (2000); Beneke and Feldmann (2001) and serves as
an indispensable ingredient for the theory description of a wide variety of
exclusive bottom-meson decays Korchemsky  _et al._ (2000); Descotes-Genon and
Sachrajda (2003); Bosch  _et al._ (2003); Beneke and Rohrwild (2011); Wang
(2016); Wang and Shen (2018); Beneke  _et al._ (2018); Shen  _et al._ (2020);
Beneke  _et al._ (2020, 2021); Wang  _et al._ (2021). It is worthwhile to
stress that we do not aim at implementing QCD resummation of the
perturbatively generated logarithms of ${m_{B}/\lambda_{B}}$, which are
nevertheless not numerically significant for the realistic $\bar{B}_{q}$-meson
mass. The resulting factorization formula (7) provides a dynamical
interpretation of the non-perturbative object $X_{A}^{i}$ introduced in Beneke
_et al._ (2001a); Beneke and Neubert (2003) (with the superscript “i”
characterizing the gluon emission from the initial-state quarks) with
$X_{A}^{i}\approx\left[1+\frac{\left(\gamma_{\rm
E}-1\right)-i\,\pi}{\ln{m_{B}\over\lambda_{B}}}\right]\,\ln{m_{B}\over\lambda_{B}}$
explicitly, which can be matched onto the very parametrization
$X_{A,\,\rm{BBNS}}^{i}=(1+\varrho_{A}\,e^{i\,\varphi_{A}})\,\ln{m_{B}\over\Lambda_{h}}$
suggested in Beneke  _et al._ (2001a); Beneke and Neubert (2003) by setting
$\Lambda_{h}=\lambda_{B}$ and
$(\varrho_{A},\varphi_{A})=\left\\{(0.97,\,-97\degree)\right.$,
$(1.17,\,-97\degree)$, $\left.(1.34,\,-97\degree)\right\\}$ for
$\lambda_{B}=\\{200,\,350,\,500\\}\,{\rm MeV}$, respectively.

We are now ready to present the factorized expressions for the tree-level
flavour amplitudes

$\displaystyle{\cal B}_{1}$ $\displaystyle=$ $\displaystyle{1\over
4}\,\left(C_{2}-{C_{1}\over 2\,N_{c}}\right)\,\widehat{\cal G}_{{\cal
B}_{1}},\,\,{\cal B}_{4}={1\over
8}\,\left[\left(C_{4}+16\,C_{6}\right)\,\widehat{\cal G}_{{\cal
B}_{1}}+\left(C_{4}+4\,C_{6}\right)\,\,\widehat{\cal G}_{{\cal
B}_{4}}\right],\,\,{\cal B}_{4,{\rm EW}}={\cal B}_{4}\left(C_{4}\to
C_{4Q},\,C_{6}\to C_{6Q}\right),\hskip 17.07182pt$ (8)

where $\widehat{\cal G}_{{\cal B}_{1}}={\cal G}_{{\cal
B}_{1}}+\int_{0}^{1}dx\,\int_{0}^{1}dy\,{\phi_{M_{2}}(x)\,\phi_{M_{1}}(y)/(y\,(1-x\,\bar{y}))}$
and $\widehat{\cal G}_{{\cal B}_{4}}$ can be obtained from $\widehat{\cal
G}_{{\cal B}_{1}}$ by performing the replacement $x\leftrightarrow\bar{y}$ for
the counterpart short-distance matching coefficient. Along the same vein we
can derive the relevant results for the NLO topological amplitudes

$\displaystyle\sum_{i=1}^{6}C_{i}\,{\cal P}_{i}^{p,(1)}$ $\displaystyle=$
$\displaystyle\left(C_{2}-{C_{1}\over 2\,N_{c}}\right)\,H_{\rm
1}(m_{p})+\left[\left(C_{3}+16\,C_{5}\right)-{1\over
2\,N_{c}}\,(C_{4}+16\,C_{6})\right]\,\left[H_{1}(m_{b})+H_{1}(0)\right]$
$\displaystyle\hskip
5.69046pt+\left(C_{4}+10\,C_{6}\right)\,\left[H_{1}(m_{b})+H_{1}(m_{c})+3\,H_{1}(0)\right]-\left[5\,C_{4}-8\,C_{5}+4\,\left({1\over
N_{c}}+5\right)C_{6}\right]H_{2}\,,$ $\displaystyle C_{8}^{\rm eff}\,{\cal
P}_{8g}^{(1)}$ $\displaystyle=$ $\displaystyle C_{8}^{\rm
eff}\,H_{3}\,,\qquad\sum_{i=1}^{6}C_{i}\,{\cal
I}_{i}^{p,(1)}=\left[\left(C_{3}+4\,C_{5}\right)+C_{F}\,(C_{4}+4\,C_{6})\right]\,H_{4}\,,$
(9)

where the explicit expressions of the $H_{i}$ functions are presented in the
Supplemental Material. In contrast to QCD factorization for the penguin
contributions $H_{1,2,3}$, the one-loop building blocks ${\cal I}_{i}^{p,(1)}$
from the triangle diagram in Figure 1 (g) (plus the one with the two virtual
gluons exchanged) are insensitive to the (anti)-hard-collinear dynamics at
leading power, thus ensuring the disappearance of the HQET $B$-meson
distribution amplitudes in the factorization formula for $H_{4}$. Adopting the
ansätz for the twist-two light-meson distribution amplitudes
$\phi_{M_{1}}(x)=\phi_{M_{2}}(x)$, the peculiar linear combinations of Wilson
coefficients $C_{i}$ in the yielding result (9) for
$\sum\limits_{i=1}^{6}C_{i}\,{\cal I}_{i}^{p,(1)}$ coincide with the emergent
patterns for the annihilation-type contributions to the axial-vector form
factors of $\bar{B}_{q}\to\gamma\ell\bar{\ell},\,\gamma\gamma$ Beneke  _et
al._ (2020); Shen  _et al._ (2020). This interesting observation, on the one
hand, stems from the Bose-Einstein statistics for the transition matrix
element

$\displaystyle\langle
g^{\ast}(p_{g},\alpha)\,g^{\ast}(\tilde{p}_{g},\beta)|{\cal H}_{\rm
eff}|\bar{B}_{q}\rangle$ $\displaystyle=i\,\epsilon_{\alpha\beta
p_{g}\tilde{p}_{g}}\,F_{V}(p_{g}^{2},\,\tilde{p}_{g}^{2})+g_{\alpha\beta}^{\perp}\,\,F_{A}(p_{g}^{2},\,\tilde{p}_{g}^{2})$
(10)

with the two transversely polarized gluons, which leads to the form-factor
relations
$F_{V}(p_{g}^{2},\,\tilde{p}_{g}^{2})=-F_{V}(\tilde{p}_{g}^{2},\,p_{g}^{2})$
and
$F_{A}(p_{g}^{2},\,\tilde{p}_{g}^{2})=F_{A}(\tilde{p}_{g}^{2},\,p_{g}^{2})$.
On the other hand, this can be attributed to the symmetry constraints of the
QCD matrix element

$\displaystyle\langle
M_{1}(p)\,M_{2}(q)|g^{\ast}(p_{g},\alpha)\,g^{\ast}(\tilde{p}_{g},\beta)\rangle$
$\displaystyle=\left\\{\begin{array}[]{l}g_{\alpha\beta}^{\perp}\,\,{\cal
S}_{\|}(M_{1}M_{2})\hskip 34.14322pt{\rm
for}\,\,M_{1}M_{2}=PP,\,V_{L}V_{L}\,,\\\ i\,\epsilon_{\alpha\beta pq}\,\,{\cal
S}_{\perp}(M_{1}M_{2})\qquad{\rm
for}\,\,M_{1}M_{2}=PV,\,VP\,,\end{array}\right.\hskip 17.07182pt$ (13)

which guarantee the transversity amplitudes ${\cal S}_{\|,\perp}$ proportional
to the products of the twist-two distribution amplitudes
$\phi_{M_{2}}(x)\,\phi_{M_{1}}(y)$. Combining together the requirements of the
hadronic matrix elements (10) and (13) as discussed above, we are then led to
conclude that only the axial-vector form factor $F_{A}$ will be in demand when
evaluating the NLO quantities displayed in (4). It is also apparent that the
observed perturbative enhancement due to the penguin contractions of
$P_{1,2}^{c}$ does not apply to the charmless $\bar{B}_{q}\to PV$ decays,
which will therefore not be discussed in the subsequent numerical analysis.

##  IV Phenomenological implications

We are now equipped to explore the phenomenological significance of the
factorized expressions for the hadronic quantities ${\cal T}^{p,\,(0)}$ and
${\cal T}^{p,\,(1)}$ dictating the exclusive $\bar{B}_{q}\to M_{1}M_{2}$ decay
amplitude (2). To this end, we will employ the three-parameter ansätz of the
leading-twist bottom-meson distribution amplitude as proposed in Beneke  _et
al._ (2018) (see Kawamura  _et al._ (2001, 2003); Huang  _et al._ (2006); Wu
and Huang (2014); Khodjamirian  _et al._ (2010); Wang  _et al._ (2017); Braun
_et al._ (2017); Lü  _et al._ (2019); Gao  _et al._ (2020, 2021) for
additional discussions on the higher-twist distribution amplitudes)
$\phi_{B}^{+}(\omega,\,\mu_{0})={\Gamma(\beta)\over\Gamma(\alpha)}\,U\left(\beta-\alpha,3-\alpha,{\omega\over\omega_{0}}\right){\omega\over\omega_{0}^{2}}{\rm
exp}\left(-{\omega\over\omega_{0}}\right)$ at the reference scale
$\mu_{0}=1\,{\rm GeV}$, where $U(a,b,z)$ represents the confluent
hypergeometric function of the second kind and the allowed intervals of the
shape parameters $\omega_{0}$, $\alpha$ and $\beta$ are determined by
reproducing the numerical values of the inverse-logarithmic moments Beneke
_et al._ (2020); Shen  _et al._ (2020) $\lambda_{B_{d}}=350\pm 150\,{\rm
MeV}$, $\lambda_{B_{s}}=400\pm 150\,{\rm MeV}$,
$\widehat{\sigma}_{B_{d,s}}^{(1)}=0.0\pm 0.7$,
$\widehat{\sigma}_{B_{d,s}}^{(2)}=0.0\pm 6.0$. We proceed to adopt the lattice
predictions for the pseudoscalar-meson decay constants ($f_{B_{d}}$,
$f_{B_{s}}$, $f_{\pi}$, $f_{K}$) as summarized in Aoki  _et al._ (2021) and
take the improved extractions of the decay constants of the longitudinally
polarized vector mesons ($f_{\rho}$, $f_{\omega}$, $f_{\phi}$) from Bharucha
_et al._ (2016). Furthermore, we will truncate the Gegenbauer expansions of
the twist-two light-meson distribution amplitudes at the next-to-next-to-
leading conformal spin accuracy and take advantage of the non-perturbative
determinations of the two lowest moments $a_{1,2}(\mu_{0})$ from Bali  _et
al._ (2019) for the energetic pseudoscalar mesons and from Bharucha  _et al._
(2016) for the light vector mesons. The default values and uncertainties for
the remaining SM parameters follow Table 1 of Shen  _et al._ (2020).

|  $\nu=m_{b}/2$ |  $\nu=m_{b}$ |  $\nu=2\,m_{b}$
---|---|---|---
${\cal T}^{u,\,(0)}$ | $2.03-17.7\,i$ | $2.00-17.4\,i$ | $1.97-17.2\,i$
${\cal T}^{c,\,(0)}$ | $-0.38+3.35\,i$ | $-0.25+2.16\,i$ | $-0.16+1.42\,i$
${\alpha_{s}/(4\,\pi)}\,\,{\cal T}^{u,\,(1)}$ | $-2.29-3.63\,i$ | $-2.41-3.61\,i$ | $-2.47-3.57\,i$
$\,\,{\alpha_{s}/(4\,\pi)}\,\,{\cal T}^{c,\,(1)}$ | $-1.65-2.34\,i$ | $-1.82-2.40\,i$ | $-1.91-2.43\,i$

Table 1: Numerical predictions of the dynamical quantities ${\cal
T}^{p,\,(0)}$ and ${\cal T}^{p,\,(1)}$ ($p=u,\,c$) for
$\bar{B}_{s}\to\pi^{+}\pi^{-}$ with three distinct renormalization scales
$\nu$.

In order to develop a transparent understanding of the numerical feature for
the identified enhancement mechanism, we present in Table 1 the obtained
results of ${\cal T}^{p,\,(0)}$ and ${\cal T}^{p,\,(1)}$ for the weak
annihilation decay $\bar{B}_{s}\to\pi^{+}\pi^{-}$ with three distinct values
of the renormalization scale $\nu=m_{b}/2,\,m_{b},\,2\,m_{b}$, while employing
the central values of the remaining theory inputs. Interestingly, the dominant
NLO contribution from the charm-loop diagrams displayed in Figure 1 (c) and
(d) (plus the two diagrams due to the crossing symmetry) will result in the
significant increase of the real part of ${\cal T}^{c,\,(0)}$ in magnitude and
simultaneously generates the considerable cancellation of its imaginary part,
keeping in mind the CKM hierarchy
$V_{ub}V_{us}^{\ast}/V_{cb}V_{cs}^{\ast}=0.021\,e^{-1.20\,i}$. In consequence,
the perturbative rescattering effect under discussion provides an important
source of the strong phases of the weak-annihilation non-leptonic
$\bar{B}_{q}$-meson decay amplitudes, which can then be probed feasibly by
investigating the CP violating observables.

|  ${\cal A}_{\rm CP}^{\rm dir}$ |  ${\cal A}_{\rm CP}^{\rm mix}$
---|---|---
$\bar{B}_{s}\to\pi^{+}\,\pi^{-},\,\pi^{0}\,\pi^{0}$ |  $-36.3^{+8.2}_{-1.3}$ ($0.0\pm 0.0$) |  $-4.2^{+21.4}_{-9.0}$ ($35.9^{+15.6}_{-11.2}$)
$\bar{B}_{s}\to\rho^{+}_{L}\,\rho^{-}_{L},\,\rho^{0}_{L}\,\rho^{0}_{L}$ |  $-36.3^{+8.3}_{-1.8}$ ($0.0\pm 0.0$) |  $-4.3^{+21.5}_{-9.0}$ ($35.9^{+15.6}_{-11.2}$)
$\bar{B}_{s}\to\omega_{L}\,\omega_{L}$ |  $-36.3^{+8.3}_{-3.1}$ ($0.0\pm 0.0$) |  $-3.8^{+21.8}_{-9.7}$ ($35.9^{+15.6}_{-11.2}$)
$\bar{B}_{s}\to\rho_{L}\,\omega_{L}$ |  $0.0\pm 0.0$ ($0.0\pm 0.0$ ) |  $-71.0^{+6.3}_{-5.4}$ ($-71.0^{+6.3}_{-5.4}$)
$\bar{B}_{d}\to K^{+}\,K^{-}$ |  $39.0^{+3.2}_{-5.6}$ ($0.0\pm 0.0$) |  $-2.2^{+19.1}_{-26.4}$ ($-47.0^{+15.7}_{-18.8}$)
$\bar{B}_{d}\to K^{\ast+}_{L}\,K^{\ast-}_{L}$ |  $39.6^{+4.9}_{-6.7}$ ($0.0\pm 0.0$) |  $-1.4^{+19.7}_{-26.9}$ ($-47.0^{+15.7}_{-18.8}$)
$\bar{B}_{d}\to\phi_{L}\,\phi_{L}$ |  $38.3^{+11.4}_{-15.8}$ ($0.0\pm 0.0$) |  $27.8^{+5.7}_{-25.9}$ ($0.0\pm 0.0$)

Table 2: Theory predictions of the CP asymmetries ${\cal A}_{\rm CP}^{\rm
dir}$ and ${\cal A}_{\rm CP}^{\rm mix}$ (in unites of $10^{-2}$) for the weak
annihilation $\bar{B}_{q}\to PP,\,V_{L}V_{L}$ decay processes, where we have
included in parentheses the corresponding LO QCD results for a comparison.

It is straightforward to derive the time-dependent CP asymmetries for the
neutral $\bar{B}_{q}$-meson decaying into CP eigenstates

$\displaystyle{\cal A}_{\rm CP}(t)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\bar{B}_{q}\to M_{1}M_{2})-\Gamma(B_{q}\to
M_{1}M_{2})}{\Gamma(\bar{B}_{q}\to M_{1}M_{2})+\Gamma(B_{q}\to M_{1}M_{2})}$
(14) $\displaystyle=$ $\displaystyle-\frac{{\cal A}_{\rm CP}^{\rm
dir}\,\cos(\Delta m_{q}\,t)+{\cal A}_{\rm CP}^{\rm mix}\,\sin(\Delta
m_{q}\,t)}{\cosh(\Delta\Gamma_{q}\,t/2)+{\cal
A}_{\Delta\Gamma}\,\sinh(\Delta\Gamma_{q}\,t/2)}\,,\hskip 14.22636pt$

where the general definitions of ${\cal A}_{\rm CP}^{\rm dir}$, ${\cal A}_{\rm
CP}^{\rm mix}$ and ${\cal A}_{\Delta\Gamma}$ can be found, for instance, in
Zyla  _et al._ (2020); Nierste (2009); Pro (2001), satisfying the exact
algebra relation $|{\cal A}_{\rm CP}^{\rm dir}|^{2}+|{\cal A}_{\rm CP}^{\rm
mix}|^{2}+|{\cal A}_{\Delta\Gamma}|^{2}=1$. The resulting predictions for the
two independent CP asymmetries displayed in Table 2 evidently reveal the
striking numerical impacts of the enhanced NLO QCD correction to the weak
annihilation amplitudes described in this Letter, without which the direct CP
violation for the entire pure annihilation channels would vanish due to the
absence of at least two separate terms with different strong phases involved
in $\bar{\cal A}(\bar{B}_{q}\to M_{1}M_{2})$. By contrast, the time-dependent
CP asymmetries of $\bar{B}_{s}\to\rho_{L}\,\omega_{L}$ turn out to be
insensitive to this particular perturbative contribution in question, on
account of the isopsin non-conservation of the QCD scattering amplitude
$\langle\rho_{L}(p)\,\omega_{L}(q)|g^{\ast}(p_{g},\alpha)\,g^{\ast}(\tilde{p}_{g},\beta)\rangle$.
We further mention in passing that the newly computed NLO quantities ${\cal
T}^{p,\,(1)}$ generate the subdominant corrections to the CP-averaged
branching fractions of the pure annihilation decay channels (numerically at
the level of ${\cal O}(10\,\%)$).

##  V Conclusions

To summarize, we have endeavored to achieve the analytical regularization of
end-point divergences in the factorization formulae of the peculiar weak
annihilation non-leptonic bottom-meson decay amplitudes, due to the gluon
radiation off the initial-state partons, by adding the missing hard-collinear
contribution on top of the hard gluon exchange effect, both of which led to
the leading power contributions in the heavy quark expansion. Applying the
improved factorization formalism, we then identified the novel perturbative
mechanism generated by the penguin contractions of the current-current
operators, yielding the significant numerical impacts on the strong phases of
the pure annihilation amplitudes. Our results are of importance for enhancing
the predictive power of the perturbative factorization approach in addressing
the entire spectrum of CP violating observables from exclusive heavy hadron
decays, accessible at the LHCb and Belle II experiments.

###### Acknowledgements.

## Acknowledgements

We would like to thank Martin Beneke for interesting discussions and comments
on the manuscript. The research of C.D.L. is supported by the National Natural
Science Foundation of China with Grant No. 11521505 and 12070131001 as well as
the National Key Research and Development Program of China under Contract No.
2020YFA0406400. C.W. is supported in part by the National Natural Science
Foundation of China with Grant No. 12105112 and the Natural Science Foundation
of Jiangsu Education Committee with Grant No. 21KJB140027. The research of
Y.L.S. is supported by the National Natural Science Foundation of China with
Grant No. 12175218 and the Natural Science Foundation of Shandong with Grant
No. ZR2020MA093. Y.M.W. acknowledges support from the National Natural Science
Foundation of China with Grant No. 11735010 and 12075125, and the Natural
Science Foundation of Tianjin with Grant No. 19JCJQJC61100.

##  Appendix A Supplemental Material: Analytic Expressions for the NLO Kernels

We collect QCD factorization formulae for the primitive kernels entering the
NLO topological amplitudes (9)

$\displaystyle H_{i\leq 3}={1\over
12}\,\int_{0}^{\infty}d\omega\,\phi_{B}^{+}(\omega)\int_{0}^{1}dx\,\phi_{M_{2}}(x)\int_{0}^{1}dy\,\phi_{M_{1}}(y)\,\left[\left({h_{i}\over\bar{x}\,y}\,\left(\frac{1}{\bar{x}-{\omega/m_{B}}+i\epsilon}+{\eta_{i}\,\bar{x}\over
1-x\,\bar{y}}\right)+\\{x\leftrightarrow\bar{y}\\}\right)+\\{x\leftrightarrow
y\\}\right],$ $\displaystyle
H_{4}=2\,\int_{0}^{1}dx\,\phi_{M_{2}}(x)\int_{0}^{1}dy\,\phi_{M_{1}}(y)\,{1\over\lambda}\,\left[h_{4}\,{\bar{x}+y\over\bar{x}\,y}+{1\over
x\,\bar{x}}\,\left({\pi\over\sqrt{3}}-{\lambda\over
2\,y\,\bar{y}}\right)+\left(1+{\lambda\over 4\,\bar{x}\,y}\,\left(2+{3\over
x\,\bar{y}}\right)\right)\,\eta_{4}+\\{x\leftrightarrow y\\}\right],$

where for brevity we have introduced the following perturbative functions

$\displaystyle h_{1}$ $\displaystyle=$
$\displaystyle\ln{m_{q}^{2}\over\mu^{2}}+{1\over 3}-{4\,z_{q}\over
x\,\bar{y}}+2\,\left(1+{2\,z_{q}\over x\bar{y}}\right)\,\sqrt{{4\,z_{q}\over
x\bar{y}}-1}\,\arctan{1\over\sqrt{{4\,z_{q}\over x\bar{y}}-1}}\,,\qquad
z_{q}={m_{q}^{2}\over m_{B}^{2}}\,,$ $\displaystyle h_{2}$ $\displaystyle=$
$\displaystyle 1,\qquad h_{3}=-{3\over\bar{y}}\,,\qquad h_{4}=-2\sqrt{{4\over
x\bar{y}}-1}\,\arctan{1\over\sqrt{{4\over
x\bar{y}}-1}}\,,\qquad\eta_{1}=\eta_{2}=1,\qquad\eta_{3}=-{y\over\bar{x}}\,,$
$\displaystyle\eta_{4}$ $\displaystyle=$
$\displaystyle{1\over\sqrt{\lambda}}\,\sum_{i=1}^{3}\left[{\rm
Li}_{2}\left({\beta_{i}-1\over\alpha_{i}+\beta_{i}}\right)+{\rm
Li}_{2}\left(-{\beta_{i}-1\over\alpha_{i}-\beta_{i}}\right)-{\rm
Li}_{2}\left({\beta_{i}+1\over\alpha_{i}+\beta_{i}}\right)-{\rm
Li}_{2}\left(-{\beta_{i}+1\over\alpha_{i}-\beta_{i}}\right)\right]\,,$ (16)

with the Källén function
$\lambda\equiv\lambda(1,\bar{x}\,y,x\,\bar{y})=(1-x-y)^{2}$ and

$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle\sqrt{1+{4\over-\bar{x}\,y-i\epsilon}}\,,\qquad\alpha_{2}=\sqrt{1+{4\over-x\,\bar{y}-i\epsilon}}\,,\qquad\alpha_{3}=\sqrt{5}\,,$
$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle{-\bar{x}\,y+x\,\bar{y}+1\over\sqrt{\lambda}}\,,\qquad\beta_{2}={-x\,\bar{y}+1+\bar{x}\,y\over\sqrt{\lambda}}\,,\qquad\beta_{3}={-1+x\,\bar{y}+\bar{x}\,y\over\sqrt{\lambda}}\,.$
(17)

## References

  * Beneke  _et al._ (1999) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, Phys. Rev. Lett. 83, 1914 (1999), arXiv:hep-ph/9905312 .
  * Beneke  _et al._ (2000) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, Nucl. Phys. B 591, 313 (2000), arXiv:hep-ph/0006124 .
  * Chay and Kim (2004) J. Chay and C. Kim, Nucl. Phys. B 680, 302 (2004), arXiv:hep-ph/0301262 .
  * Beneke  _et al._ (2002) M. Beneke, A. P. Chapovsky, M. Diehl, and T. Feldmann, Nucl. Phys. B 643, 431 (2002), arXiv:hep-ph/0206152 .
  * Bauer  _et al._ (2004) C. W. Bauer, D. Pirjol, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D 70, 054015 (2004), arXiv:hep-ph/0401188 .
  * Bauer  _et al._ (2006) C. W. Bauer, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D 74, 034010 (2006), arXiv:hep-ph/0510241 .
  * Beneke  _et al._ (2001a) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, Nucl. Phys. B 606, 245 (2001a), arXiv:hep-ph/0104110 .
  * Beneke and Neubert (2003) M. Beneke and M. Neubert, Nucl. Phys. B 675, 333 (2003), arXiv:hep-ph/0308039 .
  * Beneke  _et al._ (2007) M. Beneke, J. Rohrer, and D. Yang, Nucl. Phys. B 774, 64 (2007), arXiv:hep-ph/0612290 .
  * Bartsch  _et al._ (2008) M. Bartsch, G. Buchalla, and C. Kraus, (2008), arXiv:0810.0249 [hep-ph] .
  * Cheng and Yang (2008) H.-Y. Cheng and K.-C. Yang, Phys. Rev. D 78, 094001 (2008), [Erratum: Phys.Rev.D 79, 039903 (2009)], arXiv:0805.0329 [hep-ph] .
  * Cheng and Chua (2009a) H.-Y. Cheng and C.-K. Chua, Phys. Rev. D 80, 114008 (2009a), arXiv:0909.5229 [hep-ph] .
  * Cheng and Chua (2009b) H.-Y. Cheng and C.-K. Chua, Phys. Rev. D 80, 114026 (2009b), arXiv:0910.5237 [hep-ph] .
  * Zhu (2011) G. Zhu, Phys. Lett. B 702, 408 (2011), arXiv:1106.4709 [hep-ph] .
  * Wang and Zhu (2013) K. Wang and G. Zhu, Phys. Rev. D 88, 014043 (2013), arXiv:1304.7438 [hep-ph] .
  * Bobeth  _et al._ (2015) C. Bobeth, M. Gorbahn, and S. Vickers, Eur. Phys. J. C 75, 340 (2015), arXiv:1409.3252 [hep-ph] .
  * Chang  _et al._ (2015) Q. Chang, J. Sun, Y. Yang, and X. Li, Phys. Lett. B 740, 56 (2015), arXiv:1409.2995 [hep-ph] .
  * Chang  _et al._ (2016) Q. Chang, X.-N. Li, J.-F. Sun, and Y.-L. Yang, J. Phys. G 43, 105004 (2016), arXiv:1610.02747 [hep-ph] .
  * Chang  _et al._ (2017) Q. Chang, X. Li, X.-Q. Li, and J. Sun, Eur. Phys. J. C 77, 415 (2017), arXiv:1706.06138 [hep-ph] .
  * Arnesen  _et al._ (2008) C. M. Arnesen, Z. Ligeti, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D 77, 054006 (2008), arXiv:hep-ph/0607001 .
  * Khodjamirian  _et al._ (2005) A. Khodjamirian, T. Mannel, M. Melcher, and B. Melic, Phys. Rev. D 72, 094012 (2005), arXiv:hep-ph/0509049 .
  * Feldmann and Hurth (2004) T. Feldmann and T. Hurth, JHEP 11, 037 (2004), arXiv:hep-ph/0408188 .
  * Keum  _et al._ (2001) Y. Y. Keum, H.-N. Li, and A. I. Sanda, Phys. Rev. D 63, 054008 (2001), arXiv:hep-ph/0004173 .
  * Lu  _et al._ (2001) C.-D. Lu, K. Ukai, and M.-Z. Yang, Phys. Rev. D 63, 074009 (2001), arXiv:hep-ph/0004213 .
  * Li  _et al._ (2004) Y. Li, C.-D. Lu, Z.-J. Xiao, and X.-Q. Yu, Phys. Rev. D 70, 034009 (2004), arXiv:hep-ph/0404028 .
  * Ali  _et al._ (2007) A. Ali, G. Kramer, Y. Li, C.-D. Lu, Y.-L. Shen, W. Wang, and Y.-M. Wang, Phys. Rev. D 76, 074018 (2007), arXiv:hep-ph/0703162 .
  * Li  _et al._ (2011) H.-n. Li, Y.-L. Shen, Y.-M. Wang, and H. Zou, Phys. Rev. D 83, 054029 (2011), arXiv:1012.4098 [hep-ph] .
  * Li  _et al._ (2012) H.-n. Li, Y.-L. Shen, and Y.-M. Wang, Phys. Rev. D 85, 074004 (2012), arXiv:1201.5066 [hep-ph] .
  * Li  _et al._ (2013) H.-N. Li, Y.-L. Shen, and Y.-M. Wang, JHEP 02, 008 (2013), arXiv:1210.2978 [hep-ph] .
  * Li  _et al._ (2014) H.-N. Li, Y.-L. Shen, and Y.-M. Wang, JHEP 01, 004 (2014), arXiv:1310.3672 [hep-ph] .
  * Li and Wang (2015) H.-n. Li and Y.-M. Wang, JHEP 06, 013 (2015), arXiv:1410.7274 [hep-ph] .
  * Altmannshofer  _et al._ (2019) W. Altmannshofer _et al._ (Belle-II), PTEP 2019, 123C01 (2019), [Erratum: PTEP 2020, 029201 (2020)], arXiv:1808.10567 [hep-ex] .
  * Chetyrkin  _et al._ (1997) K. G. Chetyrkin, M. Misiak, and M. Munz, Phys. Lett. B 400, 206 (1997), [Erratum: Phys.Lett.B 425, 414 (1998)], arXiv:hep-ph/9612313 .
  * Chetyrkin  _et al._ (1998) K. G. Chetyrkin, M. Misiak, and M. Munz, Nucl. Phys. B 520, 279 (1998), arXiv:hep-ph/9711280 .
  * Bobeth  _et al._ (2000) C. Bobeth, M. Misiak, and J. Urban, Nucl. Phys. B 574, 291 (2000), arXiv:hep-ph/9910220 .
  * Bobeth  _et al._ (2004) C. Bobeth, P. Gambino, M. Gorbahn, and U. Haisch, JHEP 04, 071 (2004), arXiv:hep-ph/0312090 .
  * Huber  _et al._ (2006) T. Huber, E. Lunghi, M. Misiak, and D. Wyler, Nucl. Phys. B 740, 105 (2006), arXiv:hep-ph/0512066 .
  * Beneke  _et al._ (2001b) M. Beneke, T. Feldmann, and D. Seidel, Nucl. Phys. B 612, 25 (2001b), arXiv:hep-ph/0106067 .
  * Chernyak and Zhitnitsky (1982) V. L. Chernyak and A. R. Zhitnitsky, Nucl. Phys. B 201, 492 (1982), [Erratum: Nucl.Phys.B 214, 547 (1983)].
  * Beneke and Feldmann (2001) M. Beneke and T. Feldmann, Nucl. Phys. B 592, 3 (2001), arXiv:hep-ph/0008255 .
  * Beneke (2002) M. Beneke, Nucl. Phys. B Proc. Suppl. 111, 62 (2002), arXiv:hep-ph/0202056 .
  * Grozin and Neubert (1997) A. G. Grozin and M. Neubert, Phys. Rev. D 55, 272 (1997), arXiv:hep-ph/9607366 .
  * Bell  _et al._ (2013) G. Bell, T. Feldmann, Y.-M. Wang, and M. W. Y. Yip, JHEP 11, 191 (2013), arXiv:1308.6114 [hep-ph] .
  * Braun and Filyanov (1990) V. M. Braun and I. E. Filyanov, Z. Phys. C 48, 239 (1990).
  * Ball  _et al._ (1998) P. Ball, V. M. Braun, Y. Koike, and K. Tanaka, Nucl. Phys. B 529, 323 (1998), arXiv:hep-ph/9802299 .
  * Becher  _et al._ (2004) T. Becher, R. J. Hill, and M. Neubert, Phys. Rev. D 69, 054017 (2004), arXiv:hep-ph/0308122 .
  * Becher  _et al._ (2016) T. Becher, M. Neubert, L. Rothen, and D. Y. Shao, Phys. Rev. Lett. 116, 192001 (2016), arXiv:1508.06645 [hep-ph] .
  * Beneke  _et al._ (2020) M. Beneke, C. Bobeth, and Y.-M. Wang, JHEP 12, 148 (2020), arXiv:2008.12494 [hep-ph] .
  * Korchemsky  _et al._ (2000) G. P. Korchemsky, D. Pirjol, and T.-M. Yan, Phys. Rev. D 61, 114510 (2000), arXiv:hep-ph/9911427 .
  * Descotes-Genon and Sachrajda (2003) S. Descotes-Genon and C. T. Sachrajda, Nucl. Phys. B 650, 356 (2003), arXiv:hep-ph/0209216 .
  * Bosch  _et al._ (2003) S. W. Bosch, R. J. Hill, B. O. Lange, and M. Neubert, Phys. Rev. D 67, 094014 (2003), arXiv:hep-ph/0301123 .
  * Beneke and Rohrwild (2011) M. Beneke and J. Rohrwild, Eur. Phys. J. C 71, 1818 (2011), arXiv:1110.3228 [hep-ph] .
  * Wang (2016) Y.-M. Wang, JHEP 09, 159 (2016), arXiv:1606.03080 [hep-ph] .
  * Wang and Shen (2018) Y.-M. Wang and Y.-L. Shen, JHEP 05, 184 (2018), arXiv:1803.06667 [hep-ph] .
  * Beneke  _et al._ (2018) M. Beneke, V. M. Braun, Y. Ji, and Y.-B. Wei, JHEP 07, 154 (2018), arXiv:1804.04962 [hep-ph] .
  * Shen  _et al._ (2020) Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, JHEP 12, 169 (2020), arXiv:2009.02723 [hep-ph] .
  * Beneke  _et al._ (2021) M. Beneke, P. Böer, P. Rigatos, and K. K. Vos, Eur. Phys. J. C 81, 638 (2021), arXiv:2102.10060 [hep-ph] .
  * Wang  _et al._ (2021) C. Wang, Y.-M. Wang, and Y.-B. Wei, (2021), arXiv:2111.11811 [hep-ph] .
  * Kawamura  _et al._ (2001) H. Kawamura, J. Kodaira, C.-F. Qiao, and K. Tanaka, Phys. Lett. B 523, 111 (2001), [Erratum: Phys.Lett.B 536, 344–344 (2002)], arXiv:hep-ph/0109181 .
  * Kawamura  _et al._ (2003) H. Kawamura, J. Kodaira, C.-F. Qiao, and K. Tanaka, Int. J. Mod. Phys. A 18, 1433 (2003), arXiv:hep-ph/0112146 .
  * Huang  _et al._ (2006) T. Huang, C.-F. Qiao, and X.-G. Wu, Phys. Rev. D 73, 074004 (2006), arXiv:hep-ph/0507270 .
  * Wu and Huang (2014) X.-G. Wu and T. Huang, Chin. Sci. Bull. 59, 3801 (2014), arXiv:1312.1455 [hep-ph] .
  * Khodjamirian  _et al._ (2010) A. Khodjamirian, T. Mannel, A. A. Pivovarov, and Y. M. Wang, JHEP 09, 089 (2010), arXiv:1006.4945 [hep-ph] .
  * Wang  _et al._ (2017) Y.-M. Wang, Y.-B. Wei, Y.-L. Shen, and C.-D. Lü, JHEP 06, 062 (2017), arXiv:1701.06810 [hep-ph] .
  * Braun  _et al._ (2017) V. M. Braun, Y. Ji, and A. N. Manashov, JHEP 05, 022 (2017), arXiv:1703.02446 [hep-ph] .
  * Lü  _et al._ (2019) C.-D. Lü, Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, JHEP 01, 024 (2019), arXiv:1810.00819 [hep-ph] .
  * Gao  _et al._ (2020) J. Gao, C.-D. Lü, Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, Phys. Rev. D 101, 074035 (2020), arXiv:1907.11092 [hep-ph] .
  * Gao  _et al._ (2021) J. Gao, T. Huber, Y. Ji, C. Wang, Y.-M. Wang, and Y.-B. Wei, (2021), arXiv:2112.12674 [hep-ph] .
  * Aoki  _et al._ (2021) Y. Aoki _et al._ , (2021), arXiv:2111.09849 [hep-lat] .
  * Bharucha  _et al._ (2016) A. Bharucha, D. M. Straub, and R. Zwicky, JHEP 08, 098 (2016), arXiv:1503.05534 [hep-ph] .
  * Bali  _et al._ (2019) G. S. Bali, V. M. Braun, S. Bürger, M. Göckeler, M. Gruber, F. Hutzler, P. Korcyl, A. Schäfer, A. Sternbeck, and P. Wein (RQCD), JHEP 08, 065 (2019), [Addendum: JHEP 11, 037 (2020)], arXiv:1903.08038 [hep-lat] .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Nierste (2009) U. Nierste, in  _Helmholz International Summer School on Heavy Quark Physics_ (2009) pp. 1–38, arXiv:0904.1869 [hep-ph] .
  * Pro (2001) _$B$ physics at the Tevatron: Run II and beyond_ (2001) arXiv:hep-ph/0201071 .