# Single-photon generation at room temperature using molecular optomechanics
in a hybrid photonic-plasmonic cavity

Shabnam Abutalebi B.A.1  Seyed Mahmoud Ashrafi2,3  Hassan RanjbarAskari1
Alireza Bahrampour2,3  1 Department of Physics, Faculty of Science, Vali-e-Asr
University of Rafsanjan, Rafsanjan, Iran
2 Department of Physics, Sharif University of Technology, Tehran, Iran
3 Center for Quantum Engineering and Photonics Technology, Sharif University
of Technology, Tehran, Iran<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

We propose a novel integrated structure for single photon generation at room
temperature based on a molecular optomechanics system in a hybrid photonic-
plasmonic cavity. The proposed structure comprises a single molecule within a
plasmonic cavity, coupled to a 2D photonic crystal resonator. In this paper,
we theoretically identify the ability of the scheme through calculation second
order correlation function $g^{2}(0)$ for four different coupling regimes. We
demonstrate the quantum paths and the destructive interference mechanism
through the selection of efficient and preferred basis. Furthermore, we find
that the unconventional photon blockade effects can occurs in the weak
molecular optomechanics coupling. This structure holds the potential to serve
as an integrated single-photon source for quantum networks at room
temperature.

##  1 Introduction

In the field of quantum optics, where the manipulation of individual photons
and their interactions holds the key to quantum information processing and
communication, the concept of photon blockade (PB) emerges as a pivotal
phenomenon [1, 2, 3]. This effect pertains to the strong suppression of
subsequent photon emissions following the initial emission of a single photon,
fundamentally altering the probabilistic nature of photon emission processes
[4, 5]. Second-order correlation function ($g^{(2)}$) is key for single-photon
source quality, indicating multi-photon emission probability and mean-variance
relationship. If variance is less than mean ($g^{(2)}<1$), indicating sub-
Poissonian behavior [6]. Until now, three photon blockade with different
physical mechanisems have been arisen. first, conventional photon blockade
(CPB) arises from anharmonicity, where nonlinearity in the system prevents the
absorption of a second photon with a specific frequency. This mechanism
require strong coupling regime and explored in various systems, often involves
Kerr-type interactions or atom-resonator couplings [7, 4]. Second,
unconventional photon blockade (UCPB) arises due to quantum destructive
interference and excels in situations characterized by extremely low mean
photon numbers [8, 9, 10]. While UCPB may reduce the probability of generating
a single photon, its ability to suppress multi-photon states and induce
higher-order coherence presents both a challenge and an advantage in the
context of single-photon source applications. The last, non-Hermitian photon
blockade (NHPB) mechanism occurs when there is a significant difference in the
losses between the singly and doubly excited states of the system [11]. This
mechanism is not constrained by linewidths in contrast to conventional photon
blockade.
Photon blockade has been predicted in various optomechanical systems,
including Fabry-Perot cavities [12], microtoroids [13] and photonic crystal
cavities [14]. Recently, a new optomechanical approach has been introduced for
surface-enhanced Raman scattering (SERS) in which there is an interesting
optomechanical coupling between molecular vibrations and localized surface
plasmon resonance (LSPR) mode similar to conventional cavity optomechanical
systems [15, 16]. These new approaches could perfectly demonstrate some
unknown observation in experiments like back-action amplification of Raman
enhancement in the molecular samples [17]. In molecular optomechanics,
molecular vibrations serve as some of the smallest mechanical oscillators
amenable to nano engineering. The energy and strength of molecular
oscillations depend on specific functional groups and their chemical and
physical environments [18]. An important advantage of utilizing molecules as
mechanical oscillators in optomechanical systems is their resistance to
temperature effects, a factor that previously constrained the practical
application of optomechanical systems. To enhance the interaction between
molecule and photons, plasmonic cavity with small mode volume has been
employed. However, this cavity often has a low Q-factor. Combining plasmonic
cavity with high Q-factor optical cavity has enabled the compensation of this
limitation, making them valuable for nanoscale optical physics and integrated
photonics applications [19, 20, 21, 22, 23, 24].
In this study, we proposed a hybrid plasmonic- photonics resonator containing
a molecular optomechanics system and a 2D photonic crystal. The 2D photonic
crystal, designed with periodic dielectric structures exhibiting photonic
bandgap properties, provide high-quality resonance. This hybrid configuration,
which confines electromagnetic fields at subwavelength scales through LSPR
cavity, enhances light-matter interaction. On the other hand, plasmonic
resonance supported by metallic nanocavity enable the boosting of Raman
scattering from a single molecule. What sets our system apart is its
operational viability at room temperature, eliminating the need for the
extreme cooling methods often associated with traditional quantum systems.
This characteristic not only simplifies experimental setups but also holds
promise for practical, real-world applications. The remaining part of the
paper is organized as follows: In Sec. 2, we introduce our proposed structure.
The theoretical model of the proposed structure is presented in Sec. 3.
Numerical calculations and results can be found in Sec. 4. The paper concludes
with a summary in Sec. 5.

##  2 Proposed structure

The primary objective of this study is to explore the capacity of a molecular
optomechanical system to generate single photon at room temperature. We
proposed a molecular optomechanical system consists of a Raman active molecule
situated in the hot spot of the bowtie nano-antenna (BNA) coupled to photonic
crystal (PC) resonator. BNA is highly advantageous for confining light within
ultra-small volume, surpassing the limits of diffraction but it has low
quality factor. For compensation of this limitation, the hybrid PC-BNA nano-
resonator employed. The PC has a high quality factor and is pumped by a laser
with frequency $\omega_{l}$ and amplitude $\Omega$. Additionally, there is a
waveguide for driving the cavity and the path of a single photon that is
generated. This system shown in Fig. 1(a), and the interaction diagram of this
system is presented in Fig. 1(b).

Figure 1: (a) Schematic diagram of the photonic crystal cavity coupled to a
plasmonic nano-antenna cavity with a molecule positioned in the hot spot of
the nano-antenna to form a molecular optomechanical system. (b) Interaction
diagram of the cavities and molecule with damping rates.

##  3 Theoretical model

The Hamiltonian for the PC and BNA modes with central frequencies $\omega_{c}$
and $\omega_{p}$ is expressed as
$H_{0}=\sum_{i=c,p}\omega_{i}a_{i}^{\dagger}a_{i}$, where
$a_{c}^{\dagger}(a_{c})$ and $a_{p}^{\dagger}(a_{p})$ represent the bosonic
creation (annihilation) operators for the optical and LSPR modes,
respectively. The system of units is chosen such that $\hbar=1$. The molecule
can also be treated as quantum harmonic oscillator with the frequency
$\omega_{m}$. Consequently, its corresponding Hamiltonian is
$H_{m}=\omega_{m}b^{\dagger}b$, where $b^{\dagger}(b)$ is the bosonic creation
(annihilation) operator for the vibrational mode. It is worth to mention that
the primary assumption for coherent coupling between mechanical oscillator and
LSPR mode is $(\omega_{p}\gg\omega_{m})$ [15]. The corresponding molecular
optomechanics interaction can be written as
$H_{\text{BNA-m}}=ga_{p}^{\dagger}a_{p}(b^{\dagger}+b)$, in which the
optomechanical coupling factor is
$g=\frac{\omega_{p}}{\epsilon_{0}V_{p}}\frac{1}{\sqrt{R_{m}/2\omega_{m}}}$,
with $R_{m}$ representing the Raman activity associated with the vibration
under study and $V_{p}$ being the effective mode volume of the plasmonic
cavity mode [18]. Due to the large mode volume of the PC cavity, we can
neglect the optomechanical interaction between the PC mode and the molecule.
Meanwhile, we consider an interaction between the non-local optical mode in
the PC and LSPR modes by the Hamiltonian $H_{\text{PC-
BNA}}=J(a_{c}^{\dagger}a_{p}+a_{p}^{\dagger}a_{c})$. Here, $J$ is the coupling
strength between the PC and LSPR modes, strongly dependent on the distance
between the nano-antenna and the PC cavity [25]. By performing a rotating
transformation defined as
$U=\exp\left[-i(\omega_{c}a^{\dagger}_{c}a_{c}+\omega_{p}a^{\dagger}_{p}a_{p})t\right],$
the system Hamiltonian is transformed as:

$\displaystyle
H=\Delta_{c}a_{c}^{\dagger}a_{c}+\Delta_{p}a_{p}^{\dagger}a_{p}+\omega_{m}b^{\dagger}b+J(a_{c}^{\dagger}a_{p}+a_{p}^{\dagger}a_{c})-ga_{p}^{\dagger}a_{p}(b^{\dagger}+b)+\Omega(a_{c}^{\dagger}+a_{c}).$
(1)

Here, $\Delta_{c}=\omega_{c}-\omega_{l}$ and
$\Delta_{p}=\omega_{p}-\omega_{l}$ represent the detuning of the photonic and
plasmonic cavities from the pump frequency, respectively. The first three
terms correspond to the optical cavity, plasmon cavity, and phonons of
molecular vibration. The eigenvalues of these terms are harmonic, and their
eigenvectors are denoted as $|{nmq}\rangle$ (bare basis), where n, m, and q
represent the photon number in the optical cavity, the plasmon number, and the
phonon number of molecular vibration, respectively. The fourth term accounts
for the photon-plasmon cavity coupling. Since this term is bilinear with
respect to annihilation and creation operators, its eigenvalues are also
harmonic and do not induce the CPB effect. The corresponding eigenvectors are
denoted as $|n_{+}n_{-}q\rangle$, where $n_{+}$ and $n_{-}$ represent the
number of particles associated with resonance frequencies $\omega_{+}$ and
$\omega_{-}$, respectively. The $\omega_{+}$ and $\omega_{-}$ are introduced
in Sec. 3.2. The fifth term is non-harmonic and represents the interaction
between plasmons and phonons. This non-harmonic term results in non-harmonic
eigenvalues for the first five terms of the Hamiltonian. The eigenvectors for
the first five terms of the Hamiltonian are denoted as $|E_{nmq}\rangle$. This
term introduces nonlinear contributions in the Heisenberg-Langevin equation
and leads to the generation of new frequencies, known as the Stokes and anti-
Stokes frequencies of different orders. These frequencies can be used for
cooling and warming of cavity modes. The last term represents the laser input,
and its effect is to displace the eigenvectors and eigenvalues.
In this work, we analyze the statistical properties of photons using the
second-order correlation function at zero delay time.

$\displaystyle g^{(2)}(0)=\frac{{\langle
a^{\dagger}a^{\dagger}aa\rangle}}{{\langle a^{\dagger}a\rangle^{2}}}$ (2)

$a$ and $a^{\dagger}$ are the annihilation and creation operators of the mode
under investigation. The value of $g^{(2)}(0)<1$ corresponds to the single
photon statistics. We utilize both the Lindblad equation and the Schrödinger
equation with non-Hermitian Hamiltonian methods. To ensure the validity of our
findings, we compare the results obtained from these two approaches. To
comprehensively characterize the dynamical behavior of the system, it is
crucial to consider the interactions between phonons, plasmons, photons modes,
and their respective environments.

###  3.1 The Lindblad equation

The decoherence of the system can be described by the evolution of the density
matrix using the Lindblad master equation [26]:

$\displaystyle\frac{{d\rho}}{{dt}}=-i[H,\rho]+\mathcal{L}_{a_{c}}+\mathcal{L}_{a_{p}}+\mathcal{L}_{b}$
$\displaystyle\mathcal{L}_{a_{i}}=\frac{\kappa_{i}}{2}\mathcal{D}_{a_{i}}[\rho],\quad
i=c,p$ (3)
$\displaystyle\mathcal{L}_{b}=\frac{\gamma_{m}}{2}(n_{th}+1)\mathcal{D}_{b}[\rho]+\frac{\gamma_{m}n_{th}}{2}\mathcal{D}_{b^{\dagger}}[\rho],$

where $\rho$ is the dynamical density matrix of the system, which can reach
its own steady state $\rho_{s}$ after a long evolution time. $n_{th}$,
$\kappa_{i}$, and $\gamma_{m}$ are the mean phonon number of the molecule
(which depends on the temperature of its heat bath), the optical cavity field
decay (PC and Plasmon cavity), and the molecular damping rate, respectively.
The term $\mathcal{D}_{C}[\rho]$ is referred to as a Lindblad superoperator,
and it can be written as follows:

$\displaystyle\mathcal{D}_{C}[\rho]=2C\rho C^{\dagger}-C^{\dagger}C\rho-\rho
C^{\dagger}C.$ (4)

Eq. (2) is employed to analyze the photon blockade mechanism.

###  3.2 The non-Hermitian Schrödinger equation

To observe the dynamical behavior in our system, the non-Hermitian Schrödinger
equation can be employed. The non-Hermitian Hamiltonian is constructed by
introducing a phenomenological imaginary dissipative term into the effective
Hamiltonian:

$\displaystyle\hat{H^{\prime}}=\hat{H}-\frac{i}{2}(\kappa_{c}\hat{a}_{c}^{\dagger}\hat{a}_{c}+\kappa_{p}\hat{a}_{p}^{\dagger}\hat{a}_{p}+\gamma_{m}\hat{b}^{\dagger}\hat{b}).$
(5)

Where, $\kappa_{c}$, $\kappa_{p}$ are decay rate of photonic cavity, plasmonic
cavity and $\gamma_{m}$ is the molecular damping rate. The term
$ga_{p}^{\dagger}a_{p}(b^{\dagger}+b)$ in the Hermitian part of the effective
Hamiltonian introduces anharmonicity in the energy spectrum of the system. The
degree of anharmonicity depends on the coupling coefficient $g$, with higher
values of $g$ leading to greater anharmonicity and a more pronounced CPB
effect. In the absence of the plasmon-phonon coupling term or when $g$ has a
small value, the CPB effect is not observed. To investigate the UCPB effect,
it’s crucial to consider the quantum interference paths. For this purpose, we
obtain the solution of the non-Hermitian Schrödinger equation. The Schrödinger
equation can be expressed in an arbitrary basis, especially when the Hilbert
space has infinite dimensions. This transformation results in a system of
infinite ordinary differential equations with an infinite number of unknowns.
To solve this system of infinite differential equations, we can employ a
truncation method. The convergence speed of this method is influenced by the
choice of basis. A basis that yields fast convergence is often referred to as
an efficient basis which is similar to optimal wavelets in signal analysis
[27]. Depending on the coupling coefficients $g$ and the strength of $J$,
different preferred bases can be employed, so that four different regime are
arrived as follow:
i) Both $J$ and $g$ are relatively small compared to the cavity decay rate.
The preferred basis consists of the eigenvectors of the non-interacting part
of the Hermitian Hamiltonian, denoted as $\\{|nmq\rangle\\}$. In this case,
the eigenvalues are represented as
$E_{nmq}=n\Delta_{c}+m\Delta_{p}+q\omega_{m}$. The quantum state
$\left|\psi\right\rangle=\sum C_{nmq}|nmq\rangle$ is expressed in this basis,
and the Schrödinger equation takes the form:

$\displaystyle\begin{split}i\dot{C}_{jkl}&=E^{\prime}_{jkl}C_{jkl}+J\sum_{n,m,q}\left[\sqrt{n(m+1)}C_{n-1,m+1,q}+\sqrt{m(n+1)}C_{n+1,m-1,q}\right]\\\
&-g\sum_{n,m,q}\left[m\sqrt{q+1}C_{n,m,q+1}+m\sqrt{q}C_{n,m,q-1}\right]\\\
&+\Omega\sum_{n,m,q}\left[\sqrt{n}C_{n-1,m,q}+\sqrt{n+1}C_{n+1,m,q}\right];\qquad
j,k,l=0,1,2,...\end{split}$ (6)

Where
$E^{\prime}_{jkl}=E_{jkl}-\frac{i}{2}\left(j\kappa_{c}+k\kappa_{p}+l\gamma_{m}\right)$
are the effective Hamiltonian eigenvalues in the absence of couplings. The
eigenvalues $E_{jkl}=j\Delta_{c}+k\Delta_{p}+l\omega_{m}$ are obtained from
the Hamiltonian in the absence of loss terms. The system in Eq. (6) represents
a system of infinite ordinary differential equations with an infinite number
of unknowns. In the absence of laser input ($\Omega$=0) and coupling terms,
the sum of the photon and plasmon number operators commutes with the effective
Hamiltonian, and $N=N_{c}+N_{p}$ is a constant of motion. The anharmonicity
term facilitates energy transfer between phonons and plasmons through
different orders of Stokes and anti-Stokes effects, leading to changes in the
number of plasmons during the system’s evolution. In the weak coupling regime,
these effects are negligible. Furthermore, it is assumed that the bandwidth of
the cavities is so small that neither the Stokes nor the anti-Stokes
frequencies can excite the cavity modes. The photon-plasmon coupling does not
alter the harmonic behavior of the system but transforms the phonon-plasmon
dynamics into positive and negative polariton time evolution. For small values
of $J$, it does not have a significant effect on the conservative observables.
However, the laser input intensity changes the photon and plasmon numbers of
the system, making the conservation of $N$ ineffective for reducing the
dimension of the Hilbert space. For low-intensity laser input
($\Omega\ll\kappa_{c}$), the total number of photons and plasmons is limited
from above, allowing for the truncation of the system of governing equations.
For weak laser input and a small coupling coefficient $g$, the total number of
high-energy particles (photons and plasmons) remains constant, and the basis
$\\{|p,q\rangle;p=n+m\\}$ is efficient for our analysis. The maximum number of
high-energy particles is denoted by p, and the vector space spanned by the
eigenvectors of the Hamiltonian in the absence of input and optomechanical
coupling coefficient, with $p_{i}$ high-energy and q low-energy particles, is
denoted as $V_{p_{i}}$. Therefore, $H=\bigoplus_{p_{i}\leq p}V_{p_{i}}$ is an
efficient Hilbert space for the UCPB effect. $\bigoplus_{p_{i}\leq
p}V_{p_{i}}$ represents the direct sum of vector spaces $V_{p_{i}}$.
ii) When the coupling coefficient between the optical and plasmonic cavities
$J$ is negligible, but the phonon-plasmon coupling coefficient $g$ is not
negligible, the effect of $g$ on anharmonicity is significant. However, its
effect on the photon and plasmon numbers is negligible, making the basis
$\\{|nmq\rangle\\}$ an efficient choice for truncating the governing
equations.
iii) When the coupling coefficient $J$ is not negligible, but $g$ is small.
The transformation $a_{+}=a_{c}\cos\theta+a_{p}\sin\theta$ and
$a_{-}=-a_{c}\sin\theta+a_{p}\cos\theta$
($\tan{2\theta}=\frac{2J}{\Delta_{c}-\Delta_{p}}$) can be used to diagonalize
the cavities part of the Hamiltonian as follows:

$\displaystyle\begin{split}H^{\prime\prime}&=\Delta_{+}a^{\dagger}_{+}a_{+}+\Delta_{-}a^{\dagger}_{-}a_{-}+\omega_{m}b^{\dagger}b-g_{+}a^{\dagger}_{+}a_{+}(b^{\dagger}+b)+-g_{-}a^{\dagger}_{-}a_{-}(b^{\dagger}+b)\\\
&-g_{\pm}(a^{\dagger}_{+}a_{-}+a^{\dagger}_{-}a_{+})(b^{\dagger}+b)+i\Omega_{+}(a^{\dagger}_{+}+a_{+})+i\Omega_{-}(a^{\dagger}_{-}+a_{-}).\end{split}$
(7)

In this context, the coupling strength is defined as follows:
$g_{+}=g\cos^{2}\theta$, $g_{-}=g\sin^{2}\theta$, and
$g_{\pm}=g\sin\theta\cos\theta$. Additionally, the drive of the system is
given by $\Omega_{+}=-\Omega\sin\theta$ and $\Omega_{-}=\Omega\cos\theta$.
Here, the preferred basis consists of the eigenvectors of the Hamiltonian in
the absence of anharmonicity and laser input, denoted as
$\\{|n_{+}n_{-}q\rangle$: $n_{+}$ and $n_{-}$ represent the number of
particles with $\omega_{+}$ and $\omega_{-}$ frequencies, respectively$\\}$.
The detuning frequencies $\Delta_{\pm}=\omega_{\pm}-\omega_{l}$ can be
determined using the following relation:

$\displaystyle\Delta_{\pm}=\frac{1}{2}(\Delta_{c}+\Delta_{p})\pm\frac{1}{2}\sqrt{(\Delta_{c}-\Delta_{p})^{2}+4J^{2}}.$
(8)

When both $\Delta_{c}$ and $\Delta_{p}$ are complex, $\Delta_{+}$ and
$\Delta_{-}$ are also complex and are denoted by $\Delta_{+}^{\prime}$ and
$\Delta_{-}^{\prime}$, by replacing $\Delta_{c}$ and $\Delta_{p}$ with their
complex values $\Delta_{c}^{\prime}=\Delta_{c}-i\kappa_{c}$ and
$\Delta_{p}^{\prime}=\Delta_{p}-i\kappa_{p}$ in Eq. (8). The complex part of
$\Delta_{+}^{\prime}$ is near $\kappa_{p}$, while the complex part of
$\Delta_{-}^{\prime}$ is close to $\kappa_{c}$, hence the particles with
$\omega_{+}$ and $\omega_{-}$ frequencies are called plasmon-like and photon-
like particles. The state vector $\left|\psi\right\rangle=\sum
U_{n_{+}n_{-}q}|n_{+}n_{-}q\rangle$ is written in the new basis and is
substituted to the non-Hermitian Schrödinger equation, the following equations
are obtained:

$\displaystyle\begin{split}i\dot{U}_{jkl}&=F^{\prime}_{jkl}U_{jkl}-g_{+}\sum_{n_{+},n_{+},q}\left[n_{+}\sqrt{q+1}U_{n_{+},n_{-},q+1}+n_{+}\sqrt{q}U_{n_{+},n_{-},q-1}\right]\\\
&-g_{-}\sum_{n_{+},n_{-},q}\left[n_{-}\sqrt{q+1}U_{n_{+},n_{-},q+1}+n_{-}\sqrt{q}U_{n_{+},n_{-},q-1}\right]\\\
&-g_{\pm}\sum_{n_{+},n_{+},q}[\sqrt{(n_{+}+1)n_{-}}\sqrt{q+1}U_{n_{+}+1,n_{-}-1,q+1}+\sqrt{(n_{+}+1)n_{-}}\sqrt{q}U_{n_{+}+1,n_{-}-1,q-1}\\\
&+\sqrt{n_{+}(n_{-}+1)}\sqrt{q+1}U_{n_{+}-1,n_{-}+1,q+1}+\sqrt{n_{+}(n_{-}+1)}\sqrt{q}U_{n_{+}-1,n_{-}+1,q-1}]\\\
&+\Omega_{+}\sum_{n_{+},n_{-},q}\left[\sqrt{n_{+}+1}U_{n_{+}+1,n_{-},q}+\sqrt{n_{+}}U_{n_{+}-1,n_{-},q}\right]\\\
&+\Omega_{-}\sum_{n_{+},n_{-},q}\left[\sqrt{n_{-}+1}U_{n_{+},n_{-}+1,q}+\sqrt{n_{-}}U_{n_{+},n_{-}-1,q}\right];\qquad
j,k,l=0,1,2,....\end{split}$ (9)

Where
$F^{\prime}_{jkl}=F_{jkl}-\frac{i}{2}\left(j\kappa_{+}+k\kappa_{-}+l\gamma_{m}\right)$
are the effective Hamiltonian complex eigenvalues in the absence of couplings.
Eigenvalues $F_{jkl}=j\Delta_{+}+k\Delta_{-}+l\omega_{m}$ are obtained from
the Hamiltonian in the absence of loss terms. Eq. (9) also represents a system
of an infinite number of linear ordinary differential equations with an
infinite number of unknowns. This system of linear ordinary differential
equations can be solved using truncation methods. In the absence of
excitation, the total number of photon-like and plasmon-like particles,
denoted as $N_{p}=a^{\dagger}_{+}a_{+}+a^{\dagger}_{-}a_{-}$, remains a
constant of motion. In the presence of low input intensity
($\Omega_{\text{photon-like}}\ll\kappa_{\text{photon-like}}$), the number of
particles is finite and can be estimated by the number of photons injected
into the optical cavity. With this initial estimation, we can employ the
iteration method to determine a suitable truncation of the system of
differential equations.
iv) When the effects of both coupling coefficients $J$ and $g$ are not
negligible, the photon-plasmon coupling coefficient changes the behavior to
that of photon-like and plasmonic-like particles. While $g$ affects
harmonicity, it has a negligible effect on the total number of polaritons. In
this case, the basis $\\{|n_{+}n_{-}q\rangle\\}$ is efficient for numerical
calculations. This effective basis is independent of $g$ in high-quality
cavities and can be chosen based on the coupling coefficient $J$ between the
photon and plasmon cavities.

##  4 Numerical method and results

The system of governing Eqs. (6) or (9) in any basis is reduced to the
following system of differential equations:

$\displaystyle\dot{X}=AX,$ (10)

where $X$ is the vector of expansion coefficients of the state vector
$\left|\psi\right\rangle$ and $A$ is the matrix of coefficients on the right
hand side of governing equations. Dimensions of vector $X$ and matrix $A$ is
determined by the system parameters and chosen basis. The analytical solution
of Eq. (10) is given :

$\displaystyle X(t)=\exp(At)X(0).$ (11)

To determine the dimension of matrix $A$ and obtain numerical results, the
following algorithm is employed:
a) Based on the values of $J$, $\kappa_{c}$, and $\kappa_{p}$, select an
efficient basis.
b) Estimate the number of high-energy particles (photons and plasmons)
according to the input laser intensity.
c) Use room temperature to estimate the initial value for the phonon number
($q=q_{0}$).
d) Estimate the dimension of the truncated vector space.
e) Arrange the state vector in dictionary order.
f) Determine the elements of matrix $A$.
g) Solve the dynamic equation $\dot{X}=AX$ using analytical or numerical
methods.
h) If $q=q_{0}$, store $X$ in register $X_{0}$ and increment $q$ by 1, then
return to step $d$.
i) Compare $X$ and $X_{0}$ in the sense of the $\mathbb{C}^{N}$ metric. If the
distance between $X$ and $X_{0}$, $d(X,X_{0})$, is greater than a given
$\epsilon$ (epsilon), replace $X_{0}$ with $X$, increment $q$ by 1, and return
to step $d$.
j) If $N=N_{0}$, increment $N$ by 1 and return to step $d$.
k) If $d(X,X_{0})$ is greater than $\epsilon$, replace $X_{0}$ with $X$,
increment $N$ by 1, and return to step $d$.
l) Calculate $g^{(2)}$ and provide the values of $X$, $N$, and $q$.
Our algorithm is employed for various values of the coupling coefficient $J$
and input intensity $\Omega$. Our calculations show that for low input
intensity ($\Omega\ll\kappa_{c}$), the speed of convergence is independent of
the chosen basis. However, for high input intensity, the preferred basis
depends on the value of $J$. For large values of $J$ ($J>\kappa_{p}$), the
preferred basis is $\\{|n_{+}n_{-}q\rangle\\}$, and for small values, no
preferred basis exists.
For weak PC-BNA coupling strength (i and ii regimes), we consider a system
consists of a PC cavity with a central frequency $\omega_{c}=1.342$ eV and a
high quality factor, Q = $10^{6}$ (corresponding to $\kappa_{c}=1.342\times
10^{-3}$ meV). This PC cavity is coupled to a plasmonic nano-antenna with a
central frequency, $\omega_{p}$, of 1.36 eV and a relatively low quality
factor, $Q_{p}=40$ (implying $\kappa_{p}=34$ meV). A Raman-active molecule
with a molecular vibrational energy of $\omega_{m}=200$ meV and a damping rate
of $\gamma_{m}=0.2$ meV is placed in the hot spot of the nano-antenna to
constitute a molecular optomechanics cavity setup. We also considered an
external heat bath with a temperature of $T=300$ K for the molecular system.
Laser field coherently pumps the PC cavity, where the strength of the laser is
taken to be weak ($\Omega$ = 0.01$\kappa_{c}$). We calculated $g^{(2)}(0)$ as
a function of detuning $\Delta_{c}$ and optomechanical coupling $g$ for PC-BNA
coupling strength $J=25$ meV ($J<\kappa_{p}$). To assess the influence of the
basis on convergence speed, these calculations were performed in the bases:
$\\{|nmq\rangle\\}$ with the constraint $n+m=2$ (see Appendix A). As it shown
in Fig. 2(a), the sub-Poissonian distribution is achieved in the blue detuning
and the specific domain of optomechanical coupling which divided in the weak
and strong coupling regime.

(a)

(b)

Figure 2: The corerrelation function $g^{(2)}(0)$ versus BNA-molecule coupling
coefficient $g$ and $\Delta_{c}$ with other parameters set to (a) $J=25$ meV,
$\kappa_{c}=1.342\times 10^{-3},\kappa_{p}=34$ meV,$\omega_{c}=1.342$ eV,
$\omega_{p}=1.36$ eV, $\omega_{m}=200$ meV and $T=300$ K. (b) $J=40$ eV,
$\kappa_{c}=1.342\times 10^{-3},\kappa_{p}=27$ meV,$\omega_{c}=1.32$ eV,
$\omega_{p}=1.36$ eV, $\omega_{m}=150$ meV and $T=300$ K.

For strong PC-BNA coupling strength (iii and iv regimes), we consider almost
the same system with small different in plasmonic cavity ($\omega_{p}=1.32$
eV, $Q_{p}=50$) and molecular frequency ($\omega_{m}=150$
meV,$\gamma_{m}=0.15$ meV). We calculated $g^{(2)}(0)$ as a function of
detuning $\Delta_{c}$ and optomechanical coupling $g$ for PC-BNA coupling
strength $J=40$ meV ($J>\kappa_{p}$). It is worth to mention that, these
calculations were performed in the bases $\\{|n_{+}n_{-}q\rangle\\}$ with the
constraint $n_{+}+n_{-}=2$. As it shown in Fig. 2(b), there is a good region
for single photon generation in weak coupling regime, because of a small
$g^{(2)}(0)\approx 0.02$. We analyze these results with their underlying
physical mechanisms in the next two subsections.

###  4.1 Weak coupling regime

In the domains characterized by optomechanical weak coupling (i, iii regimes),
we plot second order correlation function $g_{c}^{(2)}(0)$ versus $\Delta_{c}$
for $J=25$ meV, $g=20$ meV and $J=40$ meV, $g=13$ meV in Fig. 3(1a) and Fig.
3(2a), respectively. Here, the analytical results (solid line) is compared
with the numerical results (dashed line), where there are almost matched. As
it depicted, there is a dip with value of $g_{c}^{(2)}(0)\approx 0.7$ and
$g_{\text{photon-like}}^{(2)}(0)\approx 0.02$ due to photon blockade effect.
We will explain that these dips occur because of a destructive interference
between two different quantum pathway transitions which called UCPB. To
analyze this effect, we obtain the probabilities of one and two-photon
transition through the following standard method [28].
In the absence of cavity field driving and interaction terms, the number
operator $N=a_{c}^{\dagger}a_{c}+a_{p}^{\dagger}a_{p}$ commutes with the free
Hamiltonian $H_{0}$. We assume that the number of high-energy particles
($p=n+m$) in the system is either less than or equal to $p$, while the number
of phonons can be any quantity. In this context, the Hilbert space
$\mathcal{H}$ of the problem is represented as a direct sum of linear spaces
$V_{i}$ ($i=0,...,p$): $V=\bigoplus_{i=0}^{p}V_{i}$, where $V_{i}$ corresponds
to the vector spaces of high-energy particles. The Hamiltonian restricted to
the space $V_{i}$ is denoted by $H^{(i)}$, where $i=1,...,p$. The eigenvectors
of the Hamiltonian restricted to the space $V$ are denoted by
$|E_{p,q}\rangle$ and can be expressed in the $\\{|nmq\rangle\\}$ basis:
$|E_{p,q}\rangle=\sum_{\begin{subarray}{c}n+m=2\\\
q\end{subarray}}C_{nmq}|nmq\rangle$. The indices $p$ ($p=n+m$) and $q$
correspond to the number of high-energy particles in PC-BNA cavities and the
number of phonons, respectively. Here, because of weak diving condition
($\Omega$=0.01$\kappa_{c}$) we can consider
$\\{C_{200},C_{110},C_{020},C_{201},C_{111},C_{021}\\}\ll\\{C_{100},C_{010},C_{101},C_{011}\\}\ll
C_{000}$, therefore it is estimated that $p=2$ is a good candidate, and the
Hilbert space can be decomposed as $V=V_{0}\oplus V_{1}\oplus V_{2}$. In the
$V_{0}$ space, $H^{(0)}|E_{00}\rangle=E^{(0)}_{00}|E_{00}\rangle$ with the
eigenstate $|E_{00}\rangle=|0,0\rangle$ and the eigenvalue $E^{(0)}_{00}=0$.
In the single-excitation subspace $V_{1}$ and single phonon space, the
Hamiltonian $H^{(1)}$ is becomes:

$\displaystyle H^{(1)}=\begin{pmatrix}\Delta_{c}&J&0&0\\\ J&\Delta_{p}&0&-g\\\
0&0&\Delta_{c}+\omega_{m}&J\\\ 0&-g&J&\Delta_{p}+\omega_{m}\end{pmatrix}.$
(12)

Four eigenvalues and eigenvectors can be obtained by solving the eigenvalue
equation of the $H^{(1)}$ matrix. Each eigenvector $|E_{1,1}\rangle_{l}$
(where $l$ is the index of the eigenvector, $l=1,2,3,4$) can be expressed in
the basis $\\{|nmq\rangle;n+m=1\\}$ as follows:

$\displaystyle|E_{1,1}\rangle_{l}=f_{100}^{(1,l)}|100\rangle+f_{010}^{(1,l)}|010\rangle+f_{101}^{(1,l)}|101\rangle+f_{011}^{(1,l)}|011\rangle,$
(13)

where $f^{(1,l)}_{nmq}$ represents the components of $|E_{1,1}\rangle_{l}$ in
the $\\{|nmq\rangle;n+m=1\\}$ basis. Similarly Hamiltonian $H^{(2)}$ in two-
excitation subspace $V_{2}$ and q=1 is:

$\displaystyle H^{(2)}=\begin{pmatrix}\sqrt{2}J&2\Delta_{c}&0&0&0&0\\\
\Delta_{c}+\Delta_{p}&\sqrt{2}J&\sqrt{2}J&-g&0&0\\\
\sqrt{2}J&0&2\Delta_{p}+2\omega_{m}&0&0&-2g\\\
0&0&0&\sqrt{2}J&2\Delta_{c}+\omega_{m}&0\\\
-g&0&0&\Delta_{c}+\Delta_{p}+\omega_{m}&\sqrt{2}J&\sqrt{2}J\\\
0&0&-2g&\sqrt{2}J&0&2\Delta_{p}+\omega_{m}\end{pmatrix}.$ (14)

Six eigenvectors and eigenvalues can be obtained by the eigenvalue equation of
$H^{(2)}$ matrix. The eigenvectors $|E_{2,1}\rangle_{l};l=1,2,...,6$ versus
the basis $\\{|nmq\rangle;n+m=2\\}$ is:

$\displaystyle|E_{2,1}\rangle_{l}=f_{200}^{(2,l)}|200\rangle+f_{110}^{(2,l)}|110\rangle+f_{020}^{(2,l)}|020\rangle+f_{201}^{(2,l)}|201\rangle+f_{111}^{(2,l)}|111\rangle+f_{021}^{(2,l)}|021\rangle.$
(15)

The state vector $|\psi(t)\rangle$ in the space $V=V_{0}\oplus V_{1}\oplus
V_{2}$ is:

$\displaystyle|\psi(t)\rangle=D_{00}|E_{0,0}\rangle+\sum_{i=1}^{4}D_{1i}|E_{1,1}\rangle_{i}+\sum_{j=1}^{6}D_{2j}|E_{2,1}\rangle_{j}.$
(16)

Figure 3: (1a) The correlation function $g^{(2)}(0)$ (1b) The state
occupations $P_{[|10\rangle]}$ and $P_{[|20\rangle]}$ for the bare states
versus the detuning $\Delta_{c}$ for $J=25$ meV and $g=20$ meV with other
parameters set to $\kappa_{c}=1.342\times 10^{-3}$ meV, $\kappa_{p}=34$ meV,
$\gamma_{m}=200$ meV, $\omega_{c}=1.342$ eV, $\omega_{p}=1.36$ eV (3a) The
correlation function $g^{(2)}(0)$. (2a) The correlation function $g^{(2)}(0)$
(2b) The state occupation probabilities $P_{[|10\rangle]}$ and
$P_{[|20\rangle]}$ for the bare states versus the detuning $\Delta_{c}$ for
$J=40$ meV and $g=13$ meV with other parameters set to $\kappa_{c}=1.32\times
10^{-3}$ meV, $\kappa_{p}=27$ meV, $\gamma_{m}=150$ meV, $\omega_{c}=1.32$ eV,
$\omega_{p}=1.36$ eV.

Eq. (16) is rewritten in the $\\{|nmq\rangle\\}$ basis, and the probability of
the state $|nmq\rangle$ is obtained as a function of the coefficients $D_{ij}$
and $f_{nmq}$. The occupation probability of state $|000\rangle$
($P_{|000\rangle}$) is nearly equal to one, indicating that almost the system
is in this state. By manipulate the system, the probability of transition to
either the one-photon or two-photon state changed. The states with given
photon number ($n$) and plasmon number ($m$), and different phonon numbers
($q_{i}$), are in the same equivalence class
($|nmq_{1}\rangle\sim|nmq_{2}\rangle$) and their equivalence class is denoted
by $[|nm\rangle]=\\{|nmq\rangle;q\in\mathbb{N}\\}$. The occupation probability
of an equivalence class is the sum of the occupation probabilities of all
their elements. Here, we consider only zero and one phonon ($q=0,1$). Each
class represents a group of states around the $p=n+m,q=0$ state. The one-
photon occupation probability can be computed as follows:

$\displaystyle
P_{|100\rangle}=\sum_{i=1}^{4}\left|D_{1i}f_{100}^{1,i}\right|^{2}+\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{4}D_{1i}D_{1j}^{*}f_{100}^{1,i}f_{100}^{1,j*}$
$\displaystyle
P_{|101\rangle}=\sum_{i=1}^{4}\left|D_{1i}f_{101}^{1,i}\right|^{2}+\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{4}D_{1i}D_{1j}^{*}f_{101}^{1,i}f_{101}^{1,j*}$ (17)
$\displaystyle P_{[|10\rangle]}=P_{|100\rangle}+P_{|101\rangle}$

As well as, two-photon occupation probability can be express as:

$\displaystyle
P_{|200\rangle}=\sum_{i=1}^{6}\left|D_{2i}f_{200}^{2,i}\right|^{2}+\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{6}D_{2i}D_{2j}^{*}f_{200}^{2,i}f_{200}^{2,j*}$
$\displaystyle
P_{|201\rangle}=\sum_{i=1}^{6}\left|D_{2i}f_{201}^{2,i}\right|^{2}+\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{6}D_{2i}D_{2j}^{*}f_{201}^{1,i}f_{201}^{1,j*}$ (18)
$\displaystyle P_{[|20\rangle]}=P_{|200\rangle}+P_{|201\rangle}$

In Fig. 3-(1b) and 3-(3b), we depict the occupation states of the bare states
$P_{[|10\rangle]}$ and $P_{[|20\rangle]}$ as functions of the $\Delta_{c}$,
with the weak coupling parameter $g=20$ meV and $g=13$ meV, respectively. The
solid lines (red color for $P_{[|10\rangle]}$ and orange color for
$P_{[|20\rangle]}$) represent the probability when accounting for the
interference terms, whereas the dashed lines (blue color for
$P_{[|10\rangle]}$ and green color for $P_{[|20\rangle]}$) represent the
probability when not considering the interference terms. A dip appears as a
result of the quantum interference effect influencing the state transitions
induced by the driving. This dip signifies the occurrence of destructive
interference, indicating unconventional photon blockade. To establish the
condition for the sub-Poissonian statistics of light, where $g^{(2)}(0)\ll 1$,
it is necessary that the numerator of Eq. (2) be almost zero
($|C_{200}|+|C_{201}|\approx 0$). The Kerr approximation has been also
proposed since the optomechanical coupling is very smaller than mechanical
frequency ($g\ll\omega_{m}$) [13]. In this approximation, the coupling
coefficient $g$ and Kerr parameter $K$ are related by
$K=\frac{g^{2}}{\omega_{m}}$, and the numerical solutions of
$P_{[|20\rangle]}=0$ yield results that closely match the approximation
solution.
Quantum interference occurs between the direct transition
$[|00\rangle]\xrightarrow{\Omega}[|10\rangle]\xrightarrow{\Omega}[|20\rangle]$
and the indirect transition
$[|00\rangle]\xrightarrow{\Omega}[|10\rangle]\xrightarrow{J}[|01\rangle]\xrightarrow{\Omega}[|11\rangle]\xrightarrow{\sqrt{2}J}[|20\rangle]$.
The class of states and corresponding transitions are presented in Fig. 4. The
problem has a strong resemblance to the coupling of a linear cavity to a Kerr
nonlinearity cavity [28]. In this context, the behavior of the molecular
optomechanical system closely resembles the Kerr effect, with the exception of
the energy levels corresponding to molecular phonons around the optical and
plasmonic cavity levels, which distinguish and group them accordingly. The
cross-product terms on the right-hand side of Eq. (17) and Eq. (18) represent
the interference probability term, which plays a crucial role in determining
the UCPB mechanism. For a large PC-BNA coupling $J$, the
${|n_{+}n_{-}q\rangle}$ basis can be employed to determine the path of quantum
interference and the PB effect. The class of states $[|n_{+}n_{-}\rangle]$ and
corresponding transitions are similar to those presented in Fig. 4, with the
replacement of $n$ by $n_{+}$ and $m$ by $n_{-}$.

Figure 4: The eigenenergy spectrum of the coupled cavity system in the
subspace with zero, one, and two excitations.

###  4.2 strong coupling regime

In strong optomechanical coupling (ii, iv regimes), we plot second order
correlation function $g_{c}^{(2)}(0)$ versus $\Delta_{c}$ for $J=25$ meV,
$g=60$ meV and $J=40$ meV, $g=60$ meV in Fig. 5(1a) and Fig. 5(2a),
respectively. Here, the analytical results (solid line) is compared with the
numerical results (dashed line). There is a dip in the regime (ii) with value
of $g^{(2)}(0)\approx 0.02$, due to UCPB photon blockade effect as it shown in
Fig. 5-(1a), in which the probabilities of one and two photons transition are
dependent to quantum interference and hence UCPB happens (see Fig 5-(1b)). In
the strong coupling regime (iv), there are two dips because of CPB mechanism
as it depicted in Fig. 5-(2a). Here, the probabilities of one and two photons
are independent of quantum interference (dips at $\Delta_{c}=-0.05$ eV and
$\Delta_{c}=0.025$) as it shown in Fig. 5-(2b).

Figure 5: (1a) The correlation function $g^{(2)}(0)$ (1b) The state
occupations $P_{[|10\rangle]}$ and $P_{[|20\rangle]}$ for the bare states
versus the detuning $\Delta_{c}$ for $J=25$ meV and $g=60$ meV with other
parameters set to $\kappa_{c}=1.342\times 10^{-3}$ meV, $\kappa_{p}=34$ meV,
$\gamma_{m}=200$ meV, $\omega_{c}=1.342$ eV, $\omega_{p}=1.36$ eV (3a) The
correlation function $g^{(2)}(0)$. (2a) The correlation function $g^{(2)}(0)$
(2b) The state occupation probabilities $P_{[|10\rangle]}$ and
$P_{[|20\rangle]}$ for the bare states versus the detuning $\Delta_{c}$ for
$J=40$ meV and $g=60$ meV with other parameters set to $\kappa_{c}=1.32\times
10^{-3}$ meV, $\kappa_{p}=27$ meV, $\gamma_{m}=150$ meV, $\omega_{c}=1.32$ eV,
$\omega_{p}=1.36$ eV.

##  5 Conclusion

In conclusion, we have studied the potential of an integrated molecular
optomechanical system in a hybrid cavity for room temperature single photon
source. For this point, We have theoretically demonstrated and calculated the
second order correlation function for different coupling regime through the
definition of efficient basis. We have also investigated the probabilities
interference in order to adopt each observed dips to CPB or UCPB mechanisms.
Meanwhile, $g^{(2)}(0)\approx 0.02$ is achieved even in the weak molecular
optomechanics coupling which is practically applicable. Our proposed model can
play a pivotal role as an integral component in photonic based quantum
computer and quantum networks.

## Appendix Appendix A

Substituting the non-Hermitian Hamiltonian given by Eq. (5) and the system
states from Eq. (6) into the Schrödinger equation, we obtain a set of linear
differential equations for the probability amplitudes, as follows:

$\displaystyle-i\dot{C}_{100}$
$\displaystyle=(\Delta_{c}-i\frac{k_{c}}{2})C_{100}+JC_{010}+\Omega(1+\sqrt{2}C_{200})$
$\displaystyle-i\dot{C}_{010}$
$\displaystyle=(\Delta_{p}-i\frac{k_{p}}{2})C_{010}+JC_{100}-gC_{011}+\Omega
C_{110}$ $\displaystyle-i\dot{C}_{110}$
$\displaystyle=(\Delta_{c}-i\frac{k_{c}}{2})C_{110}+(\Delta_{p}-i\frac{k_{p}}{2})C_{110}+\sqrt{2}J(C_{200}+C_{020})-gC_{111}+\Omega
C_{010}$ $\displaystyle-i\dot{C}_{200}$
$\displaystyle=(2\Delta_{c}-i\frac{k_{c}}{2})C_{200}+\sqrt{2}JC_{110}+\sqrt{2}\Omega
C_{100},$ $\displaystyle-i\dot{C}_{020}$
$\displaystyle=(2\Delta_{p}-i\frac{k_{p}}{2})C_{020}+\sqrt{2}JC_{110}-2gC_{021}$
$\displaystyle-i\dot{C}_{201}$
$\displaystyle=(2\Delta_{p}-i\frac{k_{p}}{2})C_{021}+\sqrt{2}JC_{111}-2gC_{020}+(\omega_{m}-i\frac{\gamma_{m}}{2})C_{021}$
(19) $\displaystyle-i\dot{C}_{021}$
$\displaystyle=(2\Delta_{c}-i\frac{k_{c}}{2})C_{201}+\sqrt{2}JC_{111}+\Omega
C_{101}+(\omega_{m}-i\frac{\gamma_{m}}{2})C_{201}$
$\displaystyle-i\dot{C}_{111}$
$\displaystyle=(\Delta_{c}-i\frac{k_{c}}{2})C_{101}+(\omega_{m}-i\frac{\gamma_{m}}{2})C_{101}+JC_{011}+\Omega
C_{111}$ $\displaystyle-i\dot{C}_{011}$
$\displaystyle=(\Delta_{p}-i\frac{k_{p}}{2})C_{011}+(\omega_{m}-i\frac{\gamma_{m}}{2})C_{011}+JC_{101}-gC_{010}+\Omega
C_{111}$ $\displaystyle-i\dot{C}_{101}$
$\displaystyle=(\Delta_{c}-i\frac{k_{c}}{2})C_{111}+(\Delta_{p}-i\frac{k_{p}}{2})C_{111}+(\omega_{m}-i\frac{\gamma_{m}}{2})C_{111}$
$\displaystyle+\sqrt{2}J(C_{021}+C_{201})-gC_{110}+\Omega C_{011}$

To study the steady-state photon statistical properties of the system, we can
disregard changes in amplitude over time.

Funding Acknowledgments Disclosures The authors declare no conflicts of
interest. Data availability Data underlying the results presented in this
paper are not publicly available at this time but may be obtained from the
authors upon reasonable request.

## References

  * [1] N. Luo, M. Zhang, A. Hu, and X. Guo, “Study of single photon source efficiency based on cavity quantum electrodynamics,” Physica B: Condensed Matter 612, 412945 (2021).
  * [2] Y. Ren, S. Duan, W. Xie, Y. Shao, and Z. Duan, “Antibunched photon-pair source based on photon blockade in a nondegenerate optical parametric oscillator,” Physical Review A 103, 053710 (2021).
  * [3] F. Flamini, N. Spagnolo, and F. Sciarrino, “Photonic quantum information processing: a review,” Reports on Progress in Physics 82, 016001 (2018).
  * [4] H. Lin, X. Wang, Z. Yao, and D. Zou, “Kerr-nonlinearity enhanced conventional photon blockade in a second-order nonlinear system,” Optics Express 28, 17643–17652 (2020).
  * [5] Y. Ren, Z. Duan, B. Fan, S. Guan, M. Xie, and C. Liu, “Antibunched single-photon/photon-pair emission with coupled jaynes-cummings model,” Optics Express 30, 21787–21796 (2022).
  * [6] A. Migdall, S. Polyakov, J. Fan, and J. Bienfang, _Single-Photon Generation and Detection: Physics and Applications_ , ISSN (Elsevier Science, 2013).
  * [7] A. Majumdar and D. Gerace, “Single-photon blockade in doubly resonant nanocavities with second-order nonlinearity,” Physical Review B 87, 235319 (2013).
  * [8] S. Abo, G. Chimczak, A. Kowalewska-Kudłaszyk, J. Peřina Jr, R. Chhajlany, and A. Miranowicz, “Hybrid photon–phonon blockade,” Scientific Reports 12, 17655 (2022).
  * [9] X. Xia, X. Zhang, J. Xu, H. Li, Z. Fu, and Y. Yang, “Improvement of nonreciprocal unconventional photon blockade by two asymmetrical arranged atoms embedded in a cavity,” Optics Express 30, 7907–7917 (2022).
  * [10] H. Snijders, J. Frey, J. Norman, H. Flayac, V. Savona, A. Gossard, J. Bowers, M. Van Exter, D. Bouwmeester, and W. Löffler, “Observation of the unconventional photon blockade,” Physical review letters 121, 043601 (2018).
  * [11] A. Ben-Asher, A. I. Fernández-Domínguez, and J. Feist, “Non-hermitian anharmonicity induces single-photon emission,” Physical Review Letters 130, 243601 (2023).
  * [12] D.-Y. Wang, C.-H. Bai, X. Han, S. Liu, S. Zhang, and H.-F. Wang, “Enhanced photon blockade in an optomechanical system with parametric amplification,” Optics Letters 45, 2604–2607 (2020).
  * [13] D.-Y. Wang, C.-H. Bai, S. Liu, S. Zhang, and H.-F. Wang, “Photon blockade in a double-cavity optomechanical system with nonreciprocal coupling,” New Journal of Physics 22, 093006 (2020).
  * [14] H. Flayac, D. Gerace, and V. Savona, “An all-silicon single-photon source by unconventional photon blockade,” Scientific reports 5, 11223 (2015).
  * [15] P. Roelli, C. Galland, N. Piro, and T. J. Kippenberg, “Molecular cavity optomechanics as a theory of plasmon-enhanced raman scattering,” Nature nanotechnology 11, 164–169 (2016).
  * [16] I. Shlesinger, K. G. Cognée, E. Verhagen, and A. F. Koenderink, “Integrated molecular optomechanics with hybrid dielectric–metallic resonators,” Acs Photonics 8, 3506–3516 (2021).
  * [17] R. Esteban, J. J. Baumberg, and J. Aizpurua, “Molecular optomechanics approach to surface-enhanced raman scattering,” Accounts of Chemical Research 55, 1889–1899 (2022).
  * [18] M. K. Dezfouli, R. Gordon, and S. Hughes, “Molecular optomechanics in the anharmonic cavity-qed regime using hybrid metal–dielectric cavity modes,” ACS Photonics 6, 1400–1408 (2019).
  * [19] S. M. Ashrafi, R. Malekfar, A. Bahrampour, and J. Feist, “Long-distance heat transfer between molecular systems through a hybrid plasmonic-photonic nanoresonator,” Journal of Optics 23, 015003 (2020).
  * [20] X. Yang, A. Ishikawa, X. Yin, and X. Zhang, “Hybrid photonic- plasmonic crystal nanocavities,” ACS nano 5, 2831–2838 (2011).
  * [21] K. Santhosh, O. Bitton, L. Chuntonov, and G. Haran, “Vacuum rabi splitting in a plasmonic cavity at the single quantum emitter limit,” Nature communications 7, ncomms11823 (2016).
  * [22] H. Zhang, Y.-C. Liu, C. Wang, N. Zhang, and C. Lu, “Hybrid photonic-plasmonic nano-cavity with ultra-high q/v,” Optics Letters 45, 4794–4797 (2020).
  * [23] Y.-F. Xiao, Y.-C. Liu, B.-B. Li, Y.-L. Chen, Y. Li, and Q. Gong, “Strongly enhanced light-matter interaction in a hybrid photonic-plasmonic resonator,” Physical Review A 85, 031805 (2012).
  * [24] A. I. Barreda, M. Zapata-Herrera, I. M. Palstra, L. Mercadé, J. Aizpurua, A. F. Koenderink, and A. Martínez, “Hybrid photonic-plasmonic cavities based on the nanoparticle-on-a-mirror configuration,” Photonics Research 9, 2398–2419 (2021).
  * [25] M. K. Dezfouli, R. Gordon, and S. Hughes, “Modal theory of modified spontaneous emission of a quantum emitter in a hybrid plasmonic photonic-crystal cavity system,” Physical Review A 95, 013846 (2017).
  * [26] D. Manzano, “A short introduction to the lindblad master equation,” Aip Advances 10 (2020).
  * [27] A. Gupta, S. D. Joshi, and S. Prasad, “A new method of estimating wavelet with desired features from a given signal,” Signal Processing 85, 147–161 (2005).
  * [28] F. Zou, D.-G. Lai, and J.-Q. Liao, “Enhancement of photon blockade effect via quantum interference,” Optics express 28, 16175–16190 (2020).


# Deciphering the long-distance penguin contribution to
$\bar{B}_{d,s}\to\gamma\gamma$ decays

Qin Qina<EMAIL_ADDRESS>Yue-Long Shenb corresponding author:
<EMAIL_ADDRESS>Chao Wangc corresponding author:
<EMAIL_ADDRESS>Yu-Ming Wangd corresponding author:
<EMAIL_ADDRESS>a School of Physics, Huazhong University of Science
and Technology, Luoyu Road 1037, Wuhan Hubei 430074, P.R. China
b College of Information Science and Engineering, Ocean University of China,
Songling Road 238, Qingdao, Shandong 266100, P.R. China
c Department of Mathematics and Physics, Huaiyin Institute of Technology,
Meicheng East Road 1, Huaian, Jiangsu 223200, P.R. China
d School of Physics, Nankai University, Weijin Road 94, Tianjin 300071, P.R.
China

###### Abstract

We compute for the first time the long-distance penguin contribution to the
double radiative $B$-meson decays due to the purely hadronic operators acting
with the electromagnetic current in the background soft-gluon field from first
field-theoretical principles by introducing a novel subleading $B$-meson
distribution amplitude. The numerically dominant penguin amplitude arises from
the soft-gluon radiation off the light up-quark loop rather than the
counterpart charm-loop effect on account of the peculiar analytical behaviour
of the short-distance hard-collinear function. Importantly the long-distance
up-quark penguin contribution brings about the substantial cancellation of the
known factorizable power correction possessing the same multiplication CKM
parameters, thus enabling $B_{d,\,s}\to\gamma\gamma$ to become new benchmark
probes of physics beyond the Standard Model.

##  I Introduction

It is widely accepted that the exclusive radiative penguin bottom-meson decays
play a central role in exploring the quark-flavour dynamics of the Standard
Model (SM) and in probing the nonstandard electroweak interactions at the LHCb
and Belle II experiments. In particular, the double radiative
$\bar{B}_{d,s}\to\gamma\gamma$ decays with non-hadronic final states offer a
remarkably clean environment to address the intricate strong interaction
mechanism of the heavy-hadron system with the perturbative factorization
technique, in comparison with the radiative decays $\bar{B}\to V\gamma$.
Phenomenologically the direct CP asymmetries of the double radiative $B$-meson
decays with the linearly polarized photon states will be also highly
beneficial for determining the CKM phase angle $\gamma$ Bosch and Buchalla
(2002). Applying the QCD factorization approach the leading-power
contributions to the exclusive $\bar{B}_{d,s}\to\gamma\gamma$ decay amplitudes
in the heavy quark expansion have been demonstrated to be factorized into the
short-distance Wilson coefficients due to the hard and hard-collinear
fluctuations as well as the leading-twist bottom-meson distribution amplitude
Descotes-Genon and Sachrajda (2003). In addition, a variety of the subleading-
power corrections to the radiative penguin $\bar{B}_{d,s}\to\gamma\gamma$
decay amplitudes (including the higher-twist off light-cone correction and the
non-leading Fock-state effect) were investigated at tree level in the strong
coupling constant with the diagrammatic factorization approach Shen  _et al._
(2020).

However, the persistent problem of evaluating the long-distance penguin
contribution to the double radiative bottom-meson decay amplitudes in the
presence of soft gluon emission remains unresolved at present. For decades the
non-local subleading power correction arising from the soft gluon radiation
off the charm-loop diagrams constitutes the longstanding obstacle to improve
theory computations for the angular observables of $B\to K^{(\ast)}\ell\ell$
at large hadronic recoil Ali  _et al._ (1991); Khodjamirian  _et al._ (2010,
2013); Hambrock  _et al._ (2015); Bobeth  _et al._ (2018); Kozachuk and
Melikhov (2018); Melikhov (2019); Gubernari  _et al._ (2021) (see Grinstein
and Pirjol (2004); Beylich  _et al._ (2011); Lyon and Zwicky (2014) for more
discussions on such non-local contribution at low hadronic recoil). Achieving
the robust predictions of the long-distance charm-loop effect in the rare
$B\to K^{(\ast)}\ell\ell$ decays will be evidently indispensable for
disentangling the genuine New Physics (NP) effect from the SM background
contribution and for advancing our understanding towards the nature of the
observed flavour anomalies (see for instance Jäger and Martin Camalich (2013);
Descotes-Genon  _et al._ (2013, 2016); Ciuchini  _et al._ (2016); Aebischer
_et al._ (2020); Ciuchini  _et al._ (2021a)). To this end, constructing the
systematic theory formalism to tackle the long-distance penguin contribution
to $\bar{B}_{d,s}\to\gamma\gamma$ will further shed new light on the model-
independent calculation of the analogous QCD corrections to the exclusive
flavour-changing neutral current (FCNC) decays $B\to
K^{(\ast)}\ell\bar{\ell}$. More generally, the newly proposed framework to
cope with the non-local power correction to the double radiative
$\bar{B}_{d,s}$-meson decays will be in the meanwhile of paramount importance
to perform the precision calculation of the radiative and electroweak penguin
decays of heavy-flavour baryons Chen and Geng (2001); He  _et al._ (2006);
Wang  _et al._ (2009a); Ball  _et al._ (2008); Wang  _et al._ (2009b); Mannel
and Wang (2011); Feldmann and Yip (2012); Wang (2012); Braun  _et al._ (2014);
Wang and Shen (2016).

According to the numerical hierarchy between the bottom and charm quark
masses, we will apply the favored power counting scheme $m_{b}\gg
m_{c}\sim{\cal O}(\sqrt{\Lambda\,m_{b}})\gg\Lambda$ as advocated in Grinstein
and Pirjol (2004); Boos  _et al._ (2006a, b); Wang  _et al._ (2017); Gao  _et
al._ (2022) in establishing the perturbative factorization formulae for the
penguin contractions of the effective four-quark operators accompanied by the
soft gluon emission, instead of the alternative counting scheme $m_{b}\sim
m_{c}\gg\Lambda$ implemented in the exclusive two-body $B$-meson decays Beneke
_et al._ (2000); Huber  _et al._ (2016). Subsequently, we will report on a
novel observation on the very hadronic matrix element responsible for the soft
gluon radiation off the penguin diagrams. Integrating out the short-distance
QCD fluctuations embedded in this hadronic quantity will give rise to the
generalized three-particle $B$-meson distribution amplitudes in heavy quark
effective theory (HQET) defined by the non-local matrix element $\langle
0|\bar{q}_{s}(\tau_{1}n)\,G_{\mu\nu}(\tau_{2}\bar{n})\,\Gamma_{i}\,h_{v}(0)|\bar{B}_{v}\rangle$
(in analogy to the subleading shape function $g_{17}(\omega,\omega_{1},\mu)$
discussed in Benzke  _et al._ (2010)) rather than the conventional light-cone
distribution amplitudes as previously introduced in Kawamura  _et al._ (2001);
Braun  _et al._ (2017). Employing the asymptotic behaviour of the generalized
$B$-meson distribution amplitudes at small quark and gluon momenta and the
model-independent theory constraints on these non-perturbative functions, we
will proceed to demonstrate that the soft-collinear convolution integrals
entering the factorized expressions of the long-distance penguin contributions
converge for both the massless-quark and massive-quark loop induced pieces, by
contrast with the corresponding mechanism in the FCNC decay processes $B\to
K^{(\ast)}\ell\ell$. Phenomenological implications of the newly computed power
correction to the double radiative bottom-meson decay observables will be
further explored with the aid of the three-parameter model for the particular
subleading $B$-meson distribution amplitude of our interest.

##  II General analysis

The effective weak Hamiltonian of the double radiative $b\to q\gamma\gamma$
transitions has been shown to be identical to the one for $b\to q\gamma$
decays Grinstein  _et al._ (1990)

$\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}^{8}C_{i}(\nu)\,P_{i}(\nu)\bigg{]}+{\rm h.c.}\,,$

by employing the classical equations of motion Politzer (1980). We will
further adopt the effective operator basis $P_{i}^{(p)}$ as advocated in
Chetyrkin  _et al._ (1997) ensuring the disappearance of Dirac traces
involving an odd number of $\gamma_{5}$ in the subsequent effective theory
computations.

Up to the lowest order in the electromagnetic interaction one can
conventionally cast the exclusive radiative decay amplitude for
$\bar{B}_{q}\to\gamma\gamma$ in the following form Shen  _et al._ (2020)

$\displaystyle{\cal\bar{A}}(\bar{B}_{q}\to\gamma\gamma)$ $\displaystyle=$
$\displaystyle-{4\,G_{F}\over\sqrt{2}}\,{\alpha_{\rm em}\over
4\pi}\,\epsilon^{\ast\alpha}(p)\,\epsilon^{\ast\beta}(q)\,$ (2)
$\displaystyle\times\,\sum_{p=u,c}\,V_{pb}V_{pq}^{\ast}\,\sum_{i=1}^{8}C_{i}\,T_{i,\,\alpha\beta}^{(p)}\,,$

where the yielding hadronic tensors $T_{i,\,\alpha\beta}^{(p)}$ can be
decomposed into the helicity form factors

$\displaystyle T_{i,\,\alpha\beta}^{(p)}$ $\displaystyle=$ $\displaystyle
i\,m_{B_{q}}^{3}\,\bigg{[}\left(g_{\alpha\beta}^{\perp}-i\,\varepsilon_{\alpha\beta}^{\perp}\right)\,F_{i,L}^{(p)}-\left(g_{\alpha\beta}^{\perp}+i\,\varepsilon_{\alpha\beta}^{\perp}\right)\,F_{i,R}^{(p)}\bigg{]},$

thanks to the QED Ward-Takahashi identities and the transversality of the on-
shell photons. Here we have introduced the shorthand notations for brevity

$\displaystyle g_{\alpha\beta}^{\perp}\equiv
g_{\alpha\beta}-{n_{\alpha}\bar{n}_{\beta}\over
2}-{\bar{n}_{\alpha}n_{\beta}\over 2},\hskip
8.5359pt\varepsilon_{\alpha\beta}^{\perp}\equiv{1\over
2}\,\varepsilon_{\alpha\beta\rho\tau}\bar{n}^{\rho}n^{\tau},\hskip 5.69046pt$
(4)

by defining two light-cone vectors $n_{\mu}$ and $\bar{n}_{\mu}$ which satisfy
the kinematic constraints $p_{\mu}=m_{B_{q}}\,\bar{n}_{\mu}/2$ and
$q_{\mu}=m_{B_{q}}\,n_{\mu}/2$. It is interesting to note that only the left-
handed form factors $F_{i,L}^{(p)}$ will survive at leading order in the heavy
quark expansion on account of the helicity conservation of the QCD interaction
at high energy. Explicitly, the resulting factorization formula for
$F_{i,L}^{(p)}$ at leading power can then be written as

$\displaystyle\sum_{i=1}^{8}C_{i}\,F_{i,L}^{(p),\,{\rm
LP}}=-{Q_{q}\,f_{B_{q}}\,\overline{m}_{b}(\nu)\over
m_{B_{q}}}\,\,V_{7,\,\rm{eff}}^{(p)}(m_{b},\mu,\nu)\,$ $\displaystyle\hskip
34.14322pt\times\,K^{-1}(m_{b},\mu)\,\int_{0}^{\infty}\,{d\omega\over\omega}\,\phi_{B}^{+}(\omega,\mu)\,J(m_{b},\omega,\mu),\hskip
22.76228pt$ (5)

where the desired expressions for the effective hard function
$V_{7,\,\rm{eff}}^{(p)}$, the perturbative matching coefficient $K$ and the
hard-collinear function $J$ at the one-loop accuracy can be found in Shen  _et
al._ (2020); Eichten and Hill (1990); Lunghi  _et al._ (2003); Bosch  _et al._
(2003).

##  III QCD factorization for the long-distance penguin contribution

Figure 1: Diagrammatical representation of the soft-gluon emission from the
factorizable quark loop generated by the effective four-quark operator and the
electromagnetic current in the double radiative $\bar{B}_{q}\to\gamma\gamma$
decay process, where the symmetric diagram due to the exchange of two on-shell
photons is not presented.

We are now in a position to explore factorization properties of the long-
distance penguin contribution to the double radiative
$\bar{B}_{q}\to\gamma\gamma$ decay amplitude by inspecting the partonic
diagram displayed in Figure 1. Integrating out the hard-collinear quark loop
one can readily derive the flavour-changing scattering amplitude of
$g(\ell)+b(v)\to q(\tilde{q})+\gamma(p)$ governed by the effective weak
Hamiltonian (1) by discarding the subleading-power terms in $\Lambda/m_{b}$
and by invoking the on-shell constraint for the external photon state

$\displaystyle{\cal M}(g+b\to q+\gamma)=i\,{4\,G_{F}\over\sqrt{2}}\,{g_{\rm
em}g_{s}\over 4\pi^{2}}\,\sum_{p=u,c}\,V_{pb}V_{pq}^{\ast}\,\bigg{\\{}$
$\displaystyle\left(C_{2}-{C_{1}\over
2N_{c}}\right)\,Q_{p}\,\left[F(z_{p})-1\right]+6\,C_{6}\,\sum_{q^{\prime}}\,Q_{q^{\prime}}\,\left[F(z_{q^{\prime}})-1\right]$
$\displaystyle+\left[\left(C_{3}-{C_{4}\over
2N_{c}}\right)+16\,\left(C_{5}-{C_{6}\over
2N_{c}}\right)\right]\,\,Q_{q}\,\left[F(z_{q})-1\right]\bigg{\\}}$
$\displaystyle\times\,\left[\bar{q}(\tilde{q})\,\gamma_{\beta}\,P_{L}\,G_{\mu\alpha}\,\,\tilde{F}^{\mu\beta}b(v)\right]\,\frac{p^{\alpha}}{(p-\ell)^{2}}\,,$
(6)

where the perturbative penguin function is given by

$\displaystyle F(x)=4\,x\,\arctan^{2}\left(\frac{1}{\sqrt{4\,x-1}}\right)\,,$
(7)

and we have further employed the conventions

$\displaystyle
z_{p}={m_{p}^{2}-i\,0^{+}\over(p-\ell)^{2}}\,,\qquad\tilde{F}^{\mu\nu}=-\left({1\over
2}\right)\,\epsilon^{\mu\nu\alpha\beta}\,F_{\alpha\beta}\,.$ (8)

Apparently, the hard-scattering kernel of the partonic amplitude displayed in
(6) depends on the unique component $\bar{n}\cdot\ell$ of the soft-gluon
momentum at leading order in the heavy quark expansion. Moreover, it is
straightforward to verify that the bottom-quark loop diagram with an insertion
of the QCD penguin operator $P_{6}$ can only bring about the subleading-power
effect, by virtue of the asymptotic behaviour of $F(x)-1\sim{\cal O}({1/x})$
at large $x$, as it should be. It remains necessary to remark that the long-
distance penguin contraction mechanism does not give rise to the non-trivial
strong phase for the radiative $\bar{B}_{q}\to\gamma\gamma$ amplitude in
consequence of the space-like four-momentum $(p-\ell)$, by contrast with the
counterpart contribution to the radiative leptonic
$\bar{B}_{q}\to\gamma\ell\bar{\ell}$ decay.

We can proceed to evaluate the five-point amplitude of
$g(\ell)+b(v)+\bar{q}(k)\to\gamma(p)+\gamma(q)$ with the diagrammatic
factorization technique by integrating out the anti-hard-collinear quark
propagator in Figure 1 subsequently

$\displaystyle\langle\gamma(p)\,\gamma(q)|\bar{q}\,\gamma_{\beta}\,P_{L}\,G_{\mu\alpha}\,\,\tilde{F}^{\mu\beta}\,b|g(\ell)\,b(v)\,\bar{q}(k)\rangle$
$\displaystyle\Rightarrow{i\,g_{\rm
em}\,e_{q}\over(q-k)^{2}}\,\epsilon^{\mu\beta\lambda\tau}\,p_{\lambda}\,\epsilon_{\tau}^{\ast}(p)\,\epsilon_{\rho}^{\ast}(q)\,$
$\displaystyle\hskip
11.38092pt\times\,\left[\bar{q}(k)\,\gamma_{\perp}^{\rho}\not\\!{q}\,\gamma_{\beta}\,P_{L}\,G_{\mu\alpha}(\ell)\,b(v)\right]\,+{\cal
O}(\alpha_{s})\,,$ (9)

where the yielding short-distance matching coefficient depends on the
component $n\cdot k$ (rather than $\bar{n}\cdot k$) of the soft-quark momentum
in the leading-power approximation. As a consequence, it becomes evident to
introduce the subleading $B$-meson distribution amplitude defined by the HQET
matrix element of the three-body non-local operator
$\bar{q}_{s}(\tau_{1}n)\,G_{\mu\nu}(\tau_{2}\bar{n})\,\Gamma_{i}\,h_{v}(0)$
for the sake of describing the soft QCD dynamics encoded in the long-distance
penguin contribution to $\bar{B}_{q}\to\gamma\gamma$. Constructing the general
parametrization of the emerged effective matrix element with the covariant
tensor formalism Falk  _et al._ (1990) (in analogy to the Lorentz
decomposition for the corresponding light-cone matrix element Braun  _et al._
(2017)) allows us to derive the soft-collinear factorization formula of the
soft-gluon radiative correction to the left-handed helicity form factor

$\displaystyle\sum_{i=1}^{8}C_{i}\,F_{i,L}^{(p),\,{\rm
soft\,4q}}=-{Q_{q}\,f_{B_{q}}\over
m_{B_{q}}}\,\int_{0}^{\infty}{d\omega_{1}\over\omega_{1}}\,\int_{0}^{\infty}{d\omega_{2}\over\omega_{2}}\,$
$\displaystyle\,\,\bigg{\\{}\left(C_{2}-{C_{1}\over
2N_{c}}\right)\,Q_{p}\,\left[F(z_{p})-1\right]+6\,C_{6}\,Q_{c}\,\,\left[F(z_{c})-1\right]$
$\displaystyle\,\,-\left[\left(C_{3}-{C_{4}\over
2N_{c}}\right)+16\,\left(C_{5}-{C_{6}\over
2N_{c}}\right)\right]\,\,Q_{q}\bigg{\\}}\,$
$\displaystyle\,\,\times\,\Phi_{\rm G}(\omega_{1},\omega_{2},\mu)+\,{\cal
O}(\alpha_{s})\,,\hskip 14.22636pt$ (10)

where the electric-charge relation for the light-flavour quarks
$Q_{u}+Q_{d}+Q_{s}=0$ and the vanishing penguin function $F(0)$ have been
applied to simplify the obtained perturbative kernel. It is customary to
define two light-cone variables $\omega_{1}=n\cdot k$ and
$\omega_{2}=\bar{n}\cdot\ell$ such that the resulting hard-collinear function
develops a peculiar dependence on $\omega_{2}$ via the dimensionless quantity
$z_{p}=-m_{p}^{2}/(m_{B}\,\omega_{2})$ (apart from an overall factor
$1/\omega_{2}$). Including the higher-order QCD corrections to the partonic
diagrams in Figure 1 will generate the non-trivial hard functions (instead of
“$1$” at tree level) from matching the effective four-quark operators
$P_{i}^{(p)}$ onto ${\rm SCET}_{\rm I}$ and simultaneously result in the
interesting impacts on the (anti)-hard-collinear matching coefficients
appeared in (10). However, the fundamental property that the determined short-
distance matching functions merely depend on the two dimensional variables
$\omega_{1}$ and $\omega_{2}$ remains valid beyond the leading-order accuracy,
thus justifying the appearance of the soft matrix element $\langle
0|\bar{q}_{s}(\tau_{1}n)\,G_{\mu\nu}(\tau_{2}\bar{n})\,\Gamma_{i}\,h_{v}(0)|\bar{B}_{v}\rangle$
(see Beneke  _et al._ (2020) for further discussions in a different context).
Schematically, the factorized expression for the long-distance penguin
correction to $\bar{B}_{q}\to\gamma\gamma$ can be cast in the form ${\cal
H}\,{\cal J}\star\bar{\cal J}\star\Phi_{\rm G}$, which resembles the very
pattern for the $Q_{1}^{q}-Q_{7\gamma}$ contribution to $\bar{B}\to
X_{s}\gamma$ Benzke  _et al._ (2010). We also mention in passing that the
newly computed long-distance soft gluon radiative correction appears to
preserve the large-recoil symmetry of the two transversality amplitudes and
does not affect the right-handed form factors $F_{i,R}^{(p)}$, in agreement
with the earlier observation on the “resolved photon” contribution to
$\bar{B}\to\gamma\ell\nu$ at twist-two Wang and Shen (2018) (see also Braun
and Khodjamirian (2013); Wang (2016) for an estimate of the soft contribution
with the dispersion approach).

Additionally, the novel subleading distribution amplitude (perhaps more
appropriately called soft function) of the $B$-meson $\Phi_{\rm G}$ in (10) is
defined in terms of the effective matrix element of the three-body non-local
operator with quark and gluon fields localized on two distinct light-cone
directions

$\displaystyle\langle
0|(\bar{q}_{s}S_{n})(\tau_{1}n)\,(S_{n}^{\dagger}\,S_{\bar{n}})(0)\,(S_{\bar{n}}^{\dagger}\,g_{s}\,G_{\mu\nu}\,S_{\bar{n}})(\tau_{2}\bar{n})\,\,\bar{n}^{\nu}\not\\!n\gamma_{\perp}^{\mu}\gamma_{5}\,(S_{\bar{n}}^{\dagger}h_{v})(0)|\bar{B}_{v}\rangle$
$\displaystyle=2\,\tilde{f}_{B}(\mu)\,m_{B}\,\int_{0}^{\infty}d\omega_{1}\,\int_{0}^{\infty}d\omega_{2}\,{\rm
exp}\left[-i(\omega_{1}\tau_{1}+\omega_{2}\tau_{2})\right]\,\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu)\,,$ (11)

where the two soft Wilson lines $S_{n}$ and $S_{\bar{n}}$ essential to
maintain gauge invariance are given by

$\displaystyle S_{n}(x)$ $\displaystyle=$ $\displaystyle{\rm
P\,exp}\left[i\,g_{s}\int_{-\infty}^{0}dt\,n\cdot A_{s}(x+t\,n)\right]\,,$
$\displaystyle S_{\bar{n}}(x)$ $\displaystyle=$ $\displaystyle{\rm
P\,exp}\left[i\,g_{s}\int_{-\infty}^{0}dt\,\bar{n}\cdot
A_{s}(x+t\,\bar{n})\right]\,.$ (12)

It is apparent that the non-local HQET matrix element on the left-hand side of
(11) can be described by the familiar three-particle light-cone distribution
amplitude when taking the limit $\tau_{1(2)}\to 0$. We are therefore led to
the following three important and model-independent normalization conditions
at tree level

$\displaystyle\int_{0}^{\infty}d\omega_{1}\,\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu)=\int_{0}^{\infty}d\omega_{1}\,\Phi_{4}(\omega_{1},\omega_{2},\mu)\,,$
$\displaystyle\int_{0}^{\infty}d\omega_{2}\,\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu)=\int_{0}^{\infty}d\omega_{2}\,\Phi_{5}(\omega_{1},\omega_{2},\mu)\,,$
$\displaystyle\int_{0}^{\infty}d\omega_{1}\,\int_{0}^{\infty}d\omega_{2}\,\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu)=\frac{\lambda_{E}^{2}+\lambda_{H}^{2}}{3}\,,$
(13)

where the explicit definitions of the twist-four and twist-five light-cone
distribution amplitudes $\Phi_{4}$ and $\Phi_{5}$ can be found in Braun  _et
al._ (2017) and the hadronic quantities ${\lambda}_{E}^{2}$ and
${\lambda}_{H}^{2}$ can be defined by the effective matrix elements of the
local chromoelectric and chromomagnetic operators Grozin and Neubert (1997).
Furthermore, the asymptotic behaviour of $\Phi_{\rm G}$ at small quark and
gluon momenta can be predicted with the dispersion technique as widely adopted
in the explorations of the two-particle and three-particle $B$-meson light-
cone distribution amplitudes Grozin and Neubert (1997); Braun  _et al._
(2004); Khodjamirian  _et al._ (2007); Lü  _et al._ (2019); Khodjamirian  _et
al._ (2020). Starting with the HQET correlation function

$\displaystyle\Pi_{\rm G}=i\,\int d^{4}x\,{\rm exp}\left(-i\,\omega\,v\cdot
x\right)\,\langle 0|{\rm T}\big{\\{}\big{[}(\bar{q}_{s}S_{n})(\tau_{1}n)$
$\displaystyle\,\,(S_{n}^{\dagger}\,S_{\bar{n}})(0)\,(S_{\bar{n}}^{\dagger}\,g_{s}\,G_{\mu\nu}\,S_{\bar{n}})(\tau_{2}\bar{n})\,\,\bar{n}^{\nu}\not\\!n\gamma_{\perp}^{\mu}\gamma_{5}\,(S_{\bar{n}}^{\dagger}h_{v})(0)\big{]},$
$\displaystyle\,\,\left[\bar{h}_{v}(x)\,\,g_{s}\,G_{\rho\lambda}(x)\,\sigma^{\rho\lambda}\,\gamma_{5}\,q_{s}(x)\right]\big{\\}}|0\rangle\,,$
(14)

we can on the one hand compute this quantity in the kinematic region
$|\omega|\gg\Lambda$ with the operator-product-expansion (OPE) method and on
the other hand derive the corresponding hadronic representation of $\Pi_{\rm
G}$ by taking advantage of analyticity with respect to the effective variable
$\omega$. Matching the above two dispersion representations with the aid of
the parton-hadron duality ansatz enables us to extract the desired asymptotic
behaviour $\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu)\sim\omega_{1}\,\omega_{2}^{2}$ at
$\omega_{1},\,\omega_{2}\to 0$ immediately, which further indicates that the
convolution integrals in the factorized expression (10) converge. We restrict
ourselves to the leading-order accuracy in this letter, however, it will be of
interest to investigate the non-trivial impact of the renormalization-group
evolution on the generalized distribution amplitude $\Phi_{\rm G}$ in the
future.

##  IV Numerical implications

We now turn to address the phenomenological implications of the soft-gluon
radiative correction to the penguin contractions of the effective four-quark
operators on the double radiative $\bar{B}_{q}\to\gamma\gamma$ decay
amplitudes. To achieve this goal, we first need to construct the acceptable
non-perturbative model for the subleading distribution amplitude $\Phi_{\rm
G}$ fulfilling the third relation in (13) as well as the obtained asymptotic
behaviour

$\displaystyle\Phi_{\rm
G}(\omega_{1},\omega_{2},\mu_{0})=\frac{\lambda_{E}^{2}+\lambda_{H}^{2}}{6}\,{\omega_{1}\omega_{2}^{2}\over\omega_{0}^{5}}\,{\rm
exp}\left(-{\omega_{1}+\omega_{2}\over\omega_{0}}\right)$ $\displaystyle\hskip
42.67912pt\frac{\Gamma(\beta+2)}{\Gamma(\alpha+2)}\,\,U\left(\beta-\alpha,4-\alpha,{\omega_{1}+\omega_{2}\over\omega_{0}}\right)\,,$
(15)

at the reference scale $\mu_{0}=1.0\,{\rm GeV}$, motivated from the suggested
three-parameter ansatz for the twist-two distribution amplitude
$\phi_{B}^{+}(\omega,\mu_{0})$ Beneke  _et al._ (2018a). The remaining shape
parameters $\omega_{0}$, $\alpha$ and $\beta$ can be further determined by
enforcing the first and second normalization relations in (13) and employing
the concrete model of the twist-four and twist-five light-cone distribution
amplitudes $\Phi_{4,\,5}(\omega_{1},\omega_{2},\mu_{0})$ as implemented in Gao
_et al._ (2022). The conventional HQET distribution amplitudes on the light-
cone appearing in the established factorization formulae of the helicity form
factors $F_{i,L(R)}^{(p)}$ from Shen  _et al._ (2020) will be also in demand
in the subsequent numerical investigations, and we will apply the same
phenomenological model as presented in this reference, with the exceptions of
updated intervals for the inverse moments $\lambda_{B_{d}}=(275\pm 75)\,{\rm
MeV}$ and $\lambda_{B_{s}}=(325\pm 75)\,{\rm MeV}$ Beneke  _et al._ (2018b)
(see Khodjamirian  _et al._ (2020) for a recent determination of
$\lambda_{B_{s}}/\lambda_{B_{d}}$ with the method of QCD sum rules). The
allowed intervals of additional theory input parameters entering our numerical
studies are identical to the ones collected in Shen  _et al._ (2020).

Figure 2: Theory predictions for the soft-gluon radiative corrections to the
left-handed helicity form factors of $\bar{B}_{q}\to\gamma\gamma$ with the
perturbative uncertainties from varying the factorization scale in the
interval $\mu\in[1.0,\,2.0]\,{\rm GeV}$ (green bands), where we further
display the numerical results for the combined factorizable power corrections
as previously derived in Shen  _et al._ (2020) (blue bands) for a comparison.
The shorthand notations for the weighted helicity form factors ${\cal
F}_{L}^{(p),X}=\sum\limits_{i}\,C_{i}\,F_{i,L}^{(p),X}$ have been introduced
here for convenience.

We present the theory predictions for the long-distance penguin contribution
to the double radiative $B$-meson decay form factors including the obtained
uncertainties from varying the factorization scale $\mu$ in Figure 2, where
the numerical predictions of the previously computed factorizable power
corrections at tree level Shen  _et al._ (2020) are also presented for the
illustration purpose. Interestingly, the newly computed soft gluon radiation
from the up-quark penguin contraction appears to generate the substantial
cancellation of the combined factorizable power corrections derived in Shen
_et al._ (2020). On the contrary, the long-distance charm-quark penguin
mechanism will only lead to the rather minor impact on the left-handed
helicity form factor $\sum\limits_{i}\,C_{i}\,F_{i,L}^{(c)}$ numerically as
indicated by Figure 2. This intriguing observation can be attributed to the
peculiar analytical behaviour of the perturbative penguin function $F(z_{p})$
entering the soft-collinear factorization formula (10) such that the yielding
result of $\left|F\left(-m_{c}^{2}/(m_{B}\,\omega_{2})\right)-1\right|$ in the
bulk of the integral domain of $\omega_{2}$ is approximately one order of
magnitude lower than the corresponding up-quark penguin contribution
$\left|F(0)-1\right|=1$. Along the same vein, one can readily observe that the
soft-gluon radiative corrections to the factorizable quark loops cannot
generate numerically important contributions to the exclusive
$\bar{B}_{s}\to\gamma\gamma$ helicity amplitudes due to the CKM suppression of
the up-quark penguin contraction and the dynamical suppression of the charm-
quark penguin contribution as discussed above. Bearing in mind the utmost
importance of understanding the charming penguin contribution in unveiling the
genuine NP effects embedded in the semileptonic $B\to
K^{(\ast)}\ell\bar{\ell}$ decays (see for instance Cerri  _et al._ (2019);
Ciuchini  _et al._ (2021b); Altmannshofer and Archilli (2022)), the achieved
robust control of such long-distance penguin contribution in the double
radiative $B$-meson decays with the diagrammatic factorization technique
evidently makes these FCNC decay processes most suitable for probing the
nonstandard four-fermion $b\to q\,f\bar{f}$ interaction at the high-luminosity
Belle II experiment Altmannshofer  _et al._ (2019) and for performing the
dedicated parton tomography of the composite heavy-quark hadron system in the
QCD framework.

We are now ready to explore the numerical impacts of the long-distance penguin
contribution on the CP-averaged branching fractions, the two polarization
fractions and the CP-violating observables for $\bar{B}_{q}\to\gamma\gamma$
(see Shen  _et al._ (2020) for their explicit definitions). It turns out that
such subleading power corrections can enhance the theory predictions for the
mixing induced CP asymmetries ${\cal A}_{\rm CP}^{\rm mix,\,\|}$ and ${\cal
A}_{\rm CP}^{\rm mix,\,\perp}$ by approximately an amount of ${\cal
O}(30\,\%)$ with the default inputs, while yielding insignificant effects in
the remaining observables numerically. Moreover, our numerical prediction for
the ratio of the two branching fractions ${\cal
BR}(B_{s}\to\gamma\gamma):{\cal BR}(B_{d}\to\gamma\gamma)$ allows for
extracting the high-profile hadronic quantity
$\lambda_{B_{s}}:\lambda_{B_{d}}$ with the improved systematic uncertainty at
the level of $(5-10)\%$, when confronting with the anticipated precision
measurements.

##  V Conclusions

In conclusion, we have presented the first computation of the long-distance
penguin contribution to the double radiative $B$-meson decay amplitudes by
applying the perturbative factorization approach. Adopting the power counting
scheme $m_{b}\gg m_{c}\sim{\cal O}(\sqrt{\Lambda\,m_{b}})\gg\Lambda$, we
demonstrated further that the novel subleading $B$-meson distribution
amplitude $\Phi_{\rm G}$ defined by the three-body HQET operator with partonic
fields localized on two different light-ray directions (instead of the
conventional light-cone distribution amplitude) emerged naturally in the
resulting factorization formula (10). Phenomenologically the soft-gluon
radiative off the factorizable up-quark loop appeared to bring about more
pronounced effect in comparison with the corresponding charm-quark penguin
mechanism thanks to the peculiar analytical behaviour of the short-distance
matching function (7). In addition, the observed destructive interference
between the long-distance penguin contribution and the available factorizable
power correction from Shen  _et al._ (2020) enabled us to determine the
pivotal ratio of the two inverse moments $\lambda_{B_{d,s}}$ with the reduced
theory uncertainty. Our analysis will be evidently beneficial for exploring
the intricate charming penguin dynamics encoded in a large variety of the
radiative and electroweak penguin decay processes including $B\to
K^{\ast}\gamma$ and $B\to K^{(\ast)}\ell\bar{\ell}$, which are generally
recognized as the flagship probes of physics beyond the SM at the LHC.

###### Acknowledgements.

## Acknowledgements

The research of Q.Q. is supported by the National Natural Science Foundation
of China with Grant No. 12005068. 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.

## References

  * Bosch and Buchalla (2002) S. W. Bosch and G. Buchalla, JHEP 08, 054 (2002), arXiv:hep-ph/0208202 .
  * Descotes-Genon and Sachrajda (2003) S. Descotes-Genon and C. T. Sachrajda, Phys. Lett. B 557, 213 (2003), arXiv:hep-ph/0212162 .
  * Shen  _et al._ (2020) Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, JHEP 12, 169 (2020), arXiv:2009.02723 [hep-ph] .
  * Ali  _et al._ (1991) A. Ali, T. Mannel, and T. Morozumi, Phys. Lett. B 273, 505 (1991).
  * Khodjamirian  _et al._ (2010) A. Khodjamirian, T. Mannel, A. A. Pivovarov, and Y. M. Wang, JHEP 09, 089 (2010), arXiv:1006.4945 [hep-ph] .
  * Khodjamirian  _et al._ (2013) A. Khodjamirian, T. Mannel, and Y. M. Wang, JHEP 02, 010 (2013), arXiv:1211.0234 [hep-ph] .
  * Hambrock  _et al._ (2015) C. Hambrock, A. Khodjamirian, and A. Rusov, Phys. Rev. D 92, 074020 (2015), arXiv:1506.07760 [hep-ph] .
  * Bobeth  _et al._ (2018) C. Bobeth, M. Chrzaszcz, D. van Dyk, and J. Virto, Eur. Phys. J. C 78, 451 (2018), arXiv:1707.07305 [hep-ph] .
  * Kozachuk and Melikhov (2018) A. Kozachuk and D. Melikhov, Phys. Lett. B 786, 378 (2018), arXiv:1805.05720 [hep-ph] .
  * Melikhov (2019) D. Melikhov, EPJ Web Conf. 222, 01007 (2019), arXiv:1911.03899 [hep-ph] .
  * Gubernari  _et al._ (2021) N. Gubernari, D. van Dyk, and J. Virto, JHEP 02, 088 (2021), arXiv:2011.09813 [hep-ph] .
  * Grinstein and Pirjol (2004) B. Grinstein and D. Pirjol, Phys. Rev. D 70, 114005 (2004), arXiv:hep-ph/0404250 .
  * Beylich  _et al._ (2011) M. Beylich, G. Buchalla, and T. Feldmann, Eur. Phys. J. C 71, 1635 (2011), arXiv:1101.5118 [hep-ph] .
  * Lyon and Zwicky (2014) J. Lyon and R. Zwicky, (2014), arXiv:1406.0566 [hep-ph] .
  * Jäger and Martin Camalich (2013) S. Jäger and J. Martin Camalich, JHEP 05, 043 (2013), arXiv:1212.2263 [hep-ph] .
  * Descotes-Genon  _et al._ (2013) S. Descotes-Genon, J. Matias, and J. Virto, Phys. Rev. D 88, 074002 (2013), arXiv:1307.5683 [hep-ph] .
  * Descotes-Genon  _et al._ (2016) S. Descotes-Genon, L. Hofer, J. Matias, and J. Virto, JHEP 06, 092 (2016), arXiv:1510.04239 [hep-ph] .
  * Ciuchini  _et al._ (2016) M. Ciuchini, M. Fedele, E. Franco, S. Mishima, A. Paul, L. Silvestrini, and M. Valli, JHEP 06, 116 (2016), arXiv:1512.07157 [hep-ph] .
  * Aebischer  _et al._ (2020) J. Aebischer, W. Altmannshofer, D. Guadagnoli, M. Reboud, P. Stangl, and D. M. Straub, Eur. Phys. J. C 80, 252 (2020), arXiv:1903.10434 [hep-ph] .
  * Ciuchini  _et al._ (2021a) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, Phys. Rev. D 103, 015030 (2021a), arXiv:2011.01212 [hep-ph] .
  * Chen and Geng (2001) C.-H. Chen and C. Q. Geng, Phys. Rev. D 64, 074001 (2001), arXiv:hep-ph/0106193 .
  * He  _et al._ (2006) X.-G. He, T. Li, X.-Q. Li, and Y.-M. Wang, Phys. Rev. D 74, 034026 (2006), arXiv:hep-ph/0606025 .
  * Wang  _et al._ (2009a) Y.-M. Wang, Y. Li, and C.-D. Lu, Eur. Phys. J. C 59, 861 (2009a), arXiv:0804.0648 [hep-ph] .
  * Ball  _et al._ (2008) P. Ball, V. M. Braun, and E. Gardi, Phys. Lett. B 665, 197 (2008), arXiv:0804.2424 [hep-ph] .
  * Wang  _et al._ (2009b) Y.-M. Wang, Y.-L. Shen, and C.-D. Lu, Phys. Rev. D 80, 074012 (2009b), arXiv:0907.4008 [hep-ph] .
  * Mannel and Wang (2011) T. Mannel and Y.-M. Wang, JHEP 12, 067 (2011), arXiv:1111.1849 [hep-ph] .
  * Feldmann and Yip (2012) T. Feldmann and M. W. Y. Yip, Phys. Rev. D 85, 014035 (2012), [Erratum: Phys.Rev.D 86, 079901 (2012)], arXiv:1111.1844 [hep-ph] .
  * Wang (2012) W. Wang, Phys. Lett. B 708, 119 (2012), arXiv:1112.0237 [hep-ph] .
  * Braun  _et al._ (2014) V. M. Braun, S. E. Derkachov, and A. N. Manashov, Phys. Lett. B 738, 334 (2014), arXiv:1406.0664 [hep-ph] .
  * Wang and Shen (2016) Y.-M. Wang and Y.-L. Shen, JHEP 02, 179 (2016), arXiv:1511.09036 [hep-ph] .
  * Boos  _et al._ (2006a) H. Boos, T. Feldmann, T. Mannel, and B. D. Pecjak, Phys. Rev. D 73, 036003 (2006a), arXiv:hep-ph/0504005 .
  * Boos  _et al._ (2006b) H. Boos, T. Feldmann, T. Mannel, and B. D. Pecjak, JHEP 05, 056 (2006b), arXiv:hep-ph/0512157 .
  * 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] .
  * Gao  _et al._ (2022) J. Gao, T. Huber, Y. Ji, C. Wang, Y.-M. Wang, and Y.-B. Wei, JHEP 05, 024 (2022), arXiv:2112.12674 [hep-ph] .
  * Beneke  _et al._ (2000) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, Nucl. Phys. B 591, 313 (2000), arXiv:hep-ph/0006124 .
  * Huber  _et al._ (2016) T. Huber, S. Kränkl, and X.-Q. Li, JHEP 09, 112 (2016), arXiv:1606.02888 [hep-ph] .
  * Benzke  _et al._ (2010) M. Benzke, S. J. Lee, M. Neubert, and G. Paz, JHEP 08, 099 (2010), arXiv:1003.5012 [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 .
  * Braun  _et al._ (2017) V. M. Braun, Y. Ji, and A. N. Manashov, JHEP 05, 022 (2017), arXiv:1703.02446 [hep-ph] .
  * Grinstein  _et al._ (1990) B. Grinstein, R. P. Springer, and M. B. Wise, Nucl. Phys. B 339, 269 (1990).
  * Politzer (1980) H. D. Politzer, Nucl. Phys. B 172, 349 (1980).
  * 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 .
  * Eichten and Hill (1990) E. Eichten and B. R. Hill, Phys. Lett. B 234, 511 (1990).
  * Lunghi  _et al._ (2003) E. Lunghi, D. Pirjol, and D. Wyler, Nucl. Phys. B 649, 349 (2003), arXiv:hep-ph/0210091 .
  * 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 .
  * Falk  _et al._ (1990) A. F. Falk, H. Georgi, B. Grinstein, and M. B. Wise, Nucl. Phys. B 343, 1 (1990).
  * Beneke  _et al._ (2020) M. Beneke, C. Bobeth, and Y.-M. Wang, JHEP 12, 148 (2020), arXiv:2008.12494 [hep-ph] .
  * Wang and Shen (2018) Y.-M. Wang and Y.-L. Shen, JHEP 05, 184 (2018), arXiv:1803.06667 [hep-ph] .
  * Braun and Khodjamirian (2013) V. M. Braun and A. Khodjamirian, Phys. Lett. B 718, 1014 (2013), arXiv:1210.4453 [hep-ph] .
  * Wang (2016) Y.-M. Wang, JHEP 09, 159 (2016), arXiv:1606.03080 [hep-ph] .
  * Grozin and Neubert (1997) A. G. Grozin and M. Neubert, Phys. Rev. D 55, 272 (1997), arXiv:hep-ph/9607366 .
  * Braun  _et al._ (2004) V. M. Braun, D. Y. Ivanov, and G. P. Korchemsky, Phys. Rev. D 69, 034014 (2004), arXiv:hep-ph/0309330 .
  * Khodjamirian  _et al._ (2007) A. Khodjamirian, T. Mannel, and N. Offen, Phys. Rev. D 75, 054013 (2007), arXiv:hep-ph/0611193 .
  * 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] .
  * Khodjamirian  _et al._ (2020) A. Khodjamirian, R. Mandal, and T. Mannel, JHEP 10, 043 (2020), arXiv:2008.03935 [hep-ph] .
  * Beneke  _et al._ (2018a) M. Beneke, V. M. Braun, Y. Ji, and Y.-B. Wei, JHEP 07, 154 (2018a), arXiv:1804.04962 [hep-ph] .
  * Beneke  _et al._ (2018b) M. Beneke, C. Bobeth, and R. Szafron, Phys. Rev. Lett. 120, 011801 (2018b), arXiv:1708.09152 [hep-ph] .
  * Cerri  _et al._ (2019) A. Cerri _et al._ , CERN Yellow Rep. Monogr. 7, 867 (2019), arXiv:1812.07638 [hep-ph] .
  * Ciuchini  _et al._ (2021b) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, (2021b), arXiv:2110.10126 [hep-ph] .
  * Altmannshofer and Archilli (2022) W. Altmannshofer and F. Archilli, in  _2022 Snowmass Summer Study_ (2022) arXiv:2206.11331 [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] .